Navier-Stokes Secondary Flows (Dean Drag Forces) in a curved channel

[Matt688]

Hey all,

I’ve imported a curved 3D microchannel in FEniCS and I am attempting to model secondary flow/Dean drag forces. I followed guidance from a COMSOL Demo using the initial FEniCS Navier Stokes IPCS FEniCS demo as a basis. The only differences are XDMF mesh import, custom code to extract a boundary faces’ length/width and co-ordinates as well as using ‘bicgstab’ & ‘hypre_amg’ for the solver and pre-conditioner in IPCS (see code below). I visualise the secondary flows in Paraview but the vectors are always biased towards the positive X axis (if Y-Z plane cut) or negative Y axis (if X-Z plane cut) and I am unsure why. Flow otherwise looks appropriate for the simulation.

I visualise secondary flows as follows:

1. Slice through 3D model
2. Apply “Surface Vector” filter. I ensure that vectors are parallel to the surface and I am using the velocity vector field
3. Apply glyph field onto the “Surface Vector” slice.
   (This is based on Secondary Flow Visualisation in Paraview.)

I’ve attempted different solvers/preconditioners with no change in the secondary flows. I have increased the density of the mesh as much as possible with the workstation RAM being the limitation. I am looking into whether this is an error with Paraview but if anyone has any thoughts or further questions, let me know!

Many thanks,
Matthew

"""
FEniCS tutorial demo program: Incompressible Navier-Stokes equations
for channel flow (Poisseuille) on the unit square using the
Incremental Pressure Correction Scheme (IPCS).

  u' + u . nabla(u)) - div(sigma(u, p)) = f
                                 div(u) = 0
"""

from __future__ import print_function
from fenics import *
# from dolfin import * 
import numpy as np
# from mpi4py import MPI

T = 3       # final time
num_steps = 3 #50    # number of time steps
dt = 0.000625 # time step size
mu = 0.001             # dynamic viscosity mPa/s
rho = 1000            # density
FlowVelocity = 0.016 # Flow velocity m/s
FlowVelocity = FlowVelocity * 1.6 # To average to the correcvt inflow

# Load mesh from file (NETGEN mesh as .grid to .xml using DOLFIN)
mesh = Mesh()
with XDMFFile(MPI.comm_world,"mesh_SpiralCrop.xdmf") as infile:
    infile.read(mesh)


# # Scale mesh from mm to m
scaling_factor = 0.001
x = mesh.coordinates()
x[:, :] *= scaling_factor
# Move mesh so co-ords always positive
xymin = x.min(axis=0) 
x[:, :] = x[:, :] - xymin
# Apply to mesh
mesh.bounding_box_tree().build(mesh) 

# Use nonzero guesses - essential for CG with non-symmetric BC
parameters['krylov_solver']['nonzero_initial_guess'] = True

# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 2, dim=3)#'P', 3)
Q = FunctionSpace(mesh, "CG", 1)#'P', 1)

# Define boundary conditions
dim = mesh.topology().dim()

if MPI.comm_world.Get_rank() == 0:
    print('Dim:',dim)

# Import face domains from XML (Created in NETGEN/GMSH, conv MESHIO)
mvc = MeshValueCollection("size_t", mesh, dim-1) 
with XDMFFile(MPI.comm_world,"mf_SpiralCrop.xdmf") as infile:
    infile.read(mvc, "name_to_read")
facet_domains = cpp.mesh.MeshFunctionSizet(mesh, mvc)

## Calculate limits so inflow parabolic can work on co-ordinates not at 0

# Create empty variable space
X = []
Y = []
Z = []
xInflowRange = 0
yInflowRange = 0
zInflowRange = 0
xInflowMin = 0
yInflowMin = 0
zInflowMin = 0

# if (facet_domains.array()[facet_domains.array() == 6] != []):
#     print("WARNING: Inlet face not detected, please re-assign")
#     break

# Retrive all co-ords as element for desired face
It_facet =  SubsetIterator(facet_domains,2)

