Skip to content
Fetching contributors…
Cannot retrieve contributors at this time
81 lines (59 sloc) 2.82 KB
// This code simulates in time the laminar, incompressible air flow around a cylinder.
// The forced input air velocity is linearly increased over time from 0 m/s to 0.15 m/s.
// Once the air velocity reaches a threshold value a von Karman vortex street appears.
// The cylinder has a diameter of 14 cm and the truncated air domain is 2 m x 0.8 m.
// The Reynolds number (rho*v*D/mu) is 1400 for the max 0.15 m/s velocity.
// - rho is the density of air [kg/m3]
// - v is the flow velocity far from the cylinder [m/s]
// - D is the cylinder diameter [m]
// - mu is the dynamic viscosity of air [Pa.s]
#include "sparselizardbase.h"
using namespace mathop;
void sparselizard(void)
// Region numbers used in this simulation as defined in the .msh file:
int fluid = 1, cylinder = 2, inlet = 3, outlet = 4;
// Load the mesh (GMSH format):
mesh mymesh("channel.msh");
// Define the cylinder skin region:
int cylskin = regionintersection({fluid,cylinder});
// Field v is the flow velocity. It uses nodal shape functions "h1" with two components in 2D.
// Field p is the relative pressure.
field v("h1xy"), p("h1");
// Force the flow velocity to 0 on the cylinder skin:
// Force a x-direction flow velocity increasing linearly over time at the inlet:
v.setconstraint(inlet, array2x1(0.01/6.0*t(),0));
// Set a 0 relative pressure at the outlet:
// Use an order 1 interpolation for p and 2 for v on the fluid region (satisfies the BB condition):
p.setorder(fluid, 1); v.setorder(fluid, 2);
// Dynamic viscosity of air [Pa.s] and density [kg/m3] at room temperature and atmospheric pressure:
double mu = 18e-6, rho = 1.2;
// Define the weak formulation for time-dependent incompressible laminar flow:
formulation laminarflow;
laminarflow += integral(fluid, predefinedlaminarflow(dof(v), tf(v), v, dof(p), tf(p), mu, rho, 0, 0, true) );
// Define the object for an implicit Euler time resolution.
// An all zero initial guess for the fields and their time derivative is set with 'vec(laminarflow)'.
impliciteuler eul(laminarflow, vec(laminarflow), vec(laminarflow));
// Set the relative tolerance on the inner nonlinear iteration:
// Run from 0sec to 90sec by steps of 0.2sec:
std::vector<vec> sols = eul.runnonlinear(0, 0.2, 90)[0];
for (int i = 0; i < sols.size(); i++)
// Transfer the solution at the ith timestep to field v:
v.setdata(fluid, sols[i]);
// Write field v with an order 2 interpolation to ParaView .vtk format:
v.write(fluid, "v" + std::to_string(1000 + i) + ".vtk", 2);
int main(void)
return 0;
You can’t perform that action at this time.