Skip to content

Taylor Couette flow using Nitsche IB

Bruno Blais edited this page Apr 28, 2021 · 11 revisions

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

Uniform mesh refinement

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.

Adaptive mesh refinement

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:

Clone this wiki locally