Goal: I wish to find a steady state solution for the flow through the cell. I need to reach a flowrate of about 0.1 ml/min.
Problem: The solver diverges for higher flow rates (above ~0.05 ml/min). I have tried increasing mesh resolution and shortening time step with no luck.
Questions:
- Is there a better way to define the boundary conditions on the face of a cylinder?
- For this type of problem (should be laminar and steady) what type of linear algebra solver/preconditoners should I be using? And are there any parameters I can adjust to find convergence.
Important information:
This example is an adaptation of the Navier-Stokes IPCS tutorial provided by FEniCS. Not much has changed apart from applying boundary conditions suited to my mesh, and attempting speed up computation with the choice of linear algebra backend and preconditioner.
Mesh:
Mesh of topological dimension 3 (tetrahedra) with 13359 vertices and 55927 cells. Generated with Tetgen via ‘generate_mesh(domain, 250, ‘tetgen’)’. The mesh is very coarse, but has been improved with the help of local refinement. The total volume of the cell is roughly 0.4 ml.
The Code:
A pressure gradient is set up by enforcing a Dirichlet condition at the inlet and outlet, and a noslip condition on the walls. The inlet is a circle in the xz-plane centered at the origin with radius 1.59 , that is R^2 = 2.52. I define the walls to be everything except a circle slightly smaller than the radius of the pipe R^2=2.4. I do not like this method for defining my boundary condition and I think that it could easily cause problems. If, for instance, the no slip condition includes some points on the face of the inlet then the velocity will have to jump from zero to an increasingly larger number an infinitesimal distance away.
#Parameters
P= 0.1 # Pressure (g mm^-1 s^-2)
μ= 8.94e-4 # Dynamic Viscosity of Water (g mm^-1 s^-1)
ρ= 1 # density of water (g/ml)
Flow_Rate = 0.00166 # 0.1 ml/ 60s
Radius = 1.59 # radius of inlet pipe (mm)
T = 120 # simulation time (s)
dt = 0.05 # Timestep
...
#Boundary Conditions
inflow = 'on_boundary && near(x[1],0,1e-4) && x[0]<10'
outflow = 'on_boundary && near(x[1],0,1e-4) && x[0]>10'
walls = 'on_boundary && !((((pow(x[0],2)+pow(x[2],2))<2.4) | ((pow(x[0]-19.05,2)+pow(x[2],2))<2.4)) && near(x[1],0) && on_boundary)'
bcu_noslip = DirichletBC(V, Constant((0,0,0)), walls)
bcp_inflow = DirichletBC(Q, Constant(P), inflow)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
...
#Solve
# Step 1: Tentative velocity step
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, u_.vector(), b1,"gmres", "default")
# Step 2: Pressure correction step
b2 = assemble(L2)
[bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2,"gmres", prec)
# Step 3: Velocity correction step
b3 = assemble(L3)
solve(A3, u_.vector(), b3,"gmres", "default")
...
Figure: (LEFT) Simple boundary color labeling. Red for the walls and blue the inlet.
(RIGHT) Pressure(surface colour) and Velocity (vector) solution.
Error:
*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_NANORINF, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** DOLFIN version: 2017.2.0
Fenics Version: 2017.2 (Installed via Ubuntu package)
OS: Ubuntu 16.04 LTS
