Taylor Couette flow using Nitsche IB
In this example, we revisit the flow that we investigated in the Taylor‐Couette flow example.
The goal of this example is to demonstrate some of the capabilities of Lethe to simulate the flow around complex geometries without meshing them explicitly with a conformal mesh, but instead by using the Nitsche immersed boundary method available within Lethe.
We simulate the same case as the regular Taylor-Couette flow where the inner cylinder rotates at a constant velocity '''Omega''', while the outer cylinder is fixed. The following figure shows the geometry of this problem and the corresponding boundary conditions. Note that the outer cylinder does not rotate, while the inner cylinder rotates anti-clockwise.
The inner cylinder is represented by an immersed boundary. Consequently, the background mesh is a circle with no holes for the inner cylinder:
subsection mesh
set type = dealii
set grid type = hyper_ball
set grid arguments = 0, 0 : 1 : true
set initial refinement = 3
end
The Nitsche Immersed Boundary (IB) uses particles located at the Gauss quadrature points of the immersed mesh to represent the immersed body. For a thorough explanation of this, we refer the reader to step-70 of deal.II.
The Nitsche IB is parametrized using the following section in the parameter file:
# --------------------------------------------------
# Nitsche
#---------------------------------------------------
subsection nitsche
set beta = 10
set number of solids = 1
subsection nitsche solid 0
subsection mesh
set type = dealii
set grid type = hyper_ball
set grid arguments = 0, 0 : 0.25 : true
set initial refinement = 6
end
subsection solid velocity
set Function expression = -y ; x
end
set enable particles motion = false
end
set calculate torques on solid = true
set verbosity = verbose
end
The beta coefficient is a parameter used to enforce the Nitsche IB. It's value is generally between 1 and 100. A value of 10 is reasonable. Then, we specify the number of Nitsche IB. For each Nitsche IB, a mesh representing the immersed solid has to be specified. Additionally, the velocity of the Nitsche IB is specified using the solid velocity subsection. Finally, the motion of the particle is disabled. This means that even if the IB has a non-zero velocity, it will not be moved. In this case, this is because our problem has rotation symmetry and we will be seeking steady-state solutions. We note that in this problem, the Nitsche solid grid has the same dimension as the background grid. This is necessary for 2D simulations. Additionally, the Nitsche solid grid is well-refined to ensure that at approximately each fluid cell contains one particle of the immersed body. Finally, we enable the calculation of the torque on the Nitsche IB by setting calculate torsque on solid to true.
The following figure illustrates the background mesh as well as the particles used to represent the IB on top of it:
Like in the first Taylor-Couette example, we add an analytical solution section to the parameter handler file. This analytical solution is more complex to define, since the simulation domain encompasses the inside of the inner cylinder as well as the gap between the cylinders. Because of this, we only specify the analytical solution for the velocity field and forego pressure. The analytical solution is only defined in the .prm file and we do not reproduce it here for the sake of brevity.
In this problem, we use first order elements for velocity and first order elements for pressure (Q1-Q1).
#---------------------------------------------------
# FEM
#---------------------------------------------------
subsection FEM
set velocity order = 1
set pressure order = 1
set qmapping all = true
end
As we will see further down the line, the Nitsche IB method reduces the global order of convergence to 1 since the boundary is only represented approximately.
Even if we output the torque on the Nitsche solid, we are still interested in calculating it on the outer cylinder. Consequently, we set the forces subsection accordingly:
#---------------------------------------------------
# Force
#---------------------------------------------------
subsection forces
set verbosity = verbose # Output force and torques in log
set calculate torques = true # Enable torque calculation
set torque name = torque # Name prefix of torque files
set calculation frequency = 1 # Frequency of the force calculation
set output frequency = 1 # Frequency of file update
end
To set up the boundary conditions, no slip condition is used for the outer cylinder.
# --------------------------------------------------
# Boundary Conditions
#---------------------------------------------------
subsection boundary conditions
set number = 1
subsection bc 0
set type = noslip
end
end
We first solve this problem on a series of 5 consequently uniformly refined meshes using the gls_nitsche_navier_stokes_22 solver.
# --------------------------------------------------
# Simulation Control
#---------------------------------------------------
subsection simulation control
set method = steady
set number mesh adapt = 4
set output name = taylor_couette_22
set output frequency = 1
set output path = ./
end
# --------------------------------------------------
# Mesh Adaptation Control
#---------------------------------------------------
subsection mesh adaptation
set type = uniform
end
The evolution of the L2 error is as follows:
320 2.6290e-02 - 1.5068e-02 -
1280 1.2266e-02 1.10 1.9538e-02 -0.37
5120 6.2622e-03 0.97 1.7759e-02 0.14
20480 3.2062e-03 0.97 1.7740e-02 0.00
81920 1.5688e-03 1.03 1.7626e-02 0.01
We discard the results for pressure since we have not specified an analytical solution. We note that as the number of cells increases, the error converges to zero at first order (error is divided by two when the mesh size decreases by a factor of two).
The torque on the inner cylinder is given in the torque_solid.00.dat file:
cells T_x T_y T_z
320 0.0000000000 0.0000000000 -0.6901522094
1280 0.0000000000 0.0000000000 -0.7673814310
5120 0.0000000000 0.0000000000 -0.8009318544
20480 0.0000000000 0.0000000000 -0.8186962282
81920 0.0000000000 0.0000000000 -0.8283917140
whereas the toque on the outer cylinder is given by the torque.00.dat file:
cells T_x T_y T_z
320 0.0000000000 0.0000000000 0.7223924685
1280 0.0000000000 0.0000000000 0.7840745866
5120 0.0000000000 0.0000000000 0.8093268556
20480 0.0000000000 0.0000000000 0.8229078025
81920 0.0000000000 0.0000000000 0.8305030116
We see that the sum of both torque converge towards zero as the mesh is refined, ensuring that Newton's third law is respected. The torque on the inner cylinder should be -0.83776 and we note that the torque on both cylinder converges close to that value. Running the simulation with finer meshes lead to this results.
Since the Nitsche IB method introduces additional error on the surface of the immersed geometry, it is pertinent to investigate the results it can produce with adaptive mesh refinement. We now consider the following options:
# --------------------------------------------------
# Simulation Control
#---------------------------------------------------
subsection simulation control
set method = steady
set number mesh adapt = 6
set output name = taylor_couette_22
set output frequency = 1
set output path = ./
end
# --------------------------------------------------
# Mesh Adaptation Control
#---------------------------------------------------
subsection mesh adaptation
set type = kelly
set variable = velocity
set fraction type = number
set max number elements = 500000
set max refinement level = 15
set min refinement level = 0
set frequency = 1
set fraction refinement = 0.3
set fraction coarsening = 0.15
end
Which produces the following output for the evolution of the L2 norm of the error on the velocity:
ells error_velocity error_pressure
320 2.6280e-02 - 1.5068e-02 -
620 1.2500e-02 1.07 1.9900e-02 -0.40
1196 6.4573e-03 0.95 1.7950e-02 0.15
2312 3.3532e-03 0.95 1.8708e-02 -0.06
4580 1.5891e-03 1.08 1.7682e-02 0.08
9056 8.4245e-04 0.92 1.7687e-02 -0.00
18284 4.3930e-04 0.94 1.8065e-02 -0.03
Interestingly, we observe an apparent second-order for the convergence of the L2 norm of the error on the velocity. This is largely in part due to the dynamic mesh adaptation and the use of the Kelly error estimator.
Correspondingly, the torque on the inner cylinder:
cells T_x T_y T_z
320 0.0000000000 0.0000000000 -0.6902135017
620 0.0000000000 0.0000000000 -0.7681788397
1196 0.0000000000 0.0000000000 -0.8020340261
2312 0.0000000000 0.0000000000 -0.8196041387
4580 0.0000000000 0.0000000000 -0.8289292869
9056 0.0000000000 0.0000000000 -0.8332883003
18284 0.0000000000 0.0000000000 -0.8353647429
We see that even for a small number of cells (~18k), the error on the torque is less than 0.5%. Finally, using paraview we can display the velocity as well as the final adapted mesh:
