Example 4 : 3D Flow around a Sphere

In this example, a flow is passing around a fixed sphere. The velocity profile is simulated at two Reynolds numbers: Re = 0.1 and Re = 150 using parameter files sphere_0.1.prm and sphere_150.prm. Later, a dynamic adaptive mesh is introduced, which can be simulated using sphere_adapt.prm This example can be solved using the gls_navier_stokes_3d solver.

The following schematic describes the simulation.

• bc = 0 (No slip boundary condition)
• bc = 1 (u=1; flow in the x direction)
• bc = 2 (Slip boundary condition)

This example can be seen as a 3D version of Example 2: 2D Flow around a cylinder, and hence the discussion concerning the boundary conditions also applies here.

The structured mesh is built using gmsh. Geometry parameters can be adapted in the "Variables" section of the .geo file, as shown below. The domain is sized so that the wake has sufficient distance to develop downstream and there is a sufficient cross-sectional area to negate any effect from the wall boundary conditions.

// Variables
gr=4; 			// Global refinement
around=10;		// Refinement around the sphere
trail=8;		// Refinement of trail of sphere
near_sphere=1.20;	// Progression of cell refinement to sphere surface
downstream=16;	        // Length of domain downstream of sphere
upstream=-6;		// Length of domain upstream of sphere (must be negative)
cross_section=10;	// Half the cross-sectional width/height
radius=0.5;		// Radius of sphere

The Reynolds number is varied by maintaining a constant inlet velocity of u=1 and varying the kinematic viscosity of the flow.

#---------------------------------------------------
# Physical Properties
#---------------------------------------------------
subsection physical properties
    set kinematic viscosity            = 0.01
end

The force acting on the sphere is also calculated.

#---------------------------------------------------
# Force
#---------------------------------------------------
subsection forces
    set verbosity             = verbose
    set calculate forces      = true
    set calculate torques     = false
    set force name            = force
    set output precision      = 10
    set calculation frequency = 1
    set output frequency      = 1
end

Re = 0.1

At Re=0.1, the flow is attached to the sphere and there is no recirculation behind the sphere. The flow regime is laminar and steady, and the wake is weak, axisymmetrical, and stable. The flow can therefore be easily solved for steady-state.

# --------------------------------------------------
# Simulation Control
#---------------------------------------------------
subsection simulation control
  set method                  = steady
  set number mesh adapt       = 0
  set output name             = sphere-output
  set output frequency        = 1
  set subdivision             = 1
end

The results are visualised by taking a z-normal slice through the centre of the mesh. The velocity is illustrated in the following image:

The pressure is also visualised:

The drag on the sphere is available in the output file force.00.dat (the other force files force.01.dat and force.02.dat give the forces on boundary conditions 1 and 2 respectively). The last line of the file shows the force calculated in the last iteration. Since the flow in the x-direction, the x-direction force f_x gives the drag force.

cells      f_x           f_y          f_z      
 5823 98.3705224612 -0.0000000785 0.0000001119

The drag coefficient of the sphere is calculated from the drag force using the equation: , where is the force obtained from force.00.dat, u is the inlet velocity and A is the cross-sectional area of the sphere.

At Re=0.1, an analytical solution of is known: = 240. The drag coefficient calculated in this example by Lethe is 250.50. Using only 6000 cells means that there is a low spatial resolution; hence, the deviation from the analytical solution is primarily due to the coarseness of the mesh.

Re = 150

At Re=150, the flow has separated, resulting in an unstable wake and recirculation. It is hence more difficult to converge to a steady-state solution. Therefore, two changes are implemented to allow a steady-state solution to be found. Firstly, the time-stepping method is changed from steady to steady_bdf.

#--------------------------------------------------
# Simulation Control
#---------------------------------------------------
subsection simulation control
  set method                  = steady_bdf
  set time step        	      = 0.1
  set adapt 		      = true
  set max cfl		      = 1000
  set stop tolerance          = 1e-5
  set adaptative time step scaling = 1.2
  set number mesh adapt       = 0
  set output name             = sphere-output
  set output frequency        = 1
  set subdivision             = 1
end

The steady_bdf method solves for a steady-state simulation using adjoint time stepping with a bdf1 scheme. An initial time step is used to complete a transient iteration, and with each iteration, the time step is increased. The simulation is considered to have reached steady-state when the L2 norm of the initial residual is lower than stop tolerance at the start of a non-linear solution step, i.e. until the time step is large enough that a pseudo-steady-state has been reached.

Secondly, an initial condition is introduced. Since the Reynolds number is varied by varying the kinematic viscosity, a viscous initial condition is set. Since the solution can easily be found at Re=10, this is used as an initial attempt to hence find the solution at Re=150:

#---------------------------------------------------
# Initial condition
#---------------------------------------------------
subsection initial conditions
    set type = viscous
    set viscosity = 0.1
end

The velocity and pressure are once again visualised:

The drag coefficient at Re=150 using this example simulation is 0.858. The coarseness of the grid can clearly be seen in the lack of clarity in the velocity profile near the sphere, and so refinement of the mesh must occur to gain a more accurate simulation.

Addition of a dynamic adaptive mesh

To increase the accuracy of the drag coefficient, the mesh must be refined in areas of interest, such as on the front face of the sphere and in the developing wake. Therefore, a dynamic adaptive mesh was introduced to refine the mesh in such regions.

# --------------------------------------------------
# Mesh Adaptation Control
#---------------------------------------------------
subsection mesh adaptation
  set type                    = kelly
  set fraction coarsening     = 0.1
  set fraction refinement     = 0.3
  set fraction type	      = number
  set max number elements     = 50000
  set min refinement level    = 0
  set max refinement level    = 3
  set variable		      = velocity
end

The mesh is dynamically adapted based on an estimate of the error of the solution for the velocity (the Kelly error estimator). The refinement is based on the number of elements. This means that the number of cells refined/coarsened per iteration is based on the fraction of the number of cells, rather than the fraction of the error (where all cells which have the fraction of the error are refined/coarsened).

The min refinement level refers to the base mesh which has been used in the previous static simulations. The mesh can only become finer than this, not coarser. The max refinement level is set at 3, giving a maximum possible number of cells at 3 million. However, the max number of elements limits the number of cells to 50,000 to keep the simulation within feasible computational expense.

The resulting velocity profile is shown without and with the underlying mesh. Refinement around the sphere and wake can be observed:

The resulting drag coefficient of 0.855 is more accurate than determined using the static mesh. To increase accuracy further, the max number of elements and max refinement level of the mesh adaption can be increased.