[samycfd]

NO ! But I still get a result. But it has not converged. I am using the conjugate gradient solver (cgs) and I get :

cgs stopped at iteration 1500 without converging to the desired tolerance 1e-006 because the maximum number of iterations was reached.
The iterate returned (number 119) has relative residual 0.067

Normally this solver is fast and if after 1500 iterations it has not converged it will never converge I guess. ( with 5000 iterations same residual )

[FMDenaro]

Originally Posted by samycfd

NO ! But I still get a result. But it has not converged. I am using the conjugate gradient solver (cgs) and I get :

cgs stopped at iteration 1500 without converging to the desired tolerance 1e-006 because the maximum number of iterations was reached.
The iterate returned (number 119) has relative residual 0.067

Normally this solver is fast and if after 1500 iterations it has not converged it will never converge I guess. ( with 5000 iterations same residual )

this is an indication that your BC.s do not satisfy the compatibility relation...
Are you using homogeneous Neumann bc.s?

[samycfd]

Ah okay
I am using homogeneous Neumann BCs, that's why I was thinking that a value in one point needed to be defined.

[FMDenaro]

Originally Posted by samycfd

Ah okay
I am using homogeneous Neumann BCs, that's why I was thinking that a value in one point needed to be defined.

if the compatibility relation is satisfied the solver converges without that...
you can check that by integrating over the volume both LHS (which will be zero) and RHS and checking if the integral of the source term (RHS) is zero

[Parth04]

need help about implementing orlanski's (non reflecting) boundary codition at outflow boundary. . .
also , what should be b.c. for pressure, when velocity b.c. is orlanski's b.c. . . .
solving stratified flow problem using FVM

thnx in advance

[FMDenaro]

Originally Posted by Parth04

need help about implementing orlanski's (non reflecting) boundary codition at outflow boundary. . .
also , what should be b.c. for pressure, when velocity b.c. is orlanski's b.c. . . .
solving stratified flow problem using FVM

thnx in advance

this topic concerns BC.s for elliptic equation, you are actaully talking about convective-type of BC.s

[samycfd]

In my case don't I integrate over the surface the RHS and the LHS ?
I get the idea even I am not sure of how doing that. Numerical integration with a Gauss quadrature approximation for example ?

[FMDenaro]

Originally Posted by samycfd

In my case don't I integrate over the surface the RHS and the LHS ?
I get the idea even I am not sure of how doing that. Numerical integration with a Gauss quadrature approximation for example ?

you have to do a volume integral, the LHS becomes a surface integral as

Int [V] Div Grad P dV = Int [SV] dp/dn dS = 0

The volume integral of the source term must be zero, even if you can write as surface integral according to the LHS.

[samycfd]

Okay.
I used the trapezoidal approximation function of Matlab to integrate the RHS and it's around 9.0.
edit: the lid velocity has been set to 1 m/s .
This is my RHS :

[FMDenaro]

when computing the RHS of the pressure equation you need the values of the normal velocity along the walls, did you consider it is zero? The previous plot of the velocity seems not reporting the wall but only interior points

[samycfd]

when I export the velocity from FLUENT the values are located on the cell-center, so near the edges I won't have 0 for the normal velocity or 1 for the lid velocity because it's not really the edges but a little bit to the interior.

I tried to add boundaries to the velocity field( 0 everywhere for v and the same for u excepted on the lid where the velocity is set to 1 ). When I am doing that I don't need to set a Dirichlet point for pressure. But the RHS integral is different from 0. (10^4)

I will post some results so you can see.

[FMDenaro]

how do you construct the discrete source term in the pressure equation?
If you set the exact normal velocity component and have the correct discretization of the source term it should actually work and converge without fixing a value.

[samycfd]

I use finite difference to compute the first derivatives of u and v along x and y. and then I construct : RHS = mu/rho * ( Ux^2 +Vy^2 + 2*Uy*Vx ).

[FMDenaro]

Originally Posted by samycfd

I use finite difference to compute the first derivatives of u and v along x and y. and then I construct : RHS = mu/rho * ( Ux^2 +Vy^2 + 2*Uy*Vx ).

but how do you work on the boundary when the dp/dn=0 BC is applied?

[samycfd]

it could sound weird but I am doing nothing because I am using FEM and the Neumann boundaries ( dp/dn=0) are naturally applied.

Take a look at the variationnal problem ( weak formulation ) of div( grad P ) = RHS

