CFD Online URL
[Sponsors]
Home > Forums > Main CFD Forum

help needed in solving Neumann peoblem

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

Reply
 
LinkBack Thread Tools Display Modes
Old   January 24, 2001, 09:49
Default help needed in solving Neumann peoblem
  #1
yf yap
Guest
 
Posts: n/a
in solving a channel flow problem, i need to solve a elliptic PED derived by introducing a velocity potential into the continuity equation(in boundary fitted coordinate), such that: U=d(phi)/dx and V=d(phi)/dy the continuity equation:

A d(U)/dx + B d(V)/dy = Source

A d^2(phi)/dx^2 + B d^2(phi)/dy^2 = Source ------------(1)

note: d( )/dx , d( )/dy , d^2( )/dx^2 , d^2( )/dy^2 are partial derivatives,

A and B are functions of x and y. boundary conditions: at the channel wall U=d(phi)/dx=0 and V=d(phi)/dy=0 (neumann boundary conditions) (1) is discretized using second order central difference, and for grid points at the wall, 3-point, one-sided second order. after substitutions are made, a linear system is obtained : [C][phi_i,j]=[Source_i,j] after ascertain that [C] is strongly diagonal dominant, it is solved using SOR. however, the solution was found diverging.

QUESTIONS 1)in handling Neumann problem, like (1), green's first integral theorem is not satisfied due to the truncation errors, can this be the problem?

2)what can be done in order to avoid divergence(SOR)?

regards, yfyap
  Reply With Quote

Old   January 24, 2001, 20:44
Default Re: help needed in solving Neumann peoblem
  #2
John C. Chien
Guest
 
Posts: n/a
(1). Because of the limited time and space available, it is not possible to get involved in your research activities. (2). Therefore, I would suggest that you first pick a method or a good technical paper, and try to follow the method presented. In this way, you are sure of the result you will be getting. (3). For new experiments and ideas, you are on your own. You sure can bring it to the technical meetings where you can discuss it in greater detail with other researchers.
  Reply With Quote

Old   January 24, 2001, 21:56
Default Re: help needed in solving Neumann peoblem
  #3
Peter Attar
Guest
 
Posts: n/a
Check the condition number of C. While I doubt it,your problem may be ill-posed for one reason or another(most like due to an error in the calculation C) and you may need to use some sort of SVD-regularization method to solve it.Other than that ask your research advisor for some help.
  Reply With Quote

Old   January 25, 2001, 07:26
Default Re: Neumann peoblem compensation scheme proposed
  #4
Dean Schrage
Guest
 
Posts: n/a
If this is similar to the pressure poissons equation in terms of coupling and response, then the incompatibility may be a problem but can be compensated for in the interation process. Note that this is only required if the entire boundary is a Neumann BC. If there is a Dirichlet BC, the iteration will naturally sink the remaining source flux to this BC. Also, note that I have applied this to the PPE for pressue with Neumann BC's however this equation is a bit different than your inhomengenous equation with functions A, B. I would need to research more to assess this last effect. Regardless these are some notes on the approach that I have used.

Recognizing that the boundary source must equal the source condition in the volume, you can compensate as follows. Each wall has a Neumann BC, the flux of this BC, integrated over all the wall area must be equal to the volume integral of the source, i.e. divergence theorem

Sum(GRAD(Phi) en * dA) = Sum(S dV)

where * signifies dot product and en is an inward facing normal into the domain. However, this equation might not be satisfied during the iteration process.

The approach is to rewrite:

G = Sum(GRAD(Phi)en * dA)- Sum(S dV)

and compute the residue G for the solution domain. Next, revise each source term S by an volume prorated amount for each CV:

S_i <= S_i - G dV_i/V

With this, the source term integrated over the total volume at the current iteration level should be consistent with the flux of phi issuing from the boundaries. You can then proceed to complete the iteration cycle.

