
[Sponsors] 
Discrepancy in the Strouhal number of a flow past circular cylinder 

LinkBack  Thread Tools  Search this Thread  Display Modes 
March 21, 2017, 17:04 
Discrepancy in the Strouhal number of a flow past circular cylinder

#1 
Senior Member
Hector Redal
Join Date: Aug 2010
Location: Madrid, Spain
Posts: 243
Rep Power: 16 
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. 

March 21, 2017, 17:19 

#2 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,781
Rep Power: 71 
You apriori 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 nonsymmetric 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. 

March 22, 2017, 00:54 

#3 
Senior Member
Uwe Pilz
Join Date: Feb 2017
Location: Leipzig, Germany
Posts: 744
Rep Power: 15 
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.
__________________
Uwe Pilz  Die der Hauptbewegung überlagerte Schwankungsbewegung ist in ihren Einzelheiten so hoffnungslos kompliziert, daß ihre theoretische Berechnung aussichtslos erscheint. (Hermann Schlichting, 1950) 

March 22, 2017, 05:57 

#4  
Senior Member
Hector Redal
Join Date: Aug 2010
Location: Madrid, Spain
Posts: 243
Rep Power: 16 
Quote:
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 nondependent stationary case and the dependent stationary case. 

March 22, 2017, 06:04 

#5  
Senior Member
Hector Redal
Join Date: Aug 2010
Location: Madrid, Spain
Posts: 243
Rep Power: 16 
Quote:
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 

March 22, 2017, 06:37 

#6 
Senior Member
Uwe Pilz
Join Date: Feb 2017
Location: Leipzig, Germany
Posts: 744
Rep Power: 15 
> 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: Code:
snGradSchemes { default corrected 1; } Code:
snGradSchemes { default corrected 0.333; }
__________________
Uwe Pilz  Die der Hauptbewegung überlagerte Schwankungsbewegung ist in ihren Einzelheiten so hoffnungslos kompliziert, daß ihre theoretische Berechnung aussichtslos erscheint. (Hermann Schlichting, 1950) 

March 22, 2017, 06:38 

#7  
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,781
Rep Power: 71 
Quote:
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. 

March 22, 2017, 06:59 

#8 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,781
Rep Power: 71 
just to add that the boundaries of the domain should be located quite far such to not interfer with the flow


March 22, 2017, 15:14 

#9 
Senior Member
Hector Redal
Join Date: Aug 2010
Location: Madrid, Spain
Posts: 243
Rep Power: 16 
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.


March 22, 2017, 15:24 

#10 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,781
Rep Power: 71 
That's not so simple! You will run much more timesteps that cumulate a certain final error. First order time integration is not suitable for transient problems.


March 22, 2017, 16:11 

#11  
Senior Member
Hector Redal
Join Date: Aug 2010
Location: Madrid, Spain
Posts: 243
Rep Power: 16 
Quote:
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 timesteps, 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). 

March 22, 2017, 16:47 

#12 
Senior Member
Hector Redal
Join Date: Aug 2010
Location: Madrid, Spain
Posts: 243
Rep Power: 16 
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) 

March 22, 2017, 16:51 

#13  
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,781
Rep Power: 71 
Quote:
No, you have to consider the approximation of the first time derivative, not the increment, and the derivative is O(dt) 

April 6, 2017, 18:20 

#14  
Senior Member
Hector Redal
Join Date: Aug 2010
Location: Madrid, Spain
Posts: 243
Rep Power: 16 
Quote:
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. 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Flow past a 2D cylinder  High Re (1E+05)  Cd too high  Pervispasco  OpenFOAM Running, Solving & CFD  4  March 14, 2022 02:19 
[snappyHexMesh] Error snappyhexmesh  Multiple outside loops  avinashjagdale  OpenFOAM Meshing & Mesh Conversion  53  March 8, 2019 09:42 
Drag force coefficient too high for a flow past a cylinder using komega sst  Scabbard  OpenFOAM Running, Solving & CFD  37  March 21, 2016 16:16 
Compressor Simulation using rhoPimpleDyMFoam  Jetfire  OpenFOAM Running, Solving & CFD  107  December 9, 2014 13:38 
[blockMesh] > foam fatal error:  lillo763  OpenFOAM Meshing & Mesh Conversion  0  March 5, 2014 10:27 