int (grad P grad phi) - int ( dp/dn )|sur le bord = int ( RHS phi )

dp/dn = 0

finally we just have :
int (grad P grad phi) = int ( RHS phi )

And we can write this as a linear system KP = RHS where K is a stiffness matrix.

But maybe I am wrong...


And I want to add that I don't know exactly how dealing with the BC during the computation of derivatives with finite differences. At the moment I am using Neumann everywhere. So I am just copying the value to the edge so the gradient on the would be zero.

[FMDenaro]

ok, I'll try to give idea of how I would proceed in a Finite Volume approach....

Assume dv/dt=a=0, so that Div a=0. Assume Div (vv)=ac and Div (2ni Gradv)=ad.
Hence the momentum equation simply writes

ac + Grad P = ad


Therefore one gets
Div a= -Div ac-Div Grad P + Div ad=0.

This leads to the pressure equationDiv Grad P = Div (ad-ac)

with the BC.s

dP/dn = n.(ad-ac)

As a consequence, if you set dP/dn =0 then you must set also n.(ad-ac)=0.
For example, assume a 1D mesh, in the node i you have to discretize the pressure equation, with second order discretization one has

[(dP/dx)i+1/2 - (dP/dx)i-1/2]/dx = [(axd-axc)i+1/2 - (axd-axc)i-1/2]/dx


Assume for example that now the section i-1/2 lies on a wall where you have to set the BC.s. Consequently the previous equation writes


[P(i+1) - P(i)]/dx^2 = [(axd-axc)(i+1) +(axd-axc)(i)]/2dx

This is the correct implementation of the BC that, as you see, modifies also the source term in a way that the compatibility relation is satisfied.

[samycfd]

In the finite volume the source term is also modified, I agree but not in the Finite Element Method referring to a lot of papers about Neumann BC. I have never used the finite volume method. Do you advise me to use it instead of Finite Element?

When I was dealing only with rectangular grid I was using the Finite differences and I was modifying the souce term and the matrix of laplacian to implement Dirichlet or Neumann BC. ( And now I know that I have to use any kind of mesh I have to adapt the method and I choosed Finite element ... )

[FMDenaro]

Originally Posted by samycfd

In the finite volume the source term is also modified, I agree but not in the Finite Element Method referring to a lot of papers about Neumann BC. I have never used the finite volume method. Do you advise me to use it instead of Finite Element?

When I was dealing only with rectangular grid I was using the Finite differences and I was modifying the souce term and the matrix of laplacian to implement Dirichlet or Neumann BC. ( And now I know that I have to use any kind of mesh I have to adapt the method and I choosed Finite element ... )

I don't see how is possible theoretically that you set dP/dn=0 without changing (congruently) also the source term near the wall. Maybe you are referring as to some method where the correction is implicitly implemented by the shape function... but remember the FV is only a special case of FEM where you use a special step-wise shape function...therefore I suppose something is implicitly modified.

However, try using the procedure I suggested instead of discretizing the pressure equation directly in FD term

[samycfd]

I understand your point and I am starting to believe that I am wrong and of course it's wrong given the results.. .

I would like to, but I don't understand the whole point

1)

[P(i+1) - P(i)]/dx^2 = [(axd-axc)(i+1) +(axd-axc)(i)]/2dx

Isn'it ? [P(i+1) - P(i)]/dx^2 = [(axd-axc)(i+1) +(axd-axc)(i-1)]/2dx

2) Finally you discretize the equation with FD in all interior nodes and apply a specific treatment for this discretization on the boundaries by set all the dp/dn value on edges to 0 ?

[FMDenaro]

1) No... the RHS is written by setting dP/dx|i-1/2= (axd-axc)i-1/2=0 at wall.
Thus the RHS becomes [(axd-axc)i+1/2]/dx. Then, assuming you know the values only at nodes in centroid (...i-1,i,i+1...) of the cells you can write a linear interpolation between i+1 and i

[(axd-axc)i+1/2]/dx=0.5*[(axd-axc)i+1 + (axd-axc)i ]/dx


2) I discretize in the interior nodes in a similar way... for example, in the 1D lines I write

[(dP/dx)i+1/2 - (dP/dx)i-1/2]/dx = [(axd-axc)i+1/2 - (axd-axc)i-1/2]/dx

as

[P(i+1)-2P(i)+P(i+1)]/dx^2 = [(axd-axc)i+1 - (axd-axc)i-1]/2dx