Now it may be possible to let the Neumann condition float and not be satisified. I have found that with SOR you don't necessarily need to satisfy the Neumann condition. Example, box with defined inlet and outlet velocities. True, the bulk pressure will continue to free-integrate, ever increasing but in this case the gradients are the desired element and the CFD simulation proceeds nicely even though the pressure is climbing.

Hope this helps

best regards

Dean
  Reply With Quote

Old   January 25, 2001, 22:42
Default Re: help needed in solving Neumann peoblem
  #5
yf yap
Guest
 
Posts: n/a
dear John, thanks for your advice. can you recommend books on neumann problems? because as far as i have found, books normally give very simple account on neumann problems, unlike dirichlet problems. thanks. regards, yfyap
  Reply With Quote

Old   January 25, 2001, 22:46
Default Re: help needed in solving Neumann peoblem
  #6
yf yap
Guest
 
Posts: n/a
dear Peter Attar, thanks for the advice. what is SVD-regularization method ? i haven't came across this method before. thanks regards, yfyap
  Reply With Quote

Old   January 25, 2001, 22:47
Default Re: Neumann peoblem compensation scheme proposed
  #7
yf yap
Guest
 
Posts: n/a
thank you, Dean. i just don't know how am i going to express my gratitude for the precious advice that you gave me.

regarding the message that you have posted, i have difficulty understanding the last paragraph of it.

"Now it may be possible to let the Neumann condition float and not be satisified. I have found that with SOR you don't necessarily need to satisfy the Neumann condition. Example, box with defined inlet and outlet velocities. True, the bulk pressure will continue to free-integrate, ever increasing but in this case the gradients are the desired element and the CFD simulation proceeds nicely even though the pressure is climbing."

questions: 1) what does it means by "let the Neumann condition float" ? 2) what does the last sentence means , ie "the bulk pressure will continue to free-integrate, ever increasing but in this case the gradients are the desired element and the CFD simulation proceeds nicely even though the pressure is climbing" ? 3) i am trying to incorporate the residue, G, into the source, at the moment. can i email you in case there is difficulty?

thanks. regards, yfyap
  Reply With Quote

Old   January 26, 2001, 02:02
Default Re: Neumann peoblem compensation scheme proposed
  #8
yf yap
Guest
 
Posts: n/a
regarding to the elliptic continuity equation, (1), it is discovered that it can actually be transformed into a simpler form. recall that: after substitutions be made to the continuity equation by defining U=d(phi)/dx and V=d(phi)/dy, it is given by:

A d(U)/dx + B d(V)/dy = Source

A d^2(phi)/dx^2 + B d^2(phi)/dy^2 = Source ------------(1)

note: d( )/dx , d( )/dy , d^2( )/dx^2 , d^2( )/dy^2 are partial derivatives,

A and B are functions of x and y. however, (1) can be reduced into a simpler form after an appropiate substitution be made, that is:

d^2(phi)/dx^2 + D d(phi)/dx + d^2(phi)/dy^2 = Source ----(2) which D=D(x)

the same boundary conditions apply to (2), ie. neumann boundary conditions. d(phi)/dx=d(phi)/dy=0 at the boundary, for a rectangular domain.

regards, yfyap
  Reply With Quote

Old   January 26, 2001, 02:32
Default Re: help needed in solving Neumann peoblem
  #9
John C. Chien
Guest
 
Posts: n/a
(1). You derived the second order equation based on the assumption the the velocity in 2-D is derived from the gradient of the potential function phi. (2). The consequence of this assumption is that the vorticity, omega, which is defined as the Curl of Velocity, will be zero. In other words, the Curl of (the gradient of phi)=0. (3). So, your assumption automatically lead to non-rotational flow. This category of flows is normally called potential flows. Most aerodynamics or gasdynamics books should cover this subject in great detail. In most cases, the linear potential flow equation is the Laplace equation. So, the boundary conditions also must be specified consistently.
  Reply With Quote

Old   January 26, 2001, 05:49
Default Re: Neumann peoblem compensation scheme proposed
  #10
