[kabz]

Hi,

I've recently coded up in C++ the SIMPLE algorithm for incompressible flow btn parallel plates (u-momentum eq'n and mass eq'n). I'm using a staggered grid and following the book by Versteeg and Malalasekera (2007).

I've run a few test cases with different BC's:
- inlet velocity (uniform) with pressure outlet (0) -> uniform velocity distribution throughout flow is attained.
- inlet velocity (parabolic) with pressure outlet (0) -> parabolic velocity distribution throughout flow is attained.


Other BCs are u=0 at top and bottom wall (u-momentum eq'n) and p'=0 for outlet. The pressure correction equation for top and bottom wall give p'=0 since the source term is always zero (u=0). And the pressure at the inlet is obtained by linear interpolation from inside the domain (Ferziger & Peric, 1997)

My question is on the first case - is it correct? Isn't the flow supposed to develop into a parabolic profile towards the outlet?

[FMDenaro]

Originally Posted by kabz

Hi,

I've recently coded up in C++ the SIMPLE algorithm for incompressible flow btn parallel plates (u-momentum eq'n and mass eq'n). I'm using a staggered grid and following the book by Versteeg and Malalasekera (2007).

I've run a few test cases with different BC's:
- inlet velocity (uniform) with pressure outlet (0) -> uniform velocity distribution throughout flow is attained.
- inlet velocity (parabolic) with pressure outlet (0) -> parabolic velocity distribution throughout flow is attained.


Other BCs are u=0 at top and bottom wall (u-momentum eq'n) and p'=0 for outlet. The pressure correction equation for top and bottom wall give p'=0 since the source term is always zero (u=0). And the pressure at the inlet is obtained by linear interpolation from inside the domain (Ferziger & Peric, 1997)

My question is on the first case - is it correct? Isn't the flow supposed to develop into a parabolic profile towards the outlet?

yes, provided that the lengh of the plate is sufficient, you should get a parobolic velocity profile. What Re number do you set?

[kabz]

Originally Posted by FMDenaro

yes, provided that the lengh of the plate is sufficient, you should get a parobolic velocity profile. What Re number do you set?

Re=10000
density=1000Kg/m^3
viscosity=0.001Kg/m s
height=0.1m
length=width=1m
velocity=0.05m/s

[FMDenaro]

at such a Re number you cannot see a developed viscous profile........you have to reduce the Re number of some order of magnutude

[kabz]

Originally Posted by FMDenaro

at such a Re number you cannot see a developed viscous profile........you have to reduce the Re number of some order of magnutude

Tried it with lower velocities (u=1e-3, u=1e-4) but I still obtain a uniform flow profile.

[FMDenaro]

Originally Posted by kabz

Tried it with lower velocities (u=1e-3, u=1e-4) but I still obtain a uniform flow profile.

Try for Re=1

[kabz]

I still obtain a uniform velocity profile.

the predicted velocity profile is okay but once i solve the pressure correction equation and correct the velocity, i obtain a uniform profile. See the attached matlab .fig files

[FMDenaro]

Originally Posted by kabz

I still obtain a uniform velocity profile.

the predicted velocity profile is okay but once i solve the pressure correction equation and correct the velocity, i obtain a uniform profile. See the attached matlab .fig files

At Re=1 you should see the viscous parablic profile...something is wrong...are you sure to have reached the steady state?

[kabz]

Originally Posted by FMDenaro

At Re=1 you should see the viscous parablic profile...something is wrong...are you sure to have reached the steady state?

Yes. At 200 iterations, only the predicted velocity changes but the corrected velocity is maintained. See attached.

[FMDenaro]

please, post directly the image... however, there is for sure something wrong in the code...check the BC.s after the correction, are they no-slip conditions?

[kabz]

Originally Posted by FMDenaro

please, post directly the image... however, there is for sure something wrong in the code...check the BC.s after the correction, are they no-slip conditions?

Yes no slip conditions are maintained

[FMDenaro]

you have a bug in the code ... at Re=1 you should see the two spatially evolving boundary layers

[kabz]

After much debugging, i have some limited success - see images attached.
I get a more or less triangular profile near the outlet

Ignore the pressure at inlet as I've not implemented the linear interpolation from inside the domain.

[FMDenaro]

you found the bug, the results seems physically reasonable now...you can compare with the analytical Blasius solution for a flat plane...then check that after some height a fully parabolic velocity profile is obtained
Of course, x=0 is a singular point.

[kabz]

Originally Posted by FMDenaro

you found the bug, the results seems physically reasonable now...you can compare with the analytical Blasius solution for a flat plane...then check that after some height a fully parabolic velocity profile is obtained
Of course, x=0 is a singular point.

My current issue is that it seems to develop a parabolic profile after entry but towards the outlet a "triangular" profile is acquired.

[FMDenaro]

Originally Posted by kabz

My current issue is that it seems to develop a parabolic profile after entry but towards the outlet a "triangular" profile is acquired.

plot the velcoity at some station and superimpose the parabolic law...anyway, do you have ensured to have a suffient lenght do let the flow develop?

[kabz]

Hopefully final question:

The momentum equations are linearized by using velocity and pressure fields from the previous outer iteration to solve for new fields. Now, at the end of an outer iteration, one has the new corrected velocity and pressure fields - u^{n},v^{n},p^{n}. Are residuals calculated based on the linearized momentum equations (ie using the old velocity fields to compute coefficients and old pressure field to compute flux and then apply the new velocity field to complete the equation) or are the residuals based on the non-linear momentum equations (ie use the new velocity fields to compute coefficients and new pressure field to compute flux and then apply new velocity field to complete the equation)?

[FMDenaro]

the residual is computed for the discrete equation you really want to satisfy...if you linearized, this discrete equation must provide a vanishing residual. You cannot use the solution vector to compute the residual in a different discrete equation

[kabz]

Originally Posted by FMDenaro

the residual is computed for the discrete equation you really want to satisfy...if you linearized, this discrete equation must provide a vanishing residual. You cannot use the solution vector to compute the residual in a different discrete equation

And if one used deferred correction terms (eg bounded higher order schemes), are these computed for residuals or is it safe to ignore them?

[FMDenaro]

Originally Posted by kabz

And if one used deferred correction terms (eg bounded higher order schemes), are these computed for residuals or is it safe to ignore them?

in the residual you must consider any discrete term you used in the computation... in symbolic term, said vn the solution at step n and A the discrete operators containing all fluxes (bounded or not), the residual at step n can be written (similar to interative methods for linear algebric systems)

A(vn) - s = rn


somehow, for steady solutions, rn can be interpreted as the time derivative you want to drive to zero