[HectorRedal]

Hi,

I am simulating a flow past a circular cilinder of an incompressible flow.
The characteristics of the flow I am using are the following:
- Re = 100
- Incompressible flow
- Isothermal flow (without energy exchange).
- Horizontal flow (perpendicular to the cilinder.
- Gravity force = 0 (No gravity is considered in the simulation).

According to the references I have checked, the vallue of the Strouhal Number for this case should be 0.164 - 0.165.

The value of the Strouhal number I am obtaining is 0.1822, which means a 10% of error.

It appears as the viscosity of the simulation is higher than what specified. Then I thought that some artificial viscosity was added by the implementation of the algorithm.
I try to use a lower time step, and even a finer grid, but the Strouhal number obtained has not differ at all.

When Re < 40, the flow is stationary (no oscilation is observed), and my results agree with any references I have checked (Eddy length at the wake of the cilinder, and angle of separation). The difference I am obtained for the Eddy length and the angle of separation are less than 0.1%.

These second results made me thought that there may not be such artificial viscosity in the implemenation of the algorithm, since if this was the case, then this artificial viscosity should be observed in this second test case.

Am I wrong?

Best regards,
Hector.

[FMDenaro]

You a-priori should know if your discretization has artificial viscosity, no matter of how many solutions you try to discover...

Artificial viscosity is a consequence of a non-symmetric discretization of the convective terms, coupled with the local truncation error of the time integration.

Thus, if you know how your code discretize the equations, you know if the solution has some amount of artificial viscosity.

[piu58]

The numerical solution depend strongly on the scheme you use for differentiation. Is is worth to try out some things and look to what amount the final solution changes. There with you get some feeling which accuracity you may expect.

In the schemes you have to blance numerical stability / boundedness vs. mathematical accuracy.

[HectorRedal]

Originally Posted by FMDenaro

You a-priori should know if your discretization has artificial viscosity, no matter of how many solutions you try to discover...

Artificial viscosity is a consequence of a non-symmetric discretization of the convective terms, coupled with the local truncation error of the time integration.

Thus, if you know how your code discretize the equations, you know if the solution has some amount of artificial viscosity.

The discretization of the convective term does not have any artificial viscosity, a-priori.
Up to my knowledge, the discretization of the convective term is second order accurate, O(h^2):
C * U + 1/2 * K U^2

Being,
U the velocity
C the convective discretization matrix
K the second order convective discretization matrix.

Regarding the time discretization is first order accurate, O(h).

My question is more related to if from the behaviour observed, one can deduce there is an error in the implementation of the algorithm, taking as input the difference between the non-dependent stationary case and the dependent stationary case.

[HectorRedal]

Originally Posted by piu58

In the schemes you have to blance numerical stability / boundedness vs. mathematical accuracy.

Your comment is quite interesting from my point of view.
You are comparing balance numerical stability vs boundedness. Can you elaborate a bit more your statement?

On the other hand, the algorithm I am using have two parameters related to the velocity and pressure value: Relaxations factors.

I am a bit puzzle about the influence of this factors in the solution obtained by the algorithm.

Basically this relaxation factors work in this way:
Velocity = Velocity (at t= tn) * (1 - theta1) + Velocity (at t=tn+1) * theta1
Pressure = Pressure (at t = tn) * (1 - theta2) + Pressure (at t=tn+1) * theta2

Being,
theta1 = the relaxation factor for the velocity
theta2 = the relaxation factor for the pressure.

My question is which value of this parameters to use in the simulations?
Based on experience? If you don't have experience, which values to use?

Right now, I am using theta1= 1.0 and theta2 = 1.0

[piu58]

> You are comparing balance numerical stability vs boundedness.

For calculation of the field you need the temporal and spatial derivation of the flux values. These are calculated by differences and schemes for this purpose. This is straightforward for a rectangular mesh. It gets more complicated for skewed meshes: The schemes contain portions for correction the geometric mesh properties.

If these correction terms are 'correct' in a mathematical sense they may lead to unbounded or in extreme infinite results. An example: It is possible to have two adjacent cell with the same cell center (if we count the center of mass as cell center). The distance between the cells comes to zero an a spatial derivation gets infinity. An example for hex cells:



Most cases are no nit sch worse, but the problem remains: Correction for geometry may lead to instabile or oscillating results. It is wise to dampen the deviations at least for the most worse parts of the mesh. The dampened result is not correct anymore, so you get numerical effects, if you increase stability.

One example for calculation of normal vectors (normal to the cell patches). The vector needs to be corrected for skewed cells:

snGradSchemes {
    default         corrected 1;  
}

This corrects mathematical correct, but may lead to very large values for the vector (unbounded).

snGradSchemes {
    default         corrected 0.333;  
}

In this case the correction is limited to the half the value of the one for orthogonal cells.

[FMDenaro]

Originally Posted by HectorRedal

The discretization of the convective term does not have any artificial viscosity, a-priori.
Up to my knowledge, the discretization of the convective term is second order accurate, O(h^2):
C * U + 1/2 * K U^2

Being,
U the velocity
C the convective discretization matrix
K the second order convective discretization matrix.

Regarding the time discretization is first order accurate, O(h).

My question is more related to if from the behaviour observed, one can deduce there is an error in the implementation of the algorithm, taking as input the difference between the non-dependent stationary case and the dependent stationary case.

I do not understand your type of discretization, it is not clear the stencil you are using. Then, first order accuracy in time ( O(dt) ) is totally inadeguate for a transient simulation and introduce large truncation error.

[FMDenaro]

just to add that the boundaries of the domain should be located quite far such to not interfer with the flow

[HectorRedal]

Originally Posted by FMDenaro

I do not understand your type of discretization, it is not clear the stencil you are using. Then, first order accuracy in time ( O(dt) ) is totally inadeguate for a transient simulation and introduce large truncation error.

This limitation can be overcome reducing the time step. This way the O(dt) becomes small. If this is the problem, the accuracy of the simulation can be imprved easily (reducing the delta t). The penalty is longer simulation times.

[FMDenaro]

Originally Posted by HectorRedal

This limitation can be overcome reducing the time step. This way the O(dt) becomes small. If this is the problem, the accuracy of the simulation can be imprved easily (reducing the delta t). The penalty is longer simulation times.

That's not so simple! You will run much more time-steps that cumulate a certain final error. First order time integration is not suitable for transient problems.

[HectorRedal]

Originally Posted by FMDenaro

That's not so simple! You will run much more time-steps that cumulate a certain final error. First order time integration is not suitable for transient problems.

The algorithm I have implemented is the CBS Scheme:
https://www.researchgate.net/publica..._CBS_algorithm

This algorithm has been used for simulating both static and transient problems (in the reference above you can find several transient problems that have been run and they provide very good agreement with validation test cases).

I can think of adapting the algorithm to implement a second order time integration scheme, but this will require time.

On the other hand, I have not as much experience as you on this subject, and I agree with you that you will run much more time-steps, but if you choose a proper time step so that the sum of O(delta t1) = O (delta t2 ^2) (where t2=sum of t1), I undersand that the error at the end of the simulation should be the same (or approximately).

[HectorRedal]

Sorry, but I have realized that the CBS Algorithm is a second order time approximation algorithm:
U(n+1) = U(n) - U(n) * delta T * (d U(n) / dx) + O (delta T^2)

[FMDenaro]

Originally Posted by HectorRedal

Sorry, but I have realized that the CBS Algorithm is a second order time approximation algorithm:
U(n+1) = U(n) - U(n) * delta T * (d U(n) / dx) + O (delta T^2)

No, you have to consider the approximation of the first time derivative, not the increment, and the derivative is O(dt)

[HectorRedal]

Originally Posted by FMDenaro

No, you have to consider the approximation of the first time derivative, not the increment, and the derivative is O(dt)

I have double-checked again the reference and I have realized that the velocity is approximated using the following formula:
U_n+1 =U_n + delta t *(-d (u_j * rho * U_i) / dx_j + d tau_ij / dx_j - rho * g_i + delta t / 2 * u_k * d(d (u_j * rho * Ui) / dx_i + rho g_i)).
I am including a copy of the formula (attached).

So, I think that the algorithm provides an aproximation of second order in time.