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

Poisson equation,Neumann BCs

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

Like Tree3Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   March 25, 2014, 17:10
Default
  #21
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,588
Rep Power: 20
FMDenaro will become famous soon enough
Quote:
Originally Posted by adrin View Post
It turns out that FMDenaro is correct. So long as the Neumann BC compatibility is satisfied _explicitly_ one need not set a Dirichlet BC at one node, and one can still get a solution. In my previous experiments in the finite volume formulation I was using a simple manufactured problem for benchmarking. It turns out that although compatibility was theoretically ensured it was _not_ satisfied numerically (in discrete form). So, a simple correction led to "a" solution. That solution shifts by a constant as a function of grid resolution, but that's not an issue. I still maintain the only reason any solution may be obtained with such an approach is sufficient numerical noise to remove the linear dependence of two arbitrary equations - this approach would fail in an "infinite accuracy" solution.

So, in summary, if the _discrete_ surface integral of the fluxes is not equal to the volume integral of the Poisson source term, perform the following pre-processing before assembling the matrix and RHS. Find the difference between the aforementioned two terms, and either (a) divide that difference by the total volume and then add the volume averaged error to the source term for every node, or (b) divide the difference by the total surface area and then add the surface averaged error to the fluxes on all surface nodes. Make sure the signs are such that compatibility is now satisfied. This yields a quick, converged solution!

Adrin
just to say that while using a FV-based flux discretization I checked in my codes that it automatically ensures the compatibility constraint (in discrete sense) without any correction..
FMDenaro is offline   Reply With Quote

Old   March 25, 2014, 18:00
Default
  #22
Member
 
adrin
Join Date: Mar 2009
Posts: 41
Rep Power: 8
adrin is on a distinguished road
The compatibility may already be incorporated in the implementation, depending on the details of the formulation. I'm just curious ... suppose you have a unit cubic box with N^3 uniformly distributed grids. The volume for this will be one, irrespective of the value of N. Now, let's perturb the interior grids (keep the surface grids uniformly distributed), say (1) randomly, or (2) sinusoidally in all three directions (A*sin(2.pi.x)*sin(2.pi.y)*sin(2.pi.z)). What is the volume you get for such grid distributions?

Adrin
adrin is offline   Reply With Quote

Old   March 25, 2014, 18:15
Default
  #23
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,588
Rep Power: 20
FMDenaro will become famous soon enough
Quote:
Originally Posted by adrin View Post
The compatibility may already be incorporated in the implementation, depending on the details of the formulation. I'm just curious ... suppose you have a unit cubic box with N^3 uniformly distributed grids. The volume for this will be one, irrespective of the value of N. Now, let's perturb the interior grids (keep the surface grids uniformly distributed), say (1) randomly, or (2) sinusoidally in all three directions (A*sin(2.pi.x)*sin(2.pi.y)*sin(2.pi.z)). What is the volume you get for such grid distributions?

Adrin

to tell the truth I don't understand that question clearly ... the key is the numerical construction of the surface integral of dp/dn that must balance the volume integral of the source term whatever the shape of the domain and its discretization is used
Now, in the actual construction of the pressure equation

Div Grad p = q

the source term q is in the form of a divergence, therefore it is numerically constructed as the same as the divergence of the pressure gradient in the left hand side. That ensures that the compatibility constraint is ensured, provided that the correct BC for Grad p are used.
FMDenaro is offline   Reply With Quote

Old   March 25, 2014, 18:26
Default
  #24
Member
 
adrin
Join Date: Mar 2009
Posts: 41
Rep Power: 8
adrin is on a distinguished road
Ah, in your case q is in divergence form. In general, though, q could be any function, and in this case it's the volume integral of q that must match the surface integral of the flux. In discrete form the volumes may be off a bit, and so the volume integral of q may not match the LHS exactly.

Adrin
adrin is offline   Reply With Quote

Old   March 25, 2014, 18:34
Default
  #25
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,588
Rep Power: 20
FMDenaro will become famous soon enough
Quote:
Originally Posted by adrin View Post
Ah, in your case q is in divergence form. In general, though, q could be any function, and in this case it's the volume integral of q that must match the surface integral of the flux. In discrete form the volumes may be off a bit, and so the volume integral of q may not match the LHS exactly.

