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

(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.

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.

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

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

dear Peter Attar, thanks for the advice. what is SVD-regularization method ? i haven't came across this method before. thanks regards, yfyap

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

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

(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.

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.

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)

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

dear Peter, thanks for the advice. regards, yfyap

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

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

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