# for mvc in It_facet:
#     i = i + 1
#     X1.append(mvc.midpoint().x())
#     Y1.append(mvc.midpoint().y())
#     Z1.append(mvc.midpoint().z())

#https://fenicsproject.org/qa/13995/print-coordinate-of-boundary-seperately-i-e-left-boundary/
#It_mesh = vertices([facet_domains.array() == 26])

# Collected all co-ords for desired face
for facet_domains in It_facet:
    for v in vertices(facet_domains):
        X.append(v.point().x())
        Y.append(v.point().y())
        Z.append(v.point().z())

# Ensure all processes collect co-ords for desired face
MPI.barrier(MPI.comm_world)

# Gather all co-ords to calc length/min
X = MPI.comm_world.gather(X, root=0)
Y = MPI.comm_world.gather(Y, root=0)
Z = MPI.comm_world.gather(Z, root=0)

# Sync all parallel processes for length/min calc
MPI.barrier(MPI.comm_world)

if MPI.comm_world.Get_rank() == 0:
    # Remove empty and combine all arrays
    X = np.concatenate(X)
    Y = np.concatenate(Y)
    Z = np.concatenate(Z)
    # Calculate length and min values
    xInflowRange = np.ptp(X,axis=0)
    yInflowRange = np.ptp(Y,axis=0)
    zInflowRange = np.ptp(Z,axis=0)
    xInflowMin = np.amin(X)
    yInflowMin = np.amin(Y)
    zInflowMin = np.amin(Z)

# END: Sync all parallel processes for length/min calc
MPI.barrier(MPI.comm_world)

# Broadcast all length/min calc to all nodes used
xInflowRange = MPI.comm_world.bcast(xInflowRange, 0)
yInflowRange = MPI.comm_world.bcast(yInflowRange, 0)
zInflowRange = MPI.comm_world.bcast(zInflowRange, 0)
xInflowMin = MPI.comm_world.bcast(xInflowMin, 0)
yInflowMin = MPI.comm_world.bcast(yInflowMin, 0)
zInflowMin = MPI.comm_world.bcast(zInflowMin, 0)

# Print result of length/min from one process
if MPI.comm_world.Get_rank() == 0:
    print("Point X range: ", xInflowRange)
    print("Point Y range: ", yInflowRange)
    print("Point Z range: ", zInflowRange) 
    print("Point X min: ", xInflowMin)
    print("Point Y min: ", yInflowMin)
    print("Point Z min: ", zInflowMin) 

# Sync all processes
MPI.barrier(MPI.comm_world)

# xInflowRange = np.amax(comm.gather(xInflowRange))
# yInflowRange = np.amax(comm.gather(yInflowRange))
# zInflowRange = np.amax(comm.gather(zInflowRange))
# xInflowMin = np.amin(comm.gather(xInflowMin))
# yInflowMin = np.amin(comm.gather(yInflowMin))
# zInflowMin = np.amin(comm.gather(zInflowMin))

# print("Gathered Point X range: ", xInflowRange)
# print("Gathered Point Y range: ", yInflowRange)
# print("Gathered Point Z range: ", zInflowRange) 
# print("Gathered Point X min: ", xInflowMin)
# print("Gathered Point Y min: ", yInflowMin)
# print("Gathered Point Z min: ", zInflowMin) 


# Where x[0] = X, x[1] = Y and x[2] = Z
inflow = Expression(("0.0","((((A*4.0*(x[1] - yMin)*(yRange - (x[1] - yMin))) / pow(yRange,2)) + ((A*4.0*(x[0] - xMin)*(xRange - (x[0] - xMin))) / pow(xRange, 2)))/2)","0.0"),
A=Constant(FlowVelocity), xRange=Constant(xInflowRange), yRange=Constant(zInflowRange), yMin=Constant(zInflowMin), xMin=Constant(xInflowMin), degree=dim)

facet_domains = cpp.mesh.MeshFunctionSizet(mesh, mvc)

#facet_domains = MeshFunction("size_t", mesh, dim-1, mesh.domains()) #Faces of mesh
print('facet domains:',facet_domains.array())
print('facet domain 2:',facet_domains.array()[facet_domains.array() == 2])

