Thank you for your reply, I am just new to the fenics software, and there are still many things I don’t understand. Thanks again for your reply, but I want to know why the flow field mesh does not change in the following code:
from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
T = 0.02 # final time
num_steps = 100 # number of time steps
dt = T / num_steps # time step size
mu = 0.001 # dynamic viscosity
rho = 1 # density
# Structure sub domain
class Structure(SubDomain):
def inside(self, x, on_boundary):
return (x[0] - 0.2) ** 2 + (x[1] - 0.2) ** 2 < 0.01 + DOLFIN_EPS
mesh = RectangleMesh(Point(0, 0), Point(2.2, 0.41), 100, 200)
# Create sub domain markers and mark everaything as 0
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim())
sub_domains.set_all(0)
# Mark structure(cylinder) domain as 1
structure = Structure()
structure.mark(sub_domains, 1)
# Extract sub meshes
fluid_mesh = SubMesh(mesh, sub_domains, 0)
structure_mesh = SubMesh(mesh, sub_domains, 1)
# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)
# Define boundaries
inflow = 'near(x[0], 0)'
outflow = 'near(x[0], 2.2)'
walls = 'near(x[1], 0) || near(x[1], 0.41)'
cylinder = Structure()
# Define inflow profile
inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')
# Define boundary conditions
bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow)
bcu_walls = DirichletBC(V, Constant((0, 0)), walls)
bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_inflow, bcu_walls, bcu_cylinder]
bcp = [bcp_outflow]
# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
print(type(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))
k = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)
# Define symmetric gradient
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)
# 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('moving_rec/velocity.xdmf')
xdmffile_p = XDMFFile('moving_rec/pressure.xdmf')
# Create time series (for use in reaction_system.py)
timeseries_u = TimeSeries('moving_rec/velocity_series')
timeseries_p = TimeSeries('moving_rec/pressure_series')
t = 0
for n in range(num_steps):
# Update current time
t += dt
# Step 1: Tentative velocity step
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')
# Step 2: Pressure correction step
b2 = assemble(L2)
[bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')
# Step 3: Velocity correction step
b3 = assemble(L3)
solve(A3, u_.vector(), b3, 'cg', 'sor')
xdmffile_u.write(u_, t)
xdmffile_p.write(p_, t)
# Save nodal values to file
timeseries_u.store(u_.vector(), t)
timeseries_p.store(p_.vector(), t)
# Plot solution
plot(u_, title='Velocity')
# plt.show()
plot(p_, title='Pressure')
for x in structure_mesh.coordinates():
x[1] += 0.1 * sin(100 * dt)
# Move fluid mesh according to structure mesh
ALE.move(fluid_mesh, structure_mesh)
fluid_mesh.smooth()
plot(fluid_mesh, title="Fluid")
plt.show()
# Update previous solution
u_n.assign(u_)
p_n.assign(p_)
print('u max:', u_.vector().max())
print(t)
thank you for your reply!