Dean Schrage
Guest
 
Posts: n/a
Hi,

Ok by float I mean that if the system is iterated without the Neumann cond satisified, the system average pressure (in the case of the PPE) will continue to increase (like float away). However, the gradients of P which are the source terms in the momentum equations are relatively preserved and you can still get a cfd result. I like the second form of your equation. Can recast as a real Poissons equation with source as a function of x, y. I have to go for now. Tracking ISS this morning. see www.d3experiment.homepage.com

Will talk more on this. Hope it's a start.
  Reply With Quote

Old   January 26, 2001, 10:04
Default Re: help needed in solving Neumann peoblem
  #11
Peter Attar
Guest
 
Posts: n/a
It has nothing to do with CFD.It is just a way to solve a linear equation if the system is ill-posed.Do a search online and you will many references to it(actually a better search would be for regularized pseudo-inverse)
  Reply With Quote

Old   January 29, 2001, 10:49
Default Re: Neumann peoblem compensation scheme proposed
  #12
yf yap
Guest
 
Posts: n/a
dear Dean, it is regretful to inform you that the method you suggested is not working. however, i would like you to know that your advice is greatly appreciated. i still have diverging results. for the time being, i re-check the formulation, who knows, maybe the problem is really ill-posed, and i fail to notice that. actually, for a neumann boundary problem akin to the one i'm solving, other than specifying all the normal derivatives at the boundary, what other things that should be specified? regards, yfyap
  Reply With Quote

Old   January 29, 2001, 10:55
Default Re: help needed in solving Neumann peoblem
  #13
yf yap
Guest
 
Posts: n/a
dear Peter, thanks for the advice. regards, yfyap
  Reply With Quote

Old   January 29, 2001, 11:12
Default Re: Neumann peoblem compensation scheme proposed
  #14
Dean Schrage
Guest
 
Posts: n/a
Ok what is the actual Poissons equation you are solving?

you mentioned the following form:

d^2(phi)/dx^2 + D d(phi)/dx + d^2(phi)/dy^2 = Source ----(2) which D=D(x) ,

rewrite this as

GRAD^2(phi) = f(x,y)

where GRAD^2() is the Laplacian and f(x,y) = S(x,y) - D d(phi)/dx which is just a source term. The boundary conditions you describe (grad(phi) = 0) are same as insulation on a heated domain. That is, the flux of Phi induced by gradients is zero at the boundaries. Are you absolutely sure of this? Do you have any boundaries defined by constants (Dirichlet BC). Thus, the source term (f(x,y)) is, by analogy, the same as the source heating on your interior domain. Since the boundaries are insulated, the integration of this function over the inside domain must add to zero. Basically all (+) content of integration must be balanced by (-) content of integration. Heat additions cancel heat extractions. With a Neumman problem, you can have local heating and cooling but they must balance to give zero heat load. They also must be perfect to within machine zero otherwise the temperature in a thermal simulation will float.

What is the function f you describe above. Can we integrate it analytically over the volume domain and check to see if the Neumann condition is upheld?

Next, I would enforce the gradient BC with one sided gradients pointing into the domain. I do this with the PPE. This can cause errors if not implemented correctly. This is tricky to do with ADI codes and care must be taken in its implementation.

Let's work from here

regards

Dean
  Reply With Quote

Old   January 29, 2001, 12:53
Default Re: Neumann peoblem compensation scheme proposed
  #15
yf yap
Guest
 
Posts: n/a
dear Dean, first of all, i think it is most suitable to restate my problem as clear as possible. the continuity equation in cyclindrical coordinate:

d(Ur)/dr + Ur/R + 1/R*[d(Ua)/da] + d(Uz)/dz = 0

note: d()/da, d()/dr, d()/dz are partial derivatives;

R=R(r)