# We know 5 and 6 are the inlet and outlet boundaries from NETGEN
#facet_domains.array()[facet_domains.array() == 30] = 2
#facet_domains.array()[facet_domains.array() == 7] = 0 #Set all boundaries above to 0 #2 and 22
#facet_domains.array()[facet_domains.array() == 6] = 0 #Set all boundaries above to 0 #2 and 22
#facet_domains.array()[facet_domains.array() == 3] = 0

# Define boundary conditions
bcu_noslip1  = DirichletBC(V, Constant((0.0, 0.0, 0.0)), facet_domains, 3)
bcu_noslip3  = DirichletBC(V, Constant((0.0, 0.0, 0.0)), facet_domains, 4)
bcu_noslip5  = DirichletBC(V, Constant((0.0, 0.0, 0.0)), facet_domains, 5)
bcu_noslip6  = DirichletBC(V, Constant((0.0, 0.0, 0.0)), facet_domains, 6)



bcu_inflow  = DirichletBC(V, inflow, facet_domains, 2) # Constant((0,0,4.2371e-7))

bcp_outflow = DirichletBC(Q, Constant(0.0), facet_domains, 7)

bcu = [bcu_inflow, bcu_noslip1,bcu_noslip3,bcu_noslip5,bcu_noslip6]
bcp = [bcp_outflow]#bcp_inflow,

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)

# Define expressions used in variational forms
U   = 0.5*(u_n + u)
n   = FacetNormal(mesh)
f   = Constant((0, 0, 0))
k   = Constant(dt)
mu  = Constant(mu)
rho = Constant(rho)

# Define strain-rate tensor
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx + \
     rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
   + inner(sigma(U, p_n), epsilon(v))*dx \
   + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
   - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Use amg preconditioner if available
#prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"
#solver_parameters={'linear_solver':'mumps'}
#solver_parameters={'krylov_solver':'nonzero_inital_guess':True}# = True
#solver = KrylovSolver("gmres", "hypre_amg")

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Create XDMF files for visualization output
xdmffile_u = XDMFFile('navier_stokes_con/velocity.xdmf')
xdmffile_p = XDMFFile('navier_stokes_con/pressure.xdmf')

# Create time series (for use in reaction_system.py)
#timeseries_u = TimeSeries('navier_stokes_con/velocity_series')
#timeseries_p = TimeSeries('navier_stokes_con/pressure_series')

# Time-stepping
t = 0
i = 0
for n in range(num_steps):

    # Update current time
    t += dt
    i += 1

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1,'bicgstab', 'hypre_amg')#,"bicgstab","default")#"gmres", "ilu") #,"bicgstab","default")#

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2,'bicgstab', 'hypre_amg')#,"bicgstab","prec")#,"cg", prec) 'cg','petsc_amg')

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3,'cg', 'sor')#,"bicgstab","default")#'gmres','ilu')#"cg", "hypre_amg")

    # Plot solution
    #plot(u_)

    # Compute error
    u_e = Expression(('4*x[1]*(1.0 - x[1])', '0', '0'), degree=3)
    u_e = interpolate(u_e, V)
    error = np.abs(u_e.vector() - u_.vector()).max()
    print('t = %.3g: error = %.3g' % (t, error))
    print('max u:', u_.vector().max())
    print("Worst possible Courant number=",(dt*(u_.vector().max()))/mesh.hmin())

    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)

    # Save solution to file (XDMF/HDF5)
    # if (i == 1):
    xdmffile_u.write_checkpoint(u_, "velocity",t, XDMFFile.Encoding.HDF5, True)
    xdmffile_p.write_checkpoint(p_, "pressure",t, XDMFFile.Encoding.HDF5, True)


    # Save nodal values to file
    #timeseries_u.store(u_.vector(), t)
    #timeseries_p.store(p_.vector(), t)

# Hold plot
#interactive()

[dokken]

Have you tested the initial demo, and checked that your method of visualization is correct?
It is very hard to help you without a mesh.
Could you try solving your problem on a simple 3D mesh (for instance the built in unit cube mesh) to make sure you understand how to set boundary conditions?

