[LAH]

Dear forum users,

I am a student and relatively new to CFD. My goal is to simulate in 2D, the flow over a square cylinder at incidence angles (alpha) from 0 to 90 degrees using increments of 5 degrees. To save some time, my plan was to simulate from alpha = 0 to alpha = 45 and 'mirror' the results to obtain results for alpha = 45 to alpha = 90.

Since this requires 10 simulations to be performed I did not want to have to create a separate mesh for each case. Instead, I change the direction of the inflow vector while using the same mesh with the help of defined variables in CFX pre (shown in one of the images). The sides of the square cylinder are 200mm and the velocity magnitude is 20m/s resulting in a Re of approximately 250,000.

I run the simulation in transient mode from 0 to 6s using adaptive timestepping. The initial timestep is 1e-5s after which I let CFX alter the timestep to reach convergence in 3 to 5 iterations. I save a cgns file containing surface pressure every 3ms for later post processing. Regardless of the angle of incidence, 'steady' oscillations are reached after roughly 1.5s so I use the cgns files from 1.5s to 6s to obtained time averaged surface pressures. Also important: I use the k-omega SST model.

The mesh has approximately 120,000 elements, a Y+ < 1 and growth ratio of 1.1. I refine the mesh towards the corners of the cylinder to x+ = 1 as I expect large gradients there.

To me, this setup seemed suitable for what I want to obtain (time averaged surface pressure from alpa = 0 to 45) and thus I ran all the required simulations. To quickly check on the results I both monitor during the simulation and calculate after the simulation the drag coefficient by integrating the pressure over the cylinder surface. For alpha = 0, 5, 10, both the Cd and the pressure distribution are close to experimental data (Cd at alpha = 0 approximately 2 for example) above these angles my found Cd values are much larger than experimental ones. At alpha = 45 degrees, I find a Cd of approximately 3.1 where it should be about 2. I use L = 0.2m as the reference length for Cd calculations and thus tought at first that the difference was due to some papers using the projected length instead, however many papers describe a Cd of approx 2 using the cylinder side length as a reference. A similar discrepancy is present for all angles between 10 and 45 degrees.

Am I clearly missing something/doing something wrong in my setup/meshing or are there any recommendations for simulating flow over this geometry. I noticed that even though I use adaptive timestepping with a minimum timestep of 1e-5s (CFL = 20) and a maximum of 1e-3s, after a while all simulations use a timestep of 1e-3s and still converge within 3 to 5 iterations. This is great as each simulation then takes 4 hours for me on the cluster however means that the CFL >= 999 could this lead to problems even though the timestepping is implicit?

[ghorrocks]

The general FAQ on accuracy is here: https://www.cfd-online.com/Wiki/Ansy..._inaccurate.3F

Specific comments for your case:
* You need to do sensitivity checks on mesh size, time step and convergence criteria. I suspect your mesh is OK (but you should check anyway), using adaptive timestepping almost always results in a good time step (but again, check it) - but you have not done a convergence criteria check. Rerun it with a convergence criteria a factor of 10 tighter and see if it makes a difference.

* Re = 2.5E5 is a tricky region to simulate as the boundary layer on the body is likely to be laminar and the wake is turbulent. You might get a better result if you run this laminar rather than turbulent. If you run it turbulent you will definitely need a turbulence model which handles low Re turbulence well such as SST, k-w. The k-e model will really struggle here. You may also need to consider a LES type approach, possible SAS or DES.

* There may also be 3D effects which are significant. I would try a 3D model to see if this is the case.

[LAH]

Thank you Glenn Horrocks for you reply,

I have run the case on many different meshes, near wall refinement from Y+ 100, 50, 25, 5 and < 1. I have also tried different boundary conditions (opening instead of outlet) and extended the domain size. All without any significant effect on the results and still the large drag coefficients. I even tried rounding the corner edges for possible errors as a result of sharp corners.

As you suggested I have now also run the case with a convergence criteria of 10^-5 instead of 10^-4, this still resulted in a drag coefficient of 3.1 for alpha 45 which I use as a test case where it should be about 2.

