COMSOL 6.3 - Flow Past a Cylinder

The flow of fluid behind a blunt body such as an automobile is difficult to compute due to the unsteady flows. The wake behind such a body consists of unordered eddies of all sizes that create large drag on the body. In contrast, the turbulence in the thin boundary layers next to the streamlined bodies of aircraft and fish create only weak disturbances of flow.

An exception to this occurs when you place a slender body at right angles to a slow flow, because the eddies organize; a von Kármán vortex street appears with a predictable frequency and involves the shedding of eddies from alternating sides. Everyday examples of this phenomenon include singing telephone wires and an automobile radio antenna vibrating in an air stream.

From an engineering standpoint, it is important to predict the frequency of vibrations at various fluid speeds and thereby avoid undesirable resonances between the vibrations of the solid structures and the vortex shedding. To help reduce such effects, plant engineers put a spiral on the upper part of high smokestacks; the resulting variation in shape prohibits the constructive interference of the vortex elements that the structure sheds from different positions.

Model Definition

To illustrate how you can study effects of the kind described above, this model examines unsteady, incompressible flow past a long cylinder placed in a channel at right angle to the oncoming fluid. With a symmetric inlet velocity profile, the flow needs some kind of asymmetry to trigger the vortex production. This can be achieved by placing the cylinder with a small offset from the center of the flow.

The simulation time necessary for a periodic flow pattern to appear is difficult to predict. A key predictor is the Reynolds number, which is based on cylinder diameter. For low values (below 100) the flow is steady. In this simulation, the Reynolds number equals 100, which gives a developed von Kármán vortex street, but the flow still is not fully turbulent.

The frequency and amplitude of oscillations are stable features, but flow details are extremely sensitive to perturbations. To gain an appreciation for this sensitivity, you can compare flow images taken at the same time but with such minor differences as are created by different tolerances for the time stepping. It is important to note that this sensitivity is a physical reality and not simply a numerical artifact.

Before calculating the time-varying forces on the cylinder, you can validate the method of computation at a lower Reynolds number using the direct nonlinear solver. This saves time because you can find and correct simple errors and mistakes before the final time-dependent simulation, which requires considerable time.

The viscous forces on the cylinder are proportional to the gradient of the velocity field at the cylinder surface. Evaluating the velocity gradient on the boundary by directly differentiating the FEM solution is possible but not very accurate. The differentiation produces first-order polynomials when second-order elements are used for the velocity field. A far better approach is to use a pair of reaction force operators to compute the integrals of the viscous forces, comparable to second-order accurate integrals of the viscous forces. An alternative approach would be to use a pair of weak-constraint variables to enforce the no-slip condition. Preferably use the reaction force operator instead of weak constraints when computing integrals of reaction forces or fluxes in results processing.

The lift and drag forces themselves are not as interesting as the dimensionless lift and drag coefficients. These depend only on the Reynolds number and an object’s shape, not its size. They are defined as

(1)

(2)

where

•

FL and FD are the lift and drag forces

•

ρ is the fluid density

•

Umean is the mean velocity

•

A is the projected area (that is, the product of the cylinder diameter and length)

Results and Discussion

Figure 1 shows the flow pattern resulting from the geometry. The points represent massive particles that are released at the inlet boundary; these are used in an animation to further help in visualizing the flow.

Figure 1: A plot of the last time step clearly shows the von Kármán path.

The flow around a cylinder is a common benchmark test for CFD algorithms. Various research teams have tried their strengths on this problem using different techniques. Results from some of these experiments have been collected by Schäfer and Turek (Ref. 1), who also used them to compute a probable value for the “real” answer.

Figure 2 shows how the lift coefficient develops a periodic variation as the von Kármán vortex structure is formed.

Figure 2: Lift coefficient, CL, as a function of time.

Figure 3: Drag coefficient, CD, as a function of time.

Notes About the COMSOL Implementation

This model uses a Particle Tracing with Mass plot (see Particle Tracing with Mass in the COMSOL Multiphysics Reference Manual for details) in an Animation feature. More advanced functionality for simulating the trajectories of particles subjected to forces of various kinds is provided with the Particle Tracing Module.

The lift and drag coefficients (see Equation 1 and Equation 2) are computed using an Integral dataset defined on the cylinder boundaries, with the lift and drag forces, FL and FD, evaluated as the reaction-force-operator expressions -reacf(v)[N/m^2] and -reacf(u)[N/m^2], respectively (the unit brackets give the correct dimension of force per unit area). Introducing a unit channel width, W, in the out-of-plane dimension, which also defines the cylinder length, the full expressions for the integrands read

-2*reacf(v)[N/m^2]*W/(spf.rho*U_mean^2*(2*R*W)) (CL)

and

-2*reacf(u)[N/m^2]*W/(spf.rho*U_mean^2*(2*R*W)) (CD)