Adrin
Since the pressure equation is derived from the continuity equation Div V = 0 the source term is necessarily in the form of divergence ...

The question for a general function q I think that depends upon the discretization of the integrals ... for example, mid-point rule for the surface integral must be associated to the use of centroid value x measure of the volume... trapezoidal rule for the surface integral must be associated to its 3D extension for the volume integral and so on
FMDenaro is offline   Reply With Quote

Old   March 29, 2014, 09:12
Default
  #26
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 368
Rep Power: 10
arjun is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
Since the pressure equation is derived from the continuity equation Div V = 0 the source term is necessarily in the form of divergence ...
Not true.

/Arjun
.......................
arjun is offline   Reply With Quote

Old   March 29, 2014, 14:48
Default
  #27
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,588
Rep Power: 20
FMDenaro will become famous soon enough
Could You elaborate?
FMDenaro is offline   Reply With Quote

Old   March 30, 2014, 02:41
Default
  #28
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 368
Rep Power: 10
arjun is on a distinguished road
In theory, what you said is fine. In practice, these are few possible scenarios.

1. You meant div(v) = 0 on cell by cell basis. Which is only possiblibly be true in fully converged case, even then for a non - staggered system we would have div( v + v_due_to_rhie_chow) = 0 (fully converged case).

So I take it that you meant globally, that is summing all equations.

2. In this case, if you have say velocity inlet and a pressure outlet, until convergence is achieved div (vel ) = 0 is not achieved.
(this scenario could be discounted because then p equation is not all neumann anyway)

3. In a case if had inlet and outflow condition. In this case div(v) = 0 is only achieved if you fixed outflow such that inflow = outflow, that too when there is no reversed flow.

In the end, in practice, div(v) = 0 is not taken to be granted.
arjun is offline   Reply With Quote

Old   March 30, 2014, 05:29
Default
  #29
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,588
Rep Power: 20
FMDenaro will become famous soon enough
Quote:
Originally Posted by arjun View Post
In theory, what you said is fine. In practice, these are few possible scenarios.

1. You meant div(v) = 0 on cell by cell basis. Which is only possiblibly be true in fully converged case, even then for a non - staggered system we would have div( v + v_due_to_rhie_chow) = 0 (fully converged case).

So I take it that you meant globally, that is summing all equations.

2. In this case, if you have say velocity inlet and a pressure outlet, until convergence is achieved div (vel ) = 0 is not achieved.
(this scenario could be discounted because then p equation is not all neumann anyway)

3. In a case if had inlet and outflow condition. In this case div(v) = 0 is only achieved if you fixed outflow such that inflow = outflow, that too when there is no reversed flow.

In the end, in practice, div(v) = 0 is not taken to be granted.

ok, now I see what you meant. However, let me address some issues.

1) Also in practice, the Div V = 0 constraint can be reached at machine precision. It depend on your choice, if you want to adopt an Exact or Approximate Projection Method (EPM/APM).
2) R-C interpolation is not necessary, you are not forced to use it as other approaches are possible.
3) In a case with inlet and outflow condition Div V = 0 is achieved because, provide correct BC are fixed, you get inflow = outflow at machine accuracy

If you are interested in, I report two of my papers about the above issues, you will find cited many other useful papers.

F.M. Denaro, A 3D second-order accurate projection-based Finite Volume code on non-staggered, non-uniform structured grids with continuity preserving properties: application to buoyancy-driven flows,Int. J. Num. Meth. Fluids, 52, 4, 393-432, 2006.A. Aprovitola, F.M. Denaro, A non-diffusive, divergence-free, Finite Volume-based double projection method on non-staggered grids, Int. J. Num. Meth. Fluids,53, 1127-1172, 2007.
FMDenaro is offline   Reply With Quote

Old   March 30, 2014, 09:03
Default
  #30
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 368
Rep Power: 10
arjun is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post