I have also run a simulation using the SAS and DES turbulence models which again resulted in similar results with drag coefficients close to 3 for alpha 45. (Note, these simulation were done with the same adaptive timestepping settings as shown before)

I am currently running a simulation with a timestep such that the CFL number does not exceed 20 to rule out if a too large time step is causing the unexpected results. This will take a while but if it doesn't work I could try a 3D simulation as a last effort.

Despite the large drag coefficients, the results do not look unphysical. The obtained time averaged pressure distributions are all smooth and show expected behaviour for varying the angle of attack, there just seems to be a change in the pressure distribution behaviour for varying alpha from alpha = 15 onwards. Also, the large drag seems to mainly be caused by much too low pressures on the 'back' of the cylinder as seen from the flow direction. The front side pressure distribution is similar to experimental results regardless of angle of attack as well as the pressure distribution on the sides relative to the airflow.

[Gert-Jan]

Originally Posted by LAH

I run the simulation in transient mode from 0 to 6s using adaptive timestepping. The initial timestep is 1e-5s after which I let CFX alter the timestep to reach convergence in 3 to 5 iterations.

- What is your criterium for convergence?
- Are you monitoring the convergence within a timestep (CFX-Pre > Output Control > Monitor Convergence Loops Coefficents)
- Do you monitor variables at monitoring points?
- If so, are these flat liners within each timestep?
- Same for the Cd

You should make sure you reach convergence within a time step. From your story, I am not convinced you have that.

[LAH]

My residual target is the default 1e-4 as shown in one of the attached images. As mentioned in the previous post, I tried running the simulation with this target a factor 10 smaller without significant effects.

I am not sure I completely understand what you mean about convergence within a timestep but I have not considered this yet.

I enabled the monitor coefficient loop convergence setting and restarted a run that previously finished after 6 seconds of simulation time. Image solver_1 shows the monitor 'variable' drag which is calculated as: (areaInt_x(Pressure) @ Cylinder)*cos(alpha [deg]) + (areaInt_y(Pressure) @ Cylinder)*cos(90 [deg] - alpha [deg]) where alpha is 45 [deg]. solver_1 shows the last 100 iterations. Image solver_2 also shows another variable dp which is the difference between the average pressure at both inlets w.r.t. both outlets. Image solver_3 shows the residuals of P-mass, U-mom and V-mom of the last 100 iterations.

If I understand correctly I should see a clear 'stairstepping effect' when enabling monitor coefficient loop convergence indicating that withing each timestep the variable is constant over the inner iterations, correct? I also attached image solver_4 which shows only the behaviour of the monitor points drag and dp over the last timestep. Within 1 timestep both drag and dp are almost flat over the inner iterations except for the slight change that naturally occurs as time progresses slightly.

Is there something off about these convergence graphs?

[ghorrocks]

A residual target is one thing, but unless the residuals actually get to the tighter target you have not converged tighter. Did the simulation achieve the tighter residuals target?

Gert-Jan's comment about convergence within a timestep is saying that you want values of importance to you to converge inside each time step. If they are still changing and you declare the time step converged then the values of importance to you have not converged and your accuracy will be poor.

It is better to use the force_x,y,z functions rather than areaInt of pressure. The force function includes wall friction which areaInt of pressure does not, and the integration is done more accurately.

Can you post an image of what the flow looks like? Also the output file would be good (to make it a manageable size just do a run with a few time steps)

Also - try running this laminar as I suggested. Also also - I note you are using first order turbulence numerics, you might try second order.

[LAH]

Thank you for all the suggestions.

I forgot to mention it in the last post but I tried a laminar case. An image of an instantaneous velocity field is attached. This velocity field looks a bit weird and there were some warnings about backflow occurring through outlets only for this case. Still the time averaged Cd was around 3 again.

I completely agree that using a sum of forces instead of pressure integral would be more accurate for the Cd calculation. I am going to repeat these simulations in a windtunnel as well where I will be estimating pressure drag as a pressure integral around the cylinder which is why I choose to do that for the simulations as well. For bluff bodies like this, the friction drag is relatively small anyway which I also confirmed when looking at pressure- vs friction drag in CFD post after the simulation.

