import mshr
from fenics import *
outer_radius = 0.6
inner_radius = 0.5
length = 10.0
resolution = 50
dt = Constant(0.01) # or any other value that is appropriate for your simulation
outer_cylinder = mshr.Cylinder(Point(0, 0, 0), Point(0, 0, length), outer_radius, outer_radius)
inner_cylinder = mshr.Cylinder(Point(0, 0, 0), Point(0, 0, length), inner_radius, inner_radius)
hollow_cylinder = outer_cylinder - inner_cylinder
mesh = mshr.generate_mesh(hollow_cylinder, resolution)
# Fluid Flow Analysis
# -------------------
# Define elements for velocity and pressure
V_element = VectorElement('P', mesh.ufl_cell(), 2)
Q_element = FiniteElement('P', mesh.ufl_cell(), 1)
# Create mixed space
W_flow = FunctionSpace(mesh, MixedElement([V_element, Q_element]))
# Split the mixed functions to access the velocity and pressure parts
(u_flow, p_flow) = TrialFunctions(W_flow)
(v_flow, q_flow) = TestFunctions(W_flow)
p_n = Function(W_flow.sub(1)) # Define in the pressure part of the mixed space
# Define boundary condition for fluid flow
def fluid_inlet(x, on_boundary):
return near(x[2], 0) and on_boundary
def fluid_outlet(x, on_boundary):
return near(x[2], length) and on_boundary
inflow_profile = Expression(('4*1.5*(1 - ((x[0]*x[0]) + (x[1]*x[1])) / pow(outer_radius, 2))', '0', '0'),
degree=2, outer_radius=outer_radius)
bc_flow_in = DirichletBC(W_flow.sub(0), inflow_profile, fluid_inlet)
bc_flow_out = DirichletBC(W_flow.sub(1), Constant(0), fluid_outlet)
bc_flow = [bc_flow_in, bc_flow_out]
# Define function space for velocity
V = VectorFunctionSpace(mesh, 'P', 2)
# Define function for solution at previous and current time steps
u_n = Function(V)
u_flow = Function(V)
# Define expressions used in variational forms
U = 0.5 * (u_n + u_flow)
n = FacetNormal(mesh)
f = Constant((0, 0, 0))
k = Constant(dt)
mu_fluid = 0.001 # Dynamic viscosity for the fluid in Pa.s (or any appropriate value)
mu = Constant(mu_fluid)
rho_fluid=1000
rho = Constant(rho_fluid)
# 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_flow - u_n) / k, v_flow) * dx \
+ rho * dot(dot(u_flow, nabla_grad(u_flow)), v_flow) * dx \
+ inner(sigma(U, p_flow), epsilon(v_flow)) * dx \
+ dot(p_flow * n, v_flow) * ds - dot(mu * nabla_grad(U) * n, v_flow) * ds \
- dot(f, v_flow) * dx
a1 = lhs(F1)
L1 = rhs(F1)
# Define variational problem for step 2
a2 = dot(nabla_grad(p_flow), nabla_grad(q_flow)) * dx
L2 = dot(nabla_grad(p_flow - p_n), nabla_grad(q_flow)) * dx - (1 / k) * div(u_) * q_flow * dx
# Define variational problem for step 3
a3 = dot(u_flow, v_flow) * dx
L3 = dot(u_, v_flow) * dx - k * dot(nabla_grad(p_ - p_n), v_flow) * dx
# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
# Thermal Analysis
# ---------------
# Define function space
V_thermal = FunctionSpace(mesh, 'P', 1)
# Define boundary condition for thermal problem
def boundary(x, on_boundary):
return on_boundary
bc_thermal = DirichletBC(V_thermal, Constant(300), 'on_boundary')
# Define variational problem
u_thermal = TrialFunction(V_thermal)
v_thermal = TestFunction(V_thermal)
f_thermal = Constant(10e3 / (outer_radius * 2 * pi * length)) # Heat flux
a_thermal = dot(grad(u_thermal), grad(v_thermal)) * dx
L_thermal = f_thermal * v_thermal * ds
# Compute solution for thermal problem
u_thermal = Function(V_thermal)
solve(a_thermal == L_thermal, u_thermal, bc_thermal)
# Mechanical Analysis
# -------------------
# Define function space
V_mechanical = VectorFunctionSpace(mesh, 'P', 1)
# Define boundary condition for mechanical problem
def left_boundary(x, on_boundary):
return near(x[2], 0)
bc_mechanical = DirichletBC(V_mechanical, Constant((0, 0, 0)), 'near(x[2], 0)')
# Define variational problem
u_mechanical = TrialFunction(V_mechanical)
v_mechanical = TestFunction(V_mechanical)
T = Constant(0) # Traction at the right end
thermal_expansion_coefficient = Constant(1e-5) # Material property
E = Constant(200e9) # Young's modulus
nu = Constant(0.3) # Poisson's ratio
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / 2 / (1 + nu)
# Thermal strain tensor
thermal_expansion_coefficient = Constant(1e-5) # Material property
reference_temperature = Constant(300) # Reference temperature (adjust as needed)
epsilon_th = thermal_expansion_coefficient * (u_thermal - reference_temperature) * Identity(3)
# Total strain tensor
epsilon = sym(grad(u_mechanical)) + epsilon_th
# Stress
sigma = lmbda * tr(epsilon) * Identity(3) + 2 * mu * epsilon
class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[2], length) and on_boundary
right_boundary = RightBoundary()
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
right_boundary.mark(boundaries, 1)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
# Define traction vector (adjust as needed)
T = Constant((0, 0, 0)) # Traction at the right end
# Weak form
n = FacetNormal(mesh) # Normal to the boundary
F_mechanical = inner(sigma, grad(v_mechanical)) * dx - dot(T, v_mechanical) * ds(1)
a_mechanical, L_mechanical = lhs(F_mechanical), rhs(F_mechanical)
# Compute solution for mechanical problem
u_mechanical = Function(V_mechanical)
solve(a_mechanical == L_mechanical, u_mechanical, bc_mechanical)
# Output results
File("temperature.pvd") << u_thermal
File("displacement.pvd") << u_mechanical
I am trying to simulate a simple case of pipe flow with hot wall to see how the walls would be deform of the pipe. Please help me debug this code. I have trying for days now.