Hi folks,
I am simulating the incompressible flow in a cylinder. Consider the axisymmetric, I created the 2D model. In my model, the geometry is a rectangle which is the axisymmetric part of a cylinder. There are left, top, bottom, and right boundaries. The left boundary is the axisymmetric boundary where dot(u, n) and shear stress are zero. On all other boundaries u=(0,0). So sigma * n * v = 0 on ds and there is no ds term in the variational form. The x-axis in the model is actually r-axis and y-axis is the z-axis. I convert the nabla_grad and div operators to cylindrical coordinates and dv should be 2pirdrdz=x[0]dx as the 2*pi are eliminated in the variational form.
The model cannot converge after a few time steps as the solved velocity is super high (~10e103). I tried a lot of different initial conditions but they were not helpful. Any suggestions would be appreciated.
My DOLFIN version is 2018.1.0 and the error message says:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** fenics-support@googlegroups.com
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the 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: 2018.1.0
*** Git changeset:
*** -------------------------------------------------------------------------
My code is:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Apr 21 10:50:16 2023
@author: xin
"""
from fenics import *
from mshr import *
import matplotlib.pyplot as plt
import numpy as np
tt = 0.01 # Simulation time
num_steps = 10 # Number of time steps
dt = tt/num_steps # Time step size
k = Constant(dt)
Lv = 0.4 # Length of the vaporizer 0.4m
La = 1.2 # Length of the adiabatic section 1.2m
Lc = 0.4 # Length of the condenser 0.4m
tot_l = Lv+La+Lc # Total length
h = 0#0.02 # Thickness of the tube
r = 0.1 # Radius of heat pipe
# Define material properties
# viscosity
mu = 1.34*10**(-5)
rho = 0.5542 # Density of vapor
# Create domain
whole_domain = Rectangle(Point(0,0),Point(h+r,tot_l+2*h)) # Whole domain
tol = 1E-14 # Tolorence used to distinguish different domain and define boundaries
# Define boundaries
class top(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], tot_l + 2*h)
class bot(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0)
class left(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0],0)
# Right boundary of the condenser
class right_con(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], r+h) and x[1]>= tot_l +h - Lc
# Right boundary of the adiabatic section
class right_adi(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], r+h) and x[1]<= tot_l -Lc + h and x[1]>=Lv + h
# Right boundary of the vaporizer
class right_vap(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], r+h) and x[1]<=Lv + h
top = top()
bot = bot()
left = left()
right_vap = right_vap()
right_adi = right_adi()
right_con = right_con()
# Generate Mesh
mesh = generate_mesh(whole_domain, 150)
x = SpatialCoordinate(mesh)
# normal vector
n = FacetNormal(mesh)
# Mark the tube domain as 0 and the fluid domain as 1
# material = MeshFunction('size_t', mesh, 2, mesh.domains())
# vtkfile = File('meshf.pvd')
# vtkfile<<material
# Mark the boundary
bf = MeshFunction('size_t', mesh,1)
bf.set_all(0)
top.mark(bf, 1)
bot.mark(bf, 2)
left.mark(bf, 3)
right_con.mark(bf, 4)
right_adi.mark(bf, 5)
right_vap.mark(bf, 6)
vtkfile = File('meshf_f.pvd')
vtkfile<<bf
# Create function space
V = VectorFunctionSpace(mesh, 'P', 2) # for velosity
Q = FunctionSpace(mesh, 'CG', 1) # for pressure, phase field, and temp
# Define test and trial functions
# For velocity
u = TrialFunction(V)
v = TestFunction(V)
# For Pressure
P = TrialFunction(Q)
q = TestFunction(Q)
#Define functions
u_n = Function(V) # save the result at n time step
u_ = Function(V) # save the tentative result
P_n = Function(Q)
P_ = Function(Q)
# define operations in cylindrical coordinates
# define gradient in cylindrical coordintes
def nabla_grad_vector_innertensor(v):
return as_tensor([[v[0].dx(0),0, v[0].dx(1)],
[0,v[0]/x[0] ,0],
[v[1].dx(0),0, v[1].dx(1)]])
# nabala grad for inner dot of tensor and vector
def nabla_grad_vector_innervector(v):
return as_tensor([[v[0].dx(0), v[0].dx(1)],
[v[1].dx(0), v[1].dx(1)]])
def nabla_grad_scalar(v):
return as_vector((v.dx(0),v.dx(1)))
def div_vector(v):
return (1/x[0])*(x[0]*v[0]).dx(0) + v[1].dx(1)
# create XDMF files for visualization output
file_u = File('heatpipe/velocity.pvd')
file_p = File('heatpipe/pressure.pvd')
# initialize the velocity field
u_n = interpolate(Expression(('0','6.0*x[0]*(0.41-x[0])/pow(0.41,2)'), degree=2), V)
# define the variational form for incompressible flow (splitting method)
#define body force
# gravity density
f = Constant((0,-9.8))
U = 1/2*(u_n+u)
# define bc for velosity
bcu_top = DirichletBC(V, Constant((0,0)),top)
bcu_bot = DirichletBC(V, Constant((0,0)),bot)
bcu_left = DirichletBC(V.sub(0), Constant(0),left)
bcu_right_con = DirichletBC(V, Constant((0,0)),right_con)
bcu_right_adi = DirichletBC(V, Constant((0,0)),right_adi)
bcu_right_vap = DirichletBC(V, Constant((0,0)),right_vap)
bcu = [bcu_top,bcu_left,bcu_bot,bcu_right_con,bcu_right_adi,bcu_right_vap]
# Define strain tensor
def epsilon(u):
return sym(nabla_grad_vector_innertensor(u))
# Define Stress tensor
def sigma(u,P):
return 2*mu*epsilon(u) - P*Identity(3)
# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*x[0]*dx \
+ rho*dot(dot(u_n, nabla_grad_vector_innervector(u_n)), v)*x[0]*dx \
+ inner(sigma(U, P_n), epsilon(v))*x[0]*dx \
- dot(f, v)*x[0]*dx
a1 = lhs(F1)
L1 = rhs(F1)
# Define variational problem for step 2
a2 = dot(nabla_grad_scalar(P), nabla_grad_scalar(q))*x[0]*dx
L2 = dot(nabla_grad_scalar(P_n), nabla_grad_scalar(q))*x[0]*dx - (1/k)*div_vector(u_)*q*x[0]*dx
# Define variational problem for step 3
a3 = dot(u, v)*x[0]*dx
L3 = dot(u_, v)*x[0]*dx - k*dot(nabla_grad_scalar(P_ - P_n), v)*x[0]*dx
# Assemble matrices
# since there is no parameters or function need to be updated at the left hand side, so we assemble it once
# step 1: solve the tentative velocity
A1 = assemble(a1) # pho need to be updated after each time step
[bc.apply(A1) for bc in bcu]
A2 = assemble(a2)
A3 = assemble(a3)
# Use amg preconditioner if available
prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"
# time steping
t = 0
for n in range(num_steps):
# for n in range(1):
print(n)
t = t+dt
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)
solve(A2, P_.vector(),b2,'cg',prec)
# step 3: velocity correction step
b3 = assemble(L3)
solve(A3, u_.vector(),b3,'gmres','default')
# update previous solution
u_n.assign(u_)
P_n.assign(P_)
# save solution to files
file_u<<(u_,t)
file_p<<(P_,t)