from numerical solution of the a-momentum (a-alpha, angle in radian), d(Ua)/da can be calculated based on the assumption of parabolized Navier-Stokes eqs. besides, by solving the r and z momentum eqs(linearized), Urp and Uzp are obtained, these are predicted values of Ur and Uz, and they do not satisfy the continuity equation. the boundary condition being Urp=Uzp=0 at the wall.

to obtain Ur and Uz that satisfy the continuity eq., i define:

Ur= Urp + Urc

Uz = Uzp + Uzc substitute in to the con. eq. d(Urp+Urc)/dr + (Urp+Urc)/R + d(Uzp+Uzc)/dz = - 1/R[d(Ua)/da] ---[1] please note that the values of Urp and Uzp are already obtained in by solving the linearized r and z momentum eqs. thus, by treating all those known derivatives as source, [1] is rearranged: d(Urc)/dr + (Urc)/R + d(Uzc)/dz = - 1/R[d(Ua)/da] -d(Urp)/dr - (Urp)/R - d(Uzp)/dz d(Urc)/dr + (Urc)/R + d(Uzc)/dz = f(r,z) ----[2] this the f that i meant in my previous discussion, and [2] is a two dimensional problem. in order to solve [2], i define:

Urc=d(phi)/dr , Uzc=d(phi)/dz -----[3]

d()/dr , d()/dz are partial derivatives. substitute [3] into [2]:

d^2(phi)/dr^2 + 1/R*d(phi)/dr + d^2(phi)/dz^2 = f(r,z) --------[4]

this is how i get the eq. for the boundary condition of [4], at the wall, i apply Urc=d(phi)/dr=0 , Uzc=d(phi)/dz=0 (Neumann boundary condition) so that the no-slip condition at the wall be satisfied: Ur=Urp + Urc = 0 + 0 =0 , similar with Uz=0

applying Green's theorem to [4]: ite--> integral of dA --> area dS --> line Dr, Dz, element of r and z ite {d^2(phi)/dr^2 + 1/R*d(phi)/dr + d^2(phi)/dz^2} dA = int {f(r,z)} dA ite { [-d(phi)/dz]Dr + [d(phi)/dr + phi/R]Dz } =int {f(r,z)} dA RHS, line integral over the area bounded by S. since the solution is obtained by numerical method, LHS will not equal RHS, due to truncation errors. define Resi= ite {f(r,z)} dA - ite { [-d(phi)/dz]Dr + [d(phi)/dr + phi/R]Dz } ---[5] however, at the boundary:

d(phi)/dz=0 , d(phi)/dr=0

then [5] reduces to: Resi= ite {f(r,z)} dA - ite { [phi/R]Dz }

define:

delta[f(r,z)]=(Resi/total Area)*(area of each grid point) ----[6]

[6] is subtracted from LHS of [4], thus we solve:

d^2(phi)/dr^2 + 1/R*d(phi)/dr + d^2(phi)/dz^2 = f(r,z) - delta[f(r,z)] -----[7] with d(phi)/dz=0 , d(phi)/dr=0 at the boundary.

This is the formulation of the problem. after several attempts to solved [7] failed(giving divergent solution), i really begin to be doubtful of it. do you agree with this formulation? i am trying to work on your suggestions, it'll takes several days because i am preparing the progress report for the research. thanks a lot for willing to spend your precious time on the matter, your help and advice is greatly appreciated. regards, yfyap
  Reply With Quote

Old   January 30, 2001, 23:59
Default Re: Neumann peoblem compensation scheme proposed
  #16
yf yap
Guest
 
Posts: n/a
dear Dean, first of all, i think it is most suitable to restate my problem as clear as possible. the continuity equation in cyclindrical coordinate:

d(Ur)/dr + Ur/R + 1/R*[d(Ua)/da] + d(Uz)/dz = 0

note: d()/da, d()/dr, d()/dz are partial derivatives;

R=R(r)

