Boundary conditions for a Navier-Stokes simulation

517 views
Skip to first unread message

Martin Lesueur

unread,
Oct 9, 2015, 3:18:25 AM10/9/15
to moose-users
Hi everyone,
I'm doing simulations of Navier-Stokes in porous media (in a granular rock matrix).
Attached is one of my working simulation.
My goal is to use the result of the simulation to calculate the permeability of the sample using Darcy Law.
This imposes my boundary condition of the fluid flow to be a pressure gradient instead of the velocity.

_ The flow is going from left to right.
_ I have a no-slip condition on the grains and on the top and bottom boundaries.

However the system defined as it is does not give me the right result.
As you can see on the failed simulation (where vel_x is visualised), the fluid does not want to move forward, it keeps coming back left.

I managed to make the simulation work but with an unrealistic boundary condition. I set vel_x=0 on the left, from where my fluid is supposed to flow ?!

I don't understand why this simulation is working and the other one isn't. One of the thing I was wondering about was the need of an inlet of velocity if we set an inlet of pressure, because of the coupling of the variables ?

Any help on the matter would be appreciated. I really can't be satisfied with a working simulation if it means to set absurd boundary conditions.

Cheers,

Martin Lesueur
working simulation.png
failed simulation.png

Jean Francois Leon

unread,
Oct 9, 2015, 8:40:21 AM10/9/15
to moose-users
Which BC do you use on the right? ( and on the left)?

The bc you want to use most probably is a Dirichlet BC on the pressure with different values to create pressure difference and flow ( if they are equal ==>no flow)

JF

Peterson, JW

unread,
Oct 9, 2015, 10:36:29 AM10/9/15
to moose-users
On Fri, Oct 9, 2015 at 6:40 AM, Jean Francois Leon <j...@galtenco.com> wrote:
Which BC do you use on the right? ( and on the left)?

The bc you want to use most probably is a Dirichlet BC on the pressure with different values to create pressure difference and flow ( if they are equal ==>no flow)

If we are talking about incompressible Navier-Stokes, a Dirichlet BC for pressure only makes sense if you are solving a pressure-Poisson equation (PPE) form of the governing equations, otherwise there are no boundary conditions on the pressure.

There's a good discussion of the various formulations in Gresho & Sani's 1987 paper, http://tinyurl.com/pyeysdn.  Basically, they advocate solving the non-PPE form of the governing equations, but if you do use a PPE, to also use a "consistent" Neumann BC as that gives the best satisfaction of the divergence-free constraint and is consistent with the related Dirichlet BCs on the velocity which they assume in the paper.

--
John

mlesueur

unread,
Oct 12, 2015, 2:43:39 AM10/12/15
to moose...@googlegroups.com
Thank you John,
It still does not tell me why my simulation is working with this weird
boundary condition.
I went through some papers on PPE after reading your answer, but I don't
think I want to change my equations.
I just need to calculate the permeability with Darcy Law, so actually my
input can be a pressure gradient or the fluid velocity, I will still be
able to calculate the permeability.
Cheers,

Martin Lesueur

Jean Francois Leon

unread,
Oct 12, 2015, 7:34:19 AM10/12/15
to moose-users
Martin

Regarding to your specific question and as mentioned previously:

 what is your BC on the right boundary of your domain?
It is impossible to answer your specific question without knowing this.

Regarding your calculation . it is not difficult to achieve. I have done it in the past with another software using either pressure differences BC as suggested in my previous post and yes you need to write the equations properly for that, an inflow + outlet pressure condition  ( which can be zero as your fluid is incompressible) or an inflow + outflow. It should not matter if you setup the rest of the calculation properly:

The real difficulty with these calculations arise if you try to resolve in a non laminar flow regime. It should not concern you as it will fall well outside of the scope of "Darcy model".

In any event I strongly recommend you setup proper BC and then  test your software first with the Stoke equation ( kill the NL tern of NS) and then use very small initial value. NS real difficulty lies in the treatment of the nonlinear term.
JF

mlesueur

unread,
Oct 12, 2015, 8:32:44 AM10/12/15
to moose...@googlegroups.com
Hi JF,
I will try tomorrow to see if solving only Stoke equation will change something and keep you updated.

My BC on the right boundary of the domain is a Dirichlet of 0 for the pressure.
(the one on the left is a Dirichlet of 1 so that I have a gradient of pressure)

As I said, the simulation is not giving a coherent result unless I put a Dirichlet of 0 for vel_x on the left.
And this "working" simulation is giving me a pressure and velocity field that are (almost) the same as when I put boundary conditions of velocity instead or pressure (a velocity Dirichlet on the left)

