CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

Adding implicit source term to momentum equation

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By fs82

Reply
 
LinkBack Thread Tools Display Modes
Old   September 18, 2009, 04:31
Default Adding implicit source term to momentum equation
  #1
Senior Member
 
Fabian Schlegel
Join Date: Apr 2009
Location: Dresden, Germany
Posts: 182
Rep Power: 8
fs82 is on a distinguished road
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
fs82 is offline   Reply With Quote

Old   September 21, 2009, 04:32
Default
  #2
Senior Member
 
dmoroian's Avatar
 
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11
dmoroian is on a distinguished road
Hello Fabian,
Did you check this post:sourceTerm ?

I hope this is helpful,
Dragos
dmoroian is offline   Reply With Quote

Old   September 21, 2009, 10:19
Default
  #3
Senior Member
 
Fabian Schlegel
Join Date: Apr 2009
Location: Dresden, Germany
Posts: 182
Rep Power: 8
fs82 is on a distinguished road
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
fs82 is offline   Reply With Quote

Old   September 21, 2009, 12:12
Default
  #4
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
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
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   September 22, 2009, 06:07
Default
  #5
Senior Member
 
Fabian Schlegel
Join Date: Apr 2009
Location: Dresden, Germany
Posts: 182
Rep Power: 8
fs82 is on a distinguished road
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
Tushar@cfd likes this.

Last edited by fs82; September 22, 2009 at 07:14.
fs82 is offline   Reply With Quote

Old   September 22, 2009, 14:15
Default
  #6
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
Quote:
Originally Posted by fs82 View Post
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,
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   September 23, 2009, 03:29
Default
  #7
Senior Member
 
Fabian Schlegel
Join Date: Apr 2009
Location: Dresden, Germany
Posts: 182
Rep Power: 8
fs82 is on a distinguished road
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
fs82 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
momentum source term zwdi FLUENT 13 December 5, 2013 18:58
Mass source, Momentum source theory. diffo FLUENT 4 August 21, 2009 10:26
Problem with Joulebs effect source term in the energy equation galaad OpenFOAM Running, Solving & CFD 0 January 19, 2006 13:01
PDE in udf for defining momentum source term shekhar FLUENT 0 April 12, 2004 01:43
adding source term confusion tedmcc CFX 4 September 27, 2001 05:05


All times are GMT -4. The time now is 10:30.