I have also tried running the case using high order turbulence numerics with no significant effect on the flowfield and drag coefficients which again oscillate around 3.

The velocity and pressure images attached are instantaneous velocity and pressure fields obtained using the SST k-omega turbulence model at a time after oscillations have started. High velocity 'eruptions' seem to occur near the trailing edge of the cylinder under this angle which clearly shows up as a very low pressure region. This low pressure region oscillates from side to side and clearly shows up in the time averaged pressure distribution (could not add an image of the pressure distribution as this message contains 5 files already), likely causing excessive drag.

I have also attached 2 output files in .txt format which hopefully can give you a better understanding about how the simulation is running. The Yp_1_alpha_45_SST.txt file contains about 30 timesteps of a run using the SST k-omega model with a convergence criteria of 10e-4 using adaptive timestepping. The Yp_alpha_45_SST_RMS.txt file also contains about 30 timesteps executed using the same settings except for the convergence criteria target of 10e-5 instead. Both these output files where obtained after initiating the simulation from a previously obtained .trn file to make sure that the vortex shedding state is reached (t > at least 1.5s).

Please let me know if there is something suspicious in the output files that I have missed. Also, if there is any other settings/figures that I should share to make things clearer.

[ghorrocks]

The plot thickens

So it appears that the source of the Cd error is the pressure on the down stream side being too low. The convergence, time step, turbulence model do not seem to affect it. I think it is time to take a step back and check some basics.

If you model flow over a circular cylinder at similar Re, do you get an accurate result? Also flow over a flat plate (with the plate perpendicular to the flow, so it is a big separation? Both of these cases require accurate modelling of the pressure in the wake.

Given that everything else seems to not affect it the thing I am suspicious of now is your mesh quality. But let's see if you can model flow over a round cylinder and a flat plate accurately first.

[LAH]

As suggested, I simulated flow over a circular cylinder at the same Reynolds number. For this I used the exact same 'meshing strategy' as for the square cylinder with even the exact same number of elements in the mesh. The only difference is that where I refined cells in x+ direction towards the corners of the square cylinder, I could not do that for the circular cylinder (did result in lower CFL numbers though). The cylinder diameter is 0.2m, equal to the square cylinder side length.

Also the CFX pre setup was kept the same as for the square cylinder simulation. Velocity = 20 m/s and even though it shouldn't matter I kept the angle of 45 degrees under which the flow enters the domain.

Attached are some figures obtained from CFD post. The results look fine and as can be seen the average drag would be between 0.3 and 0.4. This seems to be in line with results I could find on the internet. The Reynolds number range around 250,000 is tricky though as this coincindes with the 'drag crisis' region where the drag coefficient rapidly drops. Maybe a different Reynolds number might be an easier test case.

Even though the low pressure on the back of the cylinder due to high velocity 'sweeps' around the corner increase the drag coefficient, I am starting to consider the possibility of this actually being physical (or at least something CFD is supposed to produce). I tried using the same mesh for a simulation in both fluent and openfoam which also showed this low pressure near the back. Unfortunately I could not fully complete these simulations however it was clear that the transient simulations in both fluent and openfoam were heading towards similar Cd values as I found using CFX.

As mentioned before, I will be reproducing the CFD setup in a windtunnel in a week, to perform experiments and measure surface pressures and forces acting on the square cylinder. Hopefully this will clarify how accurate the CFD results are.

I will try to simulate flow around a flat plate oriented normal to the freestream as a validation case. If there is still something I haven't tried which could help, please let me know, the support is very much appreciated.

[Opaque]

https://ntrs.nasa.gov/api/citations/...9840012155.pdf

https://en.wikipedia.org/wiki/Strouhal_number



http://www.dept.aoe.vt.edu/~aborgolt...pt3/index.html

[ghorrocks]

Can you try a different meshing strategy? The diagonal approach you are currently using looks like it will result in poor mesh quality at the corners, and that is where the separations are occurring. For a square cylinder you should be able to do a mesh which is completely orthogonal, so all mesh is 90°. You can refine close to the square surface of course, but keep everything orthogonal.