# Trying to Model Thermo-Mechanical Case in a Simple Pipe

``````import mshr
from fenics import *

length = 10.0
resolution = 50
dt = Constant(0.01) # or any other value that is appropriate for your simulation

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, 0) and on_boundary

def fluid_outlet(x, on_boundary):
return near(x, length) and on_boundary

inflow_profile = Expression(('4*1.5*(1 - ((x*x) + (x*x)) / pow(outer_radius, 2))', '0', '0'),
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)

def epsilon(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
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
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, 0)

bc_mechanical = DirichletBC(V_mechanical, Constant((0, 0, 0)), 'near(x, 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
# Stress
sigma = lmbda * tr(epsilon) * Identity(3) + 2 * mu * epsilon

class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x, 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.