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.

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

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

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 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.

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.

Post Processing

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

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}