from numerical solution of the a-momentum (a-alpha, angle in radian), d(Ua)/da can be calculated based on the assumption of parabolized Navier-Stokes eqs. besides, by solving the r and z momentum eqs(linearized), Urp and Uzp are obtained, these are predicted values of Ur and Uz, and they do not satisfy the continuity equation. the boundary condition being Urp=Uzp=0 at the wall.

to obtain Ur and Uz that satisfy the continuity eq., i define:

Ur= Urp + Urc

Uz = Uzp + Uzc substitute in to the con. eq.

d(Urp+Urc)/dr + (Urp+Urc)/R + d(Uzp+Uzc)/dz = - 1/R[d(Ua)/da] ---[1] please note that the values of Urp and Uzp are already obtained in by solving the linearized r and z momentum eqs. thus, by treating all those known derivatives as source, [1] is rearranged:

d(Urc)/dr + (Urc)/R + d(Uzc)/dz = - 1/R[d(Ua)/da] -d(Urp)/dr - (Urp)/R - d(Uzp)/dz

d(Urc)/dr + (Urc)/R + d(Uzc)/dz = f(r,z) ----[2]

this the f that i meant in my previous discussion, and [2] is a two dimensional problem. in order to solve [2], i define:

Urc=d(phi)/dr , Uzc=d(phi)/dz -----[3]

d()/dr , d()/dz are partial derivatives. substitute [3] into [2]:

d^2(phi)/dr^2 + 1/R*d(phi)/dr + d^2(phi)/dz^2 = f(r,z) --------[4]

this is how i get the eq. for the boundary condition of [4], at the wall, i apply Urc=d(phi)/dr=0 , Uzc=d(phi)/dz=0 (Neumann boundary condition) so that the no-slip condition at the wall be satisfied: Ur=Urp + Urc = 0 + 0 =0 , similar with Uz=0

applying Green's theorem to [4]: ite--> integral of dA --> area dS --> line Dr, Dz, element of r and z

ite {d^2(phi)/dr^2 + 1/R*d(phi)/dr + d^2(phi)/dz^2} dA = int {f(r,z)} dA

ite { [-d(phi)/dz]Dr + [d(phi)/dr + phi/R]Dz } =int {f(r,z)} dA RHS, line integral over the area bounded by S. since the solution is obtained by numerical method, LHS will not equal RHS, due to truncation errors. define

Resi= ite {f(r,z)} dA - ite { [-d(phi)/dz]Dr + [d(phi)/dr + phi/R]Dz } ---[5]

however, at the boundary:

d(phi)/dz=0 , d(phi)/dr=0

then [5] reduces to: Resi= ite {f(r,z)} dA - ite { [phi/R]Dz }

define:

delta[f(r,z)]=(Resi/total Area)*(area of each grid point) ----[6]

[6] is subtracted from LHS of [4], thus we solve:

d^2(phi)/dr^2 + 1/R*d(phi)/dr + d^2(phi)/dz^2 = f(r,z) - delta[f(r,z)] -----[7] with d(phi)/dz=0 , d(phi)/dr=0 at the boundary.

This is the formulation of the problem. after several attempts to solved [7] failed(giving divergent solution), i really begin to be doubtful of it. do you agree with this formulation? i am trying to work on your suggestions, it'll takes several days because i am preparing the progress report for the research. thanks a lot for willing to spend your precious time on the matter, your help and advice is greatly appreciated. regards, yfyap
  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
Moving mesh Niklas Wikstrom (Wikstrom) OpenFOAM Running, Solving & CFD 122 June 15, 2014 07:20
Upgraded from Karmic Koala 9.10 to Lucid Lynx10.04.3 bookie56 OpenFOAM Installation 8 August 13, 2011 05:03
Orifice Plate with a fully developed flow - Problems with convergence jonmec OpenFOAM Running, Solving & CFD 3 July 28, 2011 06:24
Differences between serial and parallel runs carsten OpenFOAM Bugs 11 September 12, 2008 12:16
Could anybody help me see this error and give help liugx212 OpenFOAM Running, Solving & CFD 3 January 4, 2006 19:07


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