For more information about the use of reacf to compute reaction forces, see Computing Accurate Fluxes in the COMSOL Multiphysics Reference Manual. (An alternative is to use the components of the total traction, exterior boundaries auxiliary vector variable, -spf.T_stressy and -spf.T_stressx, defined in the Single-Phase Flow interface.)

Reference

1. M. Schäfer and S. Turek, “Benchmark Computations of Laminar Flow Around Cylinder”, E.H. Hirschel ed., Flow Simulation with High-Performance Computers II, Volume 52 of Notes on Numerical Fluid Mechanics, Vieweg, pp. 547–566, 1996.

COMSOL_Multiphysics/Fluid_Dynamics/cylinder_flow

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

Define geometry and velocity parameters, including a unit channel width in the out-of-plane dimension that is not included in this 2D model.

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
H	0.41[m]	0.41 m	Channel height
L	2.2[m]	2.2 m	Channel length
W	1[m]	1 m	Channel width (not modeled)
R	0.05[m]	0.05 m	Cylinder radius
U_mean	1[m/s]	1 m/s	Mean inflow velocity

Next, create a smoothed step function feature that you will use for ramping up the inflow velocity.

Step 1 (step1)

In the toolbar, click

and choose > .

In the window for , locate the section.

In the text field, type 0.05.

Geometry 1

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type L.

In the text field, type H.

Click

Circle 1 (c1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 4*R.

Locate the section. In the text field, type R.

Click

Difference 1 (dif1)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

Click to select the

toggle button for .

Select the object only.

In the toolbar, click

Click the

button in the toolbar.

In the window, click .

Materials

Specify the relevant fluid properties: the density and the dynamic viscosity.

Material 1 (mat1)

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	1	kg/m³	Basic
Dynamic viscosity	mu	1e-3	Pa·s	Basic

Laminar Flow (spf)

Inlet 1

In the toolbar, click

and choose .

Select Boundary 1 only.

Define a parabolic velocity profile using the predefined parameters and ramp up the velocity using the previously defined step function.

In the window for , locate the section.

In the U0 text field, type 6*U_mean*y*(H-y)/H^2*step1(t[1/s]).

By appending the inverse unit bracket [1/s] to the time variable t you get a dimensionless argument, as the step function expects.

Outlet 1

In the toolbar, click

and choose .

Select Boundary 4 only.

Mesh 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Click

If you zoom in on the inlet and the cylinder you can see the boundary layers that the physics-controlled mesh gives.

Click the

button in the toolbar.

Study 1

Step 1: Time Dependent

In the window, under click .

In the window for , locate the section.

In the text field, type range(0,0.2,3.4) range(3.5,0.02,7).

Solution 1 (sol1)

In the toolbar, click

In the window, expand the node, then click .

In the window for , click to expand the section.

From the list, choose .

In the toolbar, click

Results

Velocity (spf)

Click the

button in the toolbar.

Add a node to the first default plot group to reproduce the plot in Figure 1.

Particle Tracing with Mass 1

In the toolbar, click

and choose .

In the window for , click to expand the section.

In the text field, type 4*pi/3*1e-9.

Find the subsection. In the text field, type u.

In the text field, type v.

Locate the section. In the text field, type range(0.05,0.025,0.35).

In the text field, type 0.

Click to expand the section. From the list, choose .

Select the checkbox. In the associated text field, type 3.2.

In the text field, type 0.2.

Click to expand the section. Find the subsection. From the list, choose .

Find the subsection. From the list, choose .

Select the checkbox. In the associated text field, type 7.5.

From the list, choose .

Click .

Set the RGB values to 128, 128, and 128, respectively.

Click .

Click or on the cross-platform desktop.

Find the subsection. From the list, choose .

In the toolbar, click

Next, use the plot you just set up in an animation.

Animation 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type 50. Increasing the number of frames results in a smoother animation.

Click the

button in the toolbar.

To reproduce Figure 2 and Figure 3 of the lift and drag coefficients, first add an dataset for computing the total reaction force on the cylinder.

Integral 1

In the toolbar, click

and choose > .

Selection

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Select Boundaries 5–8 only.

Lift Coefficient

First, create a plot of the lift coefficient using the following steps:

In the toolbar, click

In the window for , type Lift Coefficient in the text field.

Locate the section. From the list, choose .

Point Graph 1

Right-click and choose .

In the window for , locate the section.

In the text field, type -2*reacf(v)[N/m^2]*W/(spf.rho*U_mean^2*(2*R*W)).

In this expression, reacf(v)[N/m^2] is the vertical component of the outward reaction force on the exterior cylinder boundaries; multiplying by the nonmodeled channel unit width W gives a dimensionless lift coefficient.

Select the checkbox. In the associated text field, type Lift coefficient.

In the toolbar, click