Cheers,

Martin
--
You received this message because you are subscribed to a topic in the Google Groups "moose-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/moose-users/yWi8bdfZWvU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to moose-users...@googlegroups.com.
Visit this group at http://groups.google.com/group/moose-users.
To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/cf7c684d-a49d-4df6-8b11-50b335127557%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Jean Francois Leon

unread,
Oct 12, 2015, 10:32:40 AM10/12/15
to moose-users
Hi Martin

Just keep in mind that the choice of BC is not independent of the specific formulation of the equation ( what is your primary variable P or u?) and how do you write the coupling.
this is what JW was mentioning I think with the PPE formulation. Every formulation has its upside and downside and must be chose in the specific context you want to solve them


Keep in mind as it is also a possible cause of the problem:

There are constraints on the velocity field distribution at inlet and outlet with a noslip boundary condition on wall "parallel" to the flow ( in a pure laminar flow it is classical that the velocity profile is parabolic maximimum at center at a pipe outlet).
So trying to impose non physical  flow distribution profile as BC might lead to the problem you have as well. 
When I look at the pic you posted it seems that you force a constant normal velocity through the inlet. This is impossible as unphysical and ill posed as the corner of your domain will have 2 values 0 from the noslip and non zero from the BC value you set.

Again I have never used moose for NS but I assume that  if its not already there it should not be difficult to add a mass flow rate BC ( or average velocity as it is the same) that will calculate the proper velocity profile. You will use it as inlet and outlet and compute the pressure difference with an auxiliary variable or a postprocessor depending upon the way the equations are arranged in moose..
Other open source package have it and it should be easy to add to moose once you have determined the proper formulation.

Good luck
JF

mlesueur

unread,
Oct 13, 2015, 9:46:37 PM10/13/15
to moose...@googlegroups.com
Hi JF,

My question is indeed linked to the way the system of coupled variables are being solved. I am using MOOSE’s Navier-Stokes module, in which the pressure and velocity are solved in a coupled manner simultaneously. That’s why I’m posting my question to this mailing list, hoping that another MOOSE user would share their experience with that module.

My constraint is only on vel_x, I don't know have any problem on the corner because the velocity is fully vertical.

As you say, imposing a constant velocity profile on the boundary would give a non-physical flow profile. This is also why I want a boundary condition in pressure and not velocity.
But If I only put a boundary condition on pressure, the result is wrong, fluid does not want to flow forward.
So it looks that if I want to put a boundary condition on pressure, I have to put a boundary condition on velocity as well. Is it true ?
If this is the case, I tried to put 0 on the inlet, which gives me a satisfying behaviour past the boundary, but does not seem right to me.

Cheers,

Martin

Alexander Lindsay

unread,
Oct 14, 2015, 7:12:10 AM10/14/15
to moose...@googlegroups.com
Hi Martin,

Someone else can step in here and say otherwise, but to my knowledge the NS module is not being actively developed or used by anyone on the list, so you're probably not going to get too many replies from people with experience. 

It would be great if someone took up the CFD mantle for Moose. I know we had someone else on the list a few months ago looking to use Moose for CFD, but I think he was a little discouraged by the lack of users. I think we just need to develop a critical mass. 

Alex
You received this message because you are subscribed to the Google Groups "moose-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to moose-users...@googlegroups.com.

Peterson, JW

unread,
Oct 14, 2015, 9:54:16 AM10/14/15
to moose-users
On Tue, Oct 13, 2015 at 7:46 PM, mlesueur <martin....@ecl2015.ec-lyon.fr> wrote:
Hi JF,

My question is indeed linked to the way the system of coupled variables are being solved. I am using MOOSE’s Navier-Stokes module, in which the pressure and velocity are solved in a coupled manner simultaneously. That’s why I’m posting my question to this mailing list, hoping that another MOOSE user would share their experience with that module.

And you are using the INSMomentum and INSMass kernels from that module, is that correct?  If yes, then as I mentioned before, there should not be any boundary conditions on the pressure.

Generally speaking, you need to provide Dirichlet boundary conditions for velocity on the *entire* boundary except for the outlet, where you would use e.g. Griffiths' "no bc bc" [0].  Jean Francois also mentioned a "mass flow" BC.  I am not familiar with that one or it's application, but it may also be appropriate here so it would be great if he could elaborate on that...  On solid walls, the velocity should be Dirichlet zero, on the inlet, a parabolic profile (which is appropriate for Poiseuille, aka pressure-driven, flow) should be imposed as a Dirichlet condition.
 