1) Also in practice, the Div V = 0 constraint can be reached at machine precision. It depend on your choice, if you want to adopt an Exact or Approximate Projection Method (EPM/APM).
2) R-C interpolation is not necessary, you are not forced to use it as other approaches are possible.
3) In a case with inlet and outflow condition Div V = 0 is achieved because, provide correct BC are fixed, you get inflow = outflow at machine accuracy
1) That is only after (3) that is after you force inflow = outflow, provided that no reverse flow is occuring, because in this case you can not fix flow that is going out.

2) I agree that it is matter of practice.

PS: for (1) if you meant cell by cell basis then I would not agree because div(u) = 0 would only be achieved after convergence, but the pressure equation is driving the solution to convergence div(v) = 0 cell by cell basis is not given. In fact that is source of pressure equation, without source there will be no correction.
arjun is offline   Reply With Quote

Old   March 30, 2014, 09:10
Default
  #31
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,588
Rep Power: 20
FMDenaro will become famous soon enough
Quote:
Originally Posted by arjun View Post

PS: for (1) if you meant cell by cell basis then I would not agree because div(u) = 0 would only be achieved after convergence, but the pressure equation is driving the solution to convergence div(v) = 0 cell by cell basis is not given. In fact that is source of pressure equation, without source there will be no correction.

Sorry, I do not understand clearly this point ... at each time step I will solve the pressure equation (double precision) even with a simple SOR method. With the pressure gradient I can drive the divergence-free contraint cell-by-cell until to machine precision, even on non-staggered grid...
Any source term can be used to solve an elliptic equation and the solution provide the gradient which correct the velocity...
FMDenaro is offline   Reply With Quote

Old   April 4, 2014, 19:36
Default
  #32
New Member
 
Patricio Cumsille
Join Date: Jan 2014
Posts: 4
Rep Power: 3
pcumsill is on a distinguished road
Dear Adrin:

I have implemented the trick that you adviced me, but in practice the only thing I got is to project the RHS on the orthogonal of the kernel of the matrix obtained from the discretization Poisson problem with pure Neumann BC. Unfortunately, you still cannot ensure the uniqueness of the solution and therefore you do not know which solution is computed. The experiment that I have performed is with FD and the test function is u(x,y)=x^5+y^5-r^5, with r=0.5 and the RHS is obtained by computing the Laplacian of u, and so is for the Neumann BC.

The issue is that even if you project the RHS on the orthogonal of the kernel of the matrix for this problem, you cannot ensure the uniqueness of the solution and therefore you cannot ensure the convergence of the numerical method. After all, the analytic solution is defined up to an additive constant and any method that converges have to handle correctly this issue. The fact of projecting the RHS is not enough to handle that!

Best regards,
Patricio Cumsille
pcumsill is offline   Reply With Quote

Old   April 11, 2014, 21:53
Default
  #33
Member
 
adrin
Join Date: Mar 2009
Posts: 41
Rep Power: 8
adrin is on a distinguished road
Hi Patricio,

Sorry, I hadn't checked cfd-online for a while...

