[fs82]

after spending days in unterstanding the source code and the mathematics of the oodles solver to add a velocity dependend source term to the navier stokes equation I now have a idea but I am not really sure if this is right.
So hopefully somebody could verfify my idea and here it comes:

I want to add the following term to my navier stokes equation on the right hand side: F_i = - cd * a * mag(U) * U
I want to treat this term implicit. After searching the forum and the source code of the other solver it should be possible to do this:

fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ sgsModel->divDevBeff(U)
+ fvm::Sp(cd*a*mag(U),U)
);
For my understanding OpenFoam will handle U implizit but mag(U) will be explizit? But I is it also possible and reasonable to do something like this:

fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ sgsModel->divDevBeff(U)
+ fvm::Sp(mag(fmv::Sp(cd*a, U)), U)
);

cd is a scalar and a in the first approach too but in the end it should be a scalar field. So I havent tried this yet it is still theoretical but may be someone could verify my idea :-D

thx

Fabian

[dmoroian]

Hello Fabian,
Did you check this post:sourceTerm ?

I hope this is helpful,
Dragos

[fs82]

Allright,

thx for your hint, but this is not what I'm looking for. But I spent some more time on the theory and tried the second idea in my LES Model. So as a result:

fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ sgsModel->divDevBeff(U)
+ fvm::Sp(mag(fmv::Sp(cd*a, U)), U)
);

is not resonable and do not work!
But the first approach is also not correct. This works but there is something missing.
If I apply a taylor series to liniarise my source term: F(U) = mag(U) * U
I get: F(U) = F(U°) + dF/dU * (U° - U) where ° indicates stands for explizit.
If one calculates the derivative dF/dU, one obtain:
F(U) = mag(U°)*U - U°*magSqr(U°) + U° * ( U° & U ) where & stands for the inner product (like OpenFoam)
The first term is the approach i allready mentioned but the second and third term a missing in my first idea. But I am not able to implement this liniarisation of my source term, because I am not able to access the components of my implizit treated U vector. Also I am not allowed to do a inner product with a fvVectorMatrix result from fvm::Sp(1., U). May be anybody has an idea?

thx,
Fabian

[alberto]

You are dealing with something similar to the addition of the implicit part of a drag term in multiphase flow, which looks like -K(Ur)*Ua, where Ur = Ua - Ub is the relative velocity and U is the velocity you are solving for.

This is done usually by assuming Ur lagged: in other words, it is computed with the previous iteration or time step value of Ua. Of course you might need internal iterations to obtain a converged solution (switch to unsteady SIMPLE algorithm for example).

Best,
Alberto

[fs82]

Dear Alberto,

I guess switching to SIMPLE Algorithm is not the best choice. Because SIMPLE is the same procedure like PISO but less stable. And my aim is to improve the stability when I try to implement my drag term as implicit as possible. So improving stability with an implicit drag term and reducing it by choosing SIMPLE instead of PISO would produce no benefit.
So my solution is at the moment to test my second approach:
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ sgsModel->divDevBeff(U)
+ fvm::Sp(mag(fmv::Sp(cd*a, U)), U)
);
if it is stable enough for my applications and to neglect the last two terms of the liniarisation. I think the problem is not the PISO solution technique but the solver for the equation system, which couldn't handle terms which mix the different components for one equation. The same problem should exist if you look at the diffusion term of the oodles solver:

- fvm::laplacian(nuEff(), U) - fvc::div(nuEff()*dev(fvc::grad(U)().T()))

Due to the local varying nuEff the second term isn't vanishing and you get the same problem. Parts from the other velocity components have to be taken into account for one velocity component and this is treated explicit too. I guess this is a problem of the solver for the equation system. But it is just an idea :-D

thx for your help
Fabian

[alberto]

Originally Posted by fs82

Dear Alberto,

I guess switching to SIMPLE Algorithm is not the best choice. Because SIMPLE is the same procedure like PISO but less stable. And my aim is to improve the stability when I try to implement my drag term as implicit as possible. So improving stability with an implicit drag term and reducing it by choosing SIMPLE instead of PISO would produce no benefit.

So my solution is at the moment to test my second approach:
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ sgsModel->divDevBeff(U)
+ fvm::Sp(mag(fmv::Sp(cd*a, U)), U)
);
if it is stable enough for my applications and to neglect the last two terms of the liniarisation. I think the problem is not the PISO solution technique but the solver for the equation system, which couldn't handle terms which mix the different components for one equation.

I actually did not say the problem is the PISO algorithm, but that you could linearize the source term lagging one of U terms, and iterating to obtain a converged solution.

I disagree however with your statement about the stability of a pure PISO algorithm and a SIMPLE algorithm. Unsteady SIMPLE (and corrected flavours like SIMPLEC, for example) with under-relaxation are often more efficient when the solution of the equations is problematic since you are not tied to the very small time steps a PISO algorithm requires, and this without losing in stability.

One last note. Is this syntax correct?

fvm::Sp(mag(fmv::Sp(cd*a, U)), U)

I mean, methods in fvm::Sp gives you the contribution to the matrix of a given term. I would chech that mag(fvm::Sp(cd*a,U) is actually what you want.

Best,

[fs82]

Hello Alberto,

thx for your answer. I read the book of Ferzinger and Peric about numerical fluid mechanics (I only know the german title: "Numerische Strömungsmechanik") and there is a brief discription of the so called "projection methods". In the pressure equation there will be a term which includes the velocity corrections and this term is neglected for SIMPLE. And Ferzinger also states that SIMPLE do not converge as good as improved variations of SIMPLE, like PISO or SIMPLEC. May be the under-relaxation technique improves SIMPLE, but my knowlege is not so well about this topic.
My force term looks like this: F_i = -cd * a * abs(U) * U
To handle it implicit you have to apply a taylor series for liniarising the source term:
F(U) = F(U°) + (dF/dU)° * (U-U°)
with dF_j/dU_k = U_k * U_j / sqrt(U_i*U_i) + sqrt(U_i*U_i) * kronecker_delta_ij
This leads to:
F_j = abs(U°)*U_j + sum_k [ 1/abs(U°) * U°_j * U°_k * U'_k ] with
U'_k = U_k-U°_k
The last term I am not able to get him to work in OpenFoam, so I negect him. My F_j looks:
F_j = cd*a*abs(U°)*U_j
cd is a drag coefficient and a is a function of height and has the dimension 1/m.
So my code looks like this: fvm::Sp(mag(fvc::Sp(cd_*lad_,U)),U) ...
Oh sry I see, I made a small mistake ... the second fvm::Sp I have written above should be a fvc::Sp.

thx for your help,

Fabian