My constraint is only on vel_x, I don't know have any problem on the corner because the velocity is fully vertical.

As you say, imposing a constant velocity profile on the boundary would give a non-physical flow profile. This is also why I want a boundary condition in pressure and not velocity.

Jean Francois's non-physical comment was regarding the application of non-zero velocity to a solid wall in the corners where it meets the inlet; this problem is avoided simply by using a parabolic inlet profile.
 


--
John

Peterson, JW

unread,
Oct 14, 2015, 10:09:11 AM10/14/15
to moose-users
On Wed, Oct 14, 2015 at 5:12 AM, Alexander Lindsay <adli...@ncsu.edu> wrote:
Hi Martin,

Someone else can step in here and say otherwise, but to my knowledge the NS module is not being actively developed or used by anyone on the list, so you're probably not going to get too many replies from people with experience. 

It would be great if someone took up the CFD mantle for Moose. I know we had someone else on the list a few months ago looking to use Moose for CFD, but I think he was a little discouraged by the lack of users. I think we just need to develop a critical mass. 

One issue is that "CFD" is too broad a term.  Are you referring to compressible or incompressible Navier-Stokes, or something in between?  Turbulent or laminar flows?  Enclosed or external flows?  Are there real gas effects or chemical reactions? Continuous or discontinuous formulations, or just finite volume methods?  What kind of stabilization does your problem require?  If the flow is compressible, are there shocks?  If so, how do you want to stabilize those?

One cannot use the same set of kernels and BCs to handle all these types of applications, and therefore creating a multipurpose module that handles all of them is quite challenging and time-consuming.

--
John

Jean Francois Leon

unread,
Oct 14, 2015, 11:27:24 AM10/14/15
to moose-users
JW

Let me elaborate on the "mass flow rate" integrated BC for incompressible NS but it is rather classical and boring..

Short answer:
if you have an inlet it is sufficient, in orcder to solve NS ( at least incompressible I dont know about compressible) to define Integral(u.n) over the inlet which is the mass flow rate for a density of 1. You don't need to define a point value.

Longer now:

In order to demonstrate the statement above, ones  need to manipulate weak form and I wont do it as it will be too lengthy and I don't even know how to use the markdown language for equations on the web ..
 but it is CLASSICAL in NS theory world.
Any text book that details solving NS will talk about it. ( under a number of assumptions) 

Zero tangential velocity through the inlet area is the most important one.

Why is that?
 heuristically assume that an infinite long pipe bring the fluid to your computational domain with constant shape and zero tangential velocity. calculating the velocity profile in this pipe is possible from NS assuming you know the velocity at the cross section boundary ( 0 with a no-slip BC).  It give you the classical parabolic profile if the shape of your pipe is simple enough but the calculation can be carried away numerically for complex cross section shape as well. there is one free constant in the solution and you usually fix  it by defining the mass flow.

Stated otherwise if you have a cross section with only normal velocity and known velocity at the boundary of the cross section you can write from NS a 2d pde that solve the profile distribution in the cross section assuming an infinite length pipe of constant cross section.
This is the same equation the weak form will solve. This is the generalization of the calculation you will find in textbook where the parabolic profile is calculated.

On a practical note I know of 2 codes that I have used over the years that implement this BC
dumux- open source (dumux.org)
and comsol ( commercial)
If you have access to one of this tool you can have an idea on how this work in details.
Hope this help
JF

Jean Francois Leon

unread,
Oct 14, 2015, 11:44:07 AM10/14/15
to moose-users
Forgot to mention in previous post that the 2d equation that solve the normal velocity profile through a pipe cross section is just a 2d laplace equation (or a heat equation for transient problem) as the non linear term is zero under these assumptions.
 
JF

Martin Lesueur

unread,
Mar 13, 2016, 10:19:25 PM3/13/16
to moose-users
Hi,
I'm reopening my old thread because I have some updates.

As a reminder, my problem was to run a simulation of NS with a pressure Dirichlet BC instead of a velocity one.
It is at the moment not working with the current module in MOOSE.

