Skip to content

Cylinder in a Crossflow


This problem considers a 2D cylindrical obstacle (diameter d=0.02 m) within an incompressible flow. Inflow velocity U_0 = 10 m/s, kinematic viscosity \nu = 0.001 m^2 /s, Reynolds number \mathrm{Re} \approx 200. More information can be found here.

CFD Setup

Air at 20^\circC should be defined as material for the simulation (Material parameters are inside the air.xml file.). The speed of sound of air can be assumed to be c_0 \approx 343.5 \frac{m}{s}. The computational aeroacoustic domain is depicted below.

CAA Setup

Inside the source region, a transient incompressible CFD Simulation is performed. Thereafter, aeroacoustic sourceterms can be computed from the obtained flow field. We then solve the Perturbed Convective Wave Equation (PCWE)

\frac{1}{c^2} \, \, \frac{\mathrm{D}^2\psi^{\rm a}}{\mathrm{D} t^2} - \nabla \cdot \nabla \psi^{\rm a} = - \frac{1}{\bar \rho c^2}\, \frac{\mathrm{D} p^{\rm ic}}{\mathrm{D} t}

for this low Mach number \mathrm{Ma} \approx \frac{U_0}{c_0} = 0.03 \ll 0.3 flow. The convective wave equation fully describes acoustic sound generated by incompressible flow structures and its wave propagation. The only unkown is the acoustic scalar potential \psi^{\rm a} . In order to receive the acoustic pressure, we have to derive the scalar potential and scale it with the density

p^\mathrm{a} = \rho_0 \frac{\mathrm{D} \psi^\mathrm{a}}{\mathrm{D} t} \,.

For this problem, the convective part of the operator \frac{\mathrm{D}}{\mathrm{D} t} = \frac{\partial}{\partial t} + \boldsymbol{\bar{v}} \cdot \nabla has a negligible effect on the acoustic field. Therefore, the convective part of the RHS source term \boldsymbol{v} \cdot \nabla p^{\rm ic} is neglected for simplicity.

Here you can download the full simulation setup. The convective part is taken into account in the file_PCWE files and neglected in the file_dpdt files.

Perform Transient Flow Simulation with openFOAM (openFoam v11)

In older openFoam versions (before v10), the file physicalProperties was named transportProperties

Based on the Strouhal number for flows around obstacles of \mathrm{St} \approx 0.2 we estimate a vortex shedding frequency of

f = \mathrm{St} \frac{U_0}{d} \approx 100 \,\mathrm{Hz}

Therefore we choose a timestep of \Delta t_\mathrm{CA} = 5 \cdot 10^{-4} \,\mathrm{s} \approx \frac{1}{20 f}.

  • Set up the OpenFoam Mesh
  • Define Boundary Conditions
  • Velocity inlet
  • Pressure outlet
  • Specify transient controls and solution methods
  • Time stepping based on \mathrm{CFL} < 1: \Delta t_{\mathrm{CFD}} = 5\cdot 10^{-5} \,\mathrm{s}
  • Standard solution methods
  • Incompressible transient flow solver: icoFoam using the PISO algorithm
  • Run simulation until steady oscillation is set up
  • Export pressure data as EnSight Gold Case (*.case) foamToEnsight -fields '( "p.*" )' -time 0.5:0.7
  • Standard CFD-Result file format
  • Can be exported with OpenFoam, StarCCM++, Ansys Fluent, etc.
  • Can be visulized with ParaView

The commands to be executed are on a Linux machine:

  • bash (this command initializes the run based on a precomputation stored at 0.5s, essentially sets up the MPI-run, you will get for each processor a file)
  • bash (is the actual computation, stores the data in Foam format for the requested writing timesteps and the transform the data to the Ensight format) - this computation step takes about 20 minutes

Interpolate the Flow Field and Compute Acoustic Sources

