Hi everyone.
I am trying to simulate an imcompressible flow through an axisymmetric venturi tube.
I have implemented the equations in FeniCSx. When I solve the equations the velocity looks good but the pressure does not. if I demand 0 pressure at the outlet of the pipe, I get
If I change it to pressure 1 at the outlet of the pipe
On the other hand, the pressure should increase after convergence. I have tested for different numbers of reynolds and the behavior is the same.
Anyone has a idea?
This is my full code:
import dolfinx
import gmsh
from mpi4py import MPI
from petsc4py import PETSc
import pyvista
import numpy as np
import math
import ufl
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
as_vector,as_tensor, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym, MixedElement)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
create_vector, create_matrix, set_bc)
import matplotlib.pyplot as plt
from dolfinx.fem import (Constant, Function, FunctionSpace,
assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc)
from dolfinx.mesh import create_unit_square
from petsc4py.PETSc import ScalarType
from dolfinx import geometry
from dolfinx.io import gmshio
###############################################################################
#Mesh with GMsh
################################################################################
gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0
gmsh.initialize()
#Points
n_point = 50
n_point1 = int(4*n_point)
n_point2 = int(0.8*n_point)
n_point3 = int(n_point)
n_point4 = int(1.5*n_point)
n_point5 = int(10*n_point)
n_point_axis = n_point1 + n_point2 + n_point3 + n_point4 + n_point5
bump = 0.05
a = 2.0*(bump-1)/(n_point*1.0)+1
n_bump = int(math.log(bump,a))
c = 1e-1;
gmsh.model.geo.add_point(0, 0.025, 0, c,1)
gmsh.model.geo.add_point(0.225, 0.025, 0, c,2)
gmsh.model.geo.add_point(0.252, 0.02, 0, c,3)
gmsh.model.geo.add_point(0.292, 0.02, 0, c, 4)
gmsh.model.geo.add_point(0.3737, 0.025, 0, c, 5)
gmsh.model.geo.add_point(0.945, 0.025, 0, c, 6)
gmsh.model.geo.add_point(0.945, 0, 0, c, 7)
gmsh.model.geo.add_point(0, 0, 0, c,8)
gmsh.model.occ.synchronize()
#Lines
gmsh.model.geo.addLine(1, 2 ,1)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4,3)
gmsh.model.geo.addLine(4, 5,4)
gmsh.model.geo.addLine(5, 6,5)
gmsh.model.geo.addLine(6, 7,6)
gmsh.model.geo.addLine(7, 8,7)
gmsh.model.geo.addLine(8, 1,8)
gmsh.model.occ.synchronize()
gmsh.model.geo.add_curve_loop([1, 2, 3, 4, 5, 6, 7, 8],1)
gmsh.model.geo.add_plane_surface([1],1)
gmsh.model.geo.synchronize()
volumes = gmsh.model.getEntities(dim=gdim)
gmsh.model.addPhysicalGroup(1, [8], name='inflow')
gmsh.model.addPhysicalGroup(1, [6], name = 'outflow')
gmsh.model.addPhysicalGroup(1, [1,2,3,4,5], name= 'walls')
gmsh.model.addPhysicalGroup(1, [7], name= 'axis')
gmsh.model.addPhysicalGroup(2, [1], 1)
gmsh.model.occ.synchronize()
gmsh.model.geo.mesh.setTransfiniteCurve(6, n_bump, "Progression", 1.05)
gmsh.model.geo.mesh.setTransfiniteCurve(8, n_bump, "Progression", -1.05)
gmsh.model.geo.mesh.setTransfiniteCurve(1, n_point1)
gmsh.model.geo.mesh.setTransfiniteCurve(2, n_point2)
gmsh.model.geo.mesh.setTransfiniteCurve(3, n_point3)
gmsh.model.geo.mesh.setTransfiniteCurve(4, n_point4)
gmsh.model.geo.mesh.setTransfiniteCurve(5, n_point5)
gmsh.model.geo.mesh.setTransfiniteCurve(7, n_point_axis-4)
gmsh.model.geo.mesh.setTransfiniteSurface(1, "Left", [1,6,7,8])
gmsh.model.geo.mesh.setRecombine(2, 1)
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(2)
domain, _, ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
ft.name = "Facet markers"
gmsh.finalize()
########################################################################################################
# Definition of Spaces. (dim and polinomial type)
v_cg1 = VectorElement("CG", domain.ufl_cell(), 2)
s_cg1 = FiniteElement("CG", domain.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(domain, v_cg1)
Q = dolfinx.fem.FunctionSpace(domain, s_cg1)
Z = MixedElement([v_cg1, s_cg1])
W = dolfinx.fem.FunctionSpace(domain, Z)
# Boundary subdomains and coditions
########################################################################################################
# Definition of subdomains
Re = 500
Lx =0.945
def upstream(x):
return np.isclose(x[0],0)
def downstream(x):
return np.isclose(x[0],Lx)
def walls(x):
a = np.greater(x[0],0)
b = np.less(x[0],Lx)
return np.logical_and(np.greater(x[1]**2, 0.0000001), np.logical_and(a,b))
def axis(x):
return np.isclose(x[1],0)
fdim = domain.topology.dim -1
up_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, upstream)
down_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, downstream)
wall_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, walls)
axis_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, axis)
def u_exact(x):
values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = 1
return values
def p_exact(x):
values = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)
values[:] = 0
return values
V1, _ = W.sub(0).collapse()
Q1, _ = W.sub(1).collapse()
u_icondition = Function(V1)
u_wallcondition = Function(V1)
p_icondition = Function(Q1)
u_icondition.interpolate(u_exact)
p_icondition.interpolate(p_exact)
up_f = dolfinx.fem.locate_dofs_topological((W.sub(0),V1), fdim, up_facets)
wall_f = dolfinx.fem.locate_dofs_topological((W.sub(0),V1), fdim, wall_facets)
down_f = dolfinx.fem.locate_dofs_topological((W.sub(1),Q1), fdim, down_facets)
boundaries = [(1,up_facets),(2,down_facets),(3,wall_facets), (4,axis_facets)]
facet_indices, facet_markers = [], []
for (marker, face ) in boundaries:
facet_indices.append(face)
facet_markers.append(np.full_like(face, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
#############################################################################################################
# Define ds measure
ds = Measure("ds", domain=domain, subdomain_data=facet_tag)
##############################################################################################################
# Solve Navier-Stokes equtions axis
n = FacetNormal(domain)
radius = Function(Q1)
radius.interpolate(lambda x: x[1])
# Definiton of stress tensor in cylindrical coordinates
def epsilon(u):
return sym(as_tensor([[u[0].dx(0),0,u[0].dx(1)],
[0,u[1]/radius,0],
[u[1].dx(0),0,u[1].dx(1)]]))
# Define stress tensor
def sigma(u, p):
return 2*mu*epsilon(u) - p*Identity(3)
bc_u = dirichletbc(u_icondition, up_f,W.sub(0))
bc_w = dirichletbc(u_wallcondition, wall_f,W.sub(0))
bc_p = dirichletbc(p_icondition, down_f,W.sub(1))
bc = [bc_u, bc_w, bc_p]
mu = Constant(domain, PETSc.ScalarType(1/Re))
# Define variables
w = Function(W)
u, p = ufl.split(w)
wtest = TestFunction(W)
v,q = ufl.split(wtest)
f = Function(Q1)
# rdrdz = dV in axisymmetric
a1 = dot(dot(u,nabla_grad(u)),v)*radius*dx + inner(sigma(u,p),epsilon(v))*radius*dx
a2 = dot(div(u)*radius + u[1],q)*dx
#a2 = dot(div(ua)-f,qa)*dx
a3 = dot(p*n,v)*ds(2) - mu*dot(dot(grad(u),n),v)*ds(2)
#+ dot(dot(u,n),q)*ds(4)
a0 = a1+a2+a3
problem_a = dolfinx.fem.petsc.NonlinearProblem(a0, w, bcs=bc)
solver_a = dolfinx.nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem_a)
#solver.convergence_criterion = "incremental"
solver_a.rtol = 1e-4
solver_a.report = True
dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
n_it, converged = solver_a.solve(w)
u,p = w.split()
with dolfinx.io.XDMFFile(domain.comm, "Plot_axis/u.xdmf", "w") as file:
file.write_mesh(domain)
file.write_function(u)
with dolfinx.io.XDMFFile(domain.comm, "Plot_axis/p.xdmf", "w") as file:
file.write_mesh(domain)
file.write_function(p)
Thank you much!