I managed to fix the problem by changing the weak formulation of the INSMomentum Kernel (https://github.com/idaholab/moose/blob/devel/modules/navier_stokes/src/kernels/INSMomentum.C).
This Kernel is responsible for the appearance of gradP in Navier-Stokes equation:
\frac{\partial \left( \rho \vec{v} \right)}{\partial t} + \vec{\nabla} \cdot \left(\rho \vec{v} \otimes \vec{v} \right) = - \vec{\nabla} p + \vec{\nabla} \cdot \overline{\overline {\tau}} + \rho \vec{f}

In the kernel, it is formulated as:
    Real pressure_part = -_p[_qp] * _grad_test[_i][_qp](_component);
We are using the gradient of the test function because we applied the divergence theorem on gradP.
Attached you can see how the simulation fails with this formulation.

Now what I am using is:
    Real pressure_part = _grad_p[_qp](_component) * _test[_i][_qp];
(_grad_p being the coupledGradient of p)
That means that I haven't applied the divergence theorem.
Attached is my working simulation using this new formulation.

The two formulations are both correct, only that the first one introduces in addition a natural boundary condition on gradP. So I'm wondering why this change of formulation makes my simulation working.
I was also wondering if the current MOOSE formulation had been written on purpose (stability reasons, etc).
I would do a pull request otherwise.

Cheers,

Martin
old_formulation.png
new_formulation.png

Peterson, JW

unread,
Mar 14, 2016, 10:34:21 AM3/14/16
to moose-users
On Sun, Mar 13, 2016 at 8:19 PM, Martin Lesueur <martin....@ecl2015.ec-lyon.fr> wrote:
Hi,
I'm reopening my old thread because I have some updates.

As a reminder, my problem was to run a simulation of NS with a pressure Dirichlet BC instead of a velocity one.
It is at the moment not working with the current module in MOOSE.

I managed to fix the problem by changing the weak formulation of the INSMomentum Kernel (https://github.com/idaholab/moose/blob/devel/modules/navier_stokes/src/kernels/INSMomentum.C).
This Kernel is responsible for the appearance of gradP in Navier-Stokes equation:
\frac{\partial \left( \rho \vec{v} \right)}{\partial t} + \vec{\nabla} \cdot \left(\rho \vec{v} \otimes \vec{v} \right) = - \vec{\nabla} p + \vec{\nabla} \cdot \overline{\overline {\tau}} + \rho \vec{f}

In the kernel, it is formulated as:
    Real pressure_part = -_p[_qp] * _grad_test[_i][_qp](_component);
We are using the gradient of the test function because we applied the divergence theorem on gradP.
Attached you can see how the simulation fails with this formulation.

Now what I am using is:
    Real pressure_part = _grad_p[_qp](_component) * _test[_i][_qp];
(_grad_p being the coupledGradient of p)
That means that I haven't applied the divergence theorem.
Attached is my working simulation using this new formulation.

The two formulations are both correct, only that the first one introduces in addition a natural boundary condition on gradP. So I'm wondering why this change of formulation makes my simulation working.

Hmm... for an outflow boundary I'm not sure that dp/dn=0 (the natural BC) is necessarily correct.  It should be correct on solid walls, though.

Perhaps the right approach is to actually compute the contribution of this boundary term and add it to the residual... this is similar in spirit to the "no-bc-bc" boundary condition of Griffiths linked earlier in this thread.

By the way, I still maintain that Dirichlet BCs for pressure don't make sense in any case.  I refer you again to the Gresho & Sani paper linked earlier in the thread.
 


I was also wondering if the current MOOSE formulation had been written on purpose (stability reasons, etc).

The reason that this term is normally integrated by parts in the momentum equation is so that the "B" and "B^T" parts of the resulting Jacobian:

[ A   B ]
[ B^T 0 ]

are actually symmetric.  This makes it possible (in the Stokes flow case) to take advantage of iterative solvers which are appropriate for symmetric problems, and in Schur complement solvers to avoid storing both B and B^T, but is of less importance in the current MOOSE implementation.


 
I would do a pull request otherwise.

I agree that there may be a scenario in which using the non-integrated-by-parts version of the pressure gradient term would be useful/necessary.  In fact, I have experimented with this for the compressible equations, see NSMomentumInviscidFluxWithGradP.h.

So, if you would like to add a boolean flag "integrate_p_by_parts" to the INSMomentum Kernel which defaults to true, and then compute the pressure term differently depending on that flag, that would be a useful PR.

--
John

Roma Gurung

unread,
May 9, 2016, 10:14:18 AM5/9/16
to moose-users
Hi All,
The associated update in INSMomentum.C is presented by J.W. Peterson as
https://github.com/idaholab/moose/commit/1922ccd184fc3c4f702697fe07e365907db706db

Yours Sincerely,
Anil Kunwar
Reply all
Reply to author
Forward
0 new messages