Based on the Perturbed Convective Wave Equation (PCWE) the RHS source needs to be caluclated. For this example, the convective part of the RHS source is neglected and only \frac{\partial p^\mathrm{ic}}{\partial t} needs to be calculated.

  • Define the variable $CFS_DIR pointing to your cfs version (e.g. add the line export CFS_DIR="PATH_TO_CFS_INSTALL_ROOT_DIR" to your ~/.bashrc on Linux)

    It is necessary to use a maximum CFS build (which is GPL v3 licensed) that includes CGAL, which is not included in the standard release, and has to be built manually (Instructions for building can be found here).

  • Define data input
  • Define interpolation filter
  • Define partial time derivative
  • Define data output
  • Calculate interpolated acoustic sources: $ cfsdat -p interpolation.xml -tX interpolation using X threads in parallel. (takes about 20 seconds)

    The calculation scripts for the partial time derivative of the pressure and the PCWE source are provided in the bash scripts and

Acoustic Simulation Setup

Finally, the acoustic propagation simulation is performed using the (already integrated) RHS values of

- \frac{1}{\rho_0 c_0^2} \int_\Omega \frac{\partial p^{\rm ic}}{\partial t} \, \, d \Omega \, .
  • Define I/O files
  • Define data and mesh input
  • Include material file (air.xml)
  • Define output file formats
  • Define computational domain
  • Assign material to regions (air at 20 degrees)
  • Define non-conforming interfaces between the different regions/parts (Nitsche-type Mortaring)
  • Define microphone points for evaluation
  • Define simulation
  • Specify simulation type (transient)
  • Define the acoustic PDE (Potential formulation)
    • Specify regions to use in calculation
    • Specify used NC interfaces
    • Specify damping for perfectly-matched-layer (PML) regions
    • Specify background flow for convective wave operator
    • Define RHS node values (multiply with - \frac{1}{\rho_0 c_0^2})
    • Define output results
  • Run the acoustic propagation simulation: $ cfs -p propagation.xml -tX propagation using X threads in parallel. (takes about 4 minutes)

    The calculation scripts for the partial time derivative of the pressure and the PCWE source are provided in the bash scripts and

Post Processing

  • View acoustic field result with ParaView
  • Further process point data results from history and sensorArray

    The postprocessing scripts for the partial time derivative of the pressure and the PCWE source are provided in the bash script Therefore a Python 3.X installation with the packages: numpy, scipy, and matplotlib is required.

Time signal of the acoustic pressure\label{fig:res1}

Power spectral density of the acoustic pressure\label{fig:res2}

The figure above shows the sound pressure level (SPL) of the signals in 2 microphones located downstream (microphone 1) and to the side of the cylinder (microphone 2). The SPL is calculated by

L_{\mathrm{p}_{\mathrm{a}}}=20 \log _{10} \frac{p_{\mathrm{a}, \mathrm{rms}}}{p_{\mathrm{a}, \mathrm{ref}}}

with a reference pressure p_{\mathrm{a}, \mathrm{ref}}=20 \mu \mathrm{Pa} and the RMS can be calculated from the amplitude spectrum by p_{\mathrm{a}, \mathrm{rms}} = \frac{\hat{p}_{\mathrm{a}}}{\sqrt{2}}. The amplitude spectrum is calculated as the square root of the power spectrum obtained by Welch's method. The spectrum of the acoustic pressure signal in microphone 1 shows a strong peak for the 2nd harmonic of the vortex shedding frequency, while in microphone 2 the vortex shedding frequency itself is dominant.

This behaviour can also be observed looking at the directivity plot. At 95 Hz (about vortex shedding frequency) and 190 Hz (about double of the vortex shedding frequency) the source show a dipole behaviour, facing cross-stream and in-stream respectively.

Directivity plot at 95 and 190 Hz\label{fig:res3}

The RMS of the normalized projection of the acoustic intensity

I_{r a d, i}=\frac{I_{x, i} x_i+I_{y, i} y_i}{\sqrt{x_i^2+y_i^2}}

at each microphone point gives a dipole behavior of the acoustic intensity.

Directivity plot at 95 and 190 Hz\label{fig:res3}