[alfiogn]

Hi,

I'm trying to implement an fvOption to add a new source to my solver without implementing a new solver.

In particular, what I want to do is applying a Brinkman penalization term to impose velocity or temperature in a certain cellSet like a Dirichlet condition. (I want to do it also for a Neumann condition but it is a further step).
It's something like an immersed boundary. I've to implement it for my phd studies.

Analytically the term is represented by:

where u is my unkown variable, is the viariable field that I want to impose and is the penalization coefficient (10^-6 ~ 10^-12). I also want to treat u implicitly.

I've implemented an fvOption template for scalar and vector variables. The key part is the following:

template<class Type>
void Foam::fv::PenalizationSource<Type>::addSup
(
    fvMatrix<Type>& eqn,
    const label fieldi
)
{
    const volFieldType& psi = eqn.psi();
    eqn += (fvm::Sp(Chi_/epsilon_, psi)
            - fvc::Sp(Chi_/epsilon_, Ubar_));
}

Chi_ is a characteristic function, Ubar_ is the imposed value.
volFieldType is a typedef for a GeometricField<Type>.


What happens during simulations is that the velocity field explodes but I cannot imagine how to fix this problem. Maybe my code is wrong.
I think that the problem is in the summation of the fvm and fvc terms to my equation but I didn't figured out why.

Can someone help me?

Thanks in advance
Giorgio


UPDATE:
Maybe I found out the issue: the fvc term is added well, but the fvm one is not added to the equation. In this way I have an unbalanced source term (10^10 in rhs and the normal matrix in lhs).
However I don't know how to fix it.

[alfiogn]

I feel stupid but I don't understand the reason why now the code is working perfectly. I just changed the sign of the fvm and fvc term.

Can someone explain me what is the difference between now and before?

template<class Type>
void Foam::fv::PenalizationSource<Type>::addSup
(
    fvMatrix<Type>& eqn,
    const label fieldi
)
{
    const volFieldType& psi = eqn.psi();
    eqn -= (fvm::Sp(Chi_/epsilon_, psi)
                - fvc::Sp(Chi_/epsilon_, Ubar_));

}

And it is NOT working if I write

template<class Type>
void Foam::fv::PenalizationSource<Type>::addSup
(
    fvMatrix<Type>& eqn,
    const label fieldi
)
{
    const volFieldType& psi = eqn.psi();
    eqn -= (-fvm::Sp(Chi_/epsilon_, psi)
                + fvc::Sp(Chi_/epsilon_, Ubar_));

}