Yes, you are absolutely correct that the solution is/will not be unique by just discretely ensuring the consistency - the same problem at different resolution levels will give different potential distributions (but the solution is only shifted up or down; the shape is the same). But, as you suggested, that is to be expected because the original problem is itself not unique. The only way to do that properly would be to remove one of the equations (any one would do, which basically means you're setting its value to zero) and replace it with the consistency equation. That is essentially what I did in my boundary element methods based solution, where a unique solution is obtained irrespective of the resolution.

Nevertheless, since we're generally interested in the gradient and not the potential itself, it really doesn't matter if the value shifts up/down - the only real concern is whether we're guaranteed to get a solution. I suspect in an infinite precision computation we won't get a converged solution. So far, I haven't experienced a problem in my finite precision tests.

If you want a "pinned" solution, that means you know the potential value at a particular point (no other option available, by definition). In that case, just get the difference between the value that you have obtained and the value that you want at that node, and then adjust the value at all nodes by this difference. This will give you the unique values you want (but it's a simple post-processing and has nothing to do with your real concern of convergence guarantee)

Adrin
adrin is offline   Reply With Quote

Old   April 23, 2015, 17:32
Default imposing pressure gradient at free surface of a channel flow
  #34
New Member
 
marjan
Join Date: Aug 2013
Posts: 11
Rep Power: 3
mdallali is on a distinguished road
Dear All,
I am trying to simulate the free-surface shape, by imposing a vertical pressure gradient at uppest ghost cell of a rectangular channel.
I m using a LES code, which solves the Navier -Stockes Equation, and the flow is driven by a stream- wise pressure gradient while the periodicity condition is applied in stream -wise and span- wise direction.
To do that I modify the Poisson equation by changing the Boundary condition to desired pressure gradient which is do=- g*dh/dz.(dh /dz is presenting the variation of free-surface level along the span wise direction) at upper grid and a hydrostatic pressure gradient along all the vertical direction dp/dh=-g.

After Running my code, The residual of poisson equation is not converging and mass-balance get apart from zero.
does anyone knows anything related to this issue?

Thanks in advance
mdallali is offline   Reply With Quote

Old   April 23, 2015, 17:46
Default
  #35
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,588
Rep Power: 20
FMDenaro will become famous soon enough
Quote:
Originally Posted by mdallali View Post
Dear All,
I am trying to simulate the free-surface shape, by imposing a vertical pressure gradient at uppest ghost cell of a rectangular channel.
I m using a LES code, which solves the Navier -Stockes Equation, and the flow is driven by a stream- wise pressure gradient while the periodicity condition is applied in stream -wise and span- wise direction.
To do that I modify the Poisson equation by changing the Boundary condition to desired pressure gradient which is do=- g*dh/dz.(dh /dz is presenting the variation of free-surface level along the span wise direction) at upper grid and a hydrostatic pressure gradient along all the vertical direction dp/dh=-g.

After Running my code, The residual of poisson equation is not converging and mass-balance get apart from zero.
does anyone knows anything related to this issue?

Thanks in advance

could you better elaborate your model? I think you are not correctly setting the BC.s. You can check if the compatibility relation is numerically satisfied
FMDenaro is offline   Reply With Quote

Old   April 23, 2015, 18:03
Default
  #36
New Member
 
marjan
Join Date: Aug 2013
Posts: 11
Rep Power: 3
mdallali is on a distinguished road
Dear Fillipo ,
The mass balance is also getting far form zero, but I can not find where it rise up from.
mdallali is offline   Reply With Quote

Old   April 23, 2015, 18:43
Default
  #37
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,588
Rep Power: 20
FMDenaro will become famous soon enough
Quote:
Originally Posted by mdallali View Post
Dear Fillipo ,
The mass balance is also getting far form zero, but I can not find where it rise up from.
Without details is not possibile to help you... I suppose you are not correctly setting the pressure gradient at the boundary.
Are you using the correct normal component inserted in the Hodge decomposition?
FMDenaro is offline   Reply With Quote

Old   April 23, 2015, 22:46
Default
  #38
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 368
Rep Power: 10
arjun is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
Sorry, I do not understand clearly this point ... at each time step I will solve the pressure equation (double precision) even with a simple SOR method. With the pressure gradient I can drive the divergence-free contraint cell-by-cell until to machine precision, even on non-staggered grid...
Any source term can be used to solve an elliptic equation and the solution provide the gradient which correct the velocity...





Edited to add:

I went back and read again what you wrote. Sorry I misunderstood you. You were saying the same thing (div (V) = 0 is used to construct pressure equation rather than it is a given condition. Discard my comment above.
FMDenaro likes this.

Last edited by arjun; April 24, 2015 at 00:19.
arjun is offline   Reply With Quote

Reply

Tags
neumann bcs, poisson equation

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
Dealing with BC's in OF 1.6 vkrastev OpenFOAM Running, Solving & CFD 5 September 4, 2012 11:58
BCs for Pressure Correction Equation (SIMPLE) Bharath Somayaji Main CFD Forum 1 March 1, 2006 07:12
Recommendation of a good poisson solver Quarkz Main CFD Forum 2 December 2, 2005 10:12
Poisson Equation in CFD Maciej Matyka Main CFD Forum 9 November 10, 2004 12:30
Poisson BC's Peter Main CFD Forum 2 May 18, 2001 04:40


All times are GMT -4. The time now is 03:09.