[Matt688]

Hey dokken,

I have tested the initial 2D demo and it performs as expected. The secondary flow visualisation requires a 3D model I believe. Presented is a similar 3D model which uses the same code with minor changes where inflow occurs from top right and bottom right and reaches the centre of the ‘H’ structure. Secondary flow within the ‘H’ structure appears to be correct (see images). Secondary flows were performed using the Paraview visualisation detailed above and tested with a chemical transport equation.

I’ll have a look at built in unit cube mesh for a spiral mesh. Additionally I’ll revisit the boundary conditions as they have caused issues in the past with a different model. Currently all 6 boundaries have been defined (2 - inlet, 7 - outlet, 3-6 - walls). There is not a ‘1’ boundary face label.

Apologies for not including meshes in the initial post. Please find the “stp”, “geo”, “msh”, Meshio script and XDMF meshes available at this OneDrive link. I believe the mesh was generated with a max cell size of 0.02 (when translated to mesh units is 20 um).

Images below are:

• Mid-plane Y slice of ‘H’ channel with glyphs applied.
• Cross-section of Y slice with Z slice of ‘H’ channel with glyphs applied. Note secondary flow.
• Z slice of ‘H’ channel with glyphs applied. Note secondary flow.
• Cross-section of Y slice for the diffusion of species over time within the ‘H’ channel.
• Visualisation of diffusion within the ‘H’ channel with surface edges shown.
• Visualisation of curved microchannel surface edges (current modelled mesh).

Many thanks,
Matthew.

[dokken]

You could Also consider the Turek 3D benchmark as a reference solution: http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_flow3d.html

I have a code for a slightly newer version of dolfin, dolfinx available here: https://github.com/jorgensd/dolfinx_ipcs/tree/master
There is a mesh there, along with an IPCS solver.

[Matt688]

Hey dokken and all,

The problem is solved, see the solution below.

I’ve solved the issues which were mainly my wrong assumptions about how I visualise secondary flows and the time step required to produce these secondary flow effects.

Firstly, the reason why I had biased towards the positive X axis (if Y-Z plane cut) or negative Y axis (if X-Z plane cut) was that I took a plane which was perpendicular to the X, Y or Z axis. I did not take a plane which was perpendicular to my model and flow. This meant, if I was slightly off perpendicular to the model, I would start becoming parallel to flow and generate flow bias. The flow bias was raised in the comments section of the COMSOL Demo. Ensuring planes are taken perpendicular to the model, in my case the centre of the X axis, resulted in no flow bias (see images below).

Secondly, if the time step is too large, the simulation makes assumptions regarding the vectors and velocity through several cells of the mesh. Equally, too small of a time step, the vectors may not form correctly. These are obvious factors for simulations but they become more apparent when studying secondary flows. Thinking about how some systems decide on an appropriate time step, I came across Courant-Friedrichs-Lewy (CFL) condition. The CFL condition takes into account the velocity of the system, the smallest cell size of the mesh and a Courant number to aim for, in my case 0.5 as it had to be lower than 1 due to the explicit simulation. A CFL example was available in FEniCS in this forum which I adapted to make it velocity dependent than viscosity. The value calculated was dt = 9.73e^-5 seconds which was a little smaller dt = 6.25e^-4 seconds in the example I presented earlier.

    # Calculate dt as per CFL (with CFLmax = 0.5)
    # dt = 0.5*deltaX**2/mu - I believe is a Taylor series expansion including numerical diffusion, or numerical viscosity
    # Courant Number (<1) * smallest cell size in mesh / velocity magnitude
    dt = (0.5*deltaX)/mu 
    # dt = 1E-2*dt # Just for testing

After testing, the system produced Dean Drag forces in a stable manner.

More information about CFL can be found here. Images below show plane in the middle of the curved microchannel and of the Y-Z plane presenting Dean Drag forces. Also, I have implemented this problem within BERNAISE, a FEniCS based CFD simulator that can do phase field and electrochemistry. The example in BERNAISE is available here, with the mesh available here and a Paraview state which should load the data similar to the images presented below.