I have issur with the thermal boundary condition i need to define it T depend of t (from 0 to 60s) T(t)= 298K and T(60)= 350K, in this code the temperature stay stable at 298 K, i idea is to heat the lower part of the domain to reach the 350K to see the final temperature distribution on the domain
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
# Function to create a mesh and function space
def create_simulation_environment():
mesh = RectangleMesh(Point(0, 0), Point(0.06, 0.001), 120, 20) # Refined mesh
V_u = VectorFunctionSpace(mesh, 'P', 1) # Displacement function space
V = FunctionSpace(mesh, 'DG', 0)
V_T = FunctionSpace(mesh, 'P', 1) # Temperature function space
return mesh, V_u, V, V_T
# Function to initialize property functions
def initialize_properties(V):
properties = {
'young_modulus': Function(V),
'thermal_expansion': Function(V),
'conductivity_x': Function(V),
'conductivity_y': Function(V),
'density': Function(V),
'specific_heat': Function(V),
'temperature': Function(V),
'Poisson': Function(V),
'beta': Function(V)
}
return properties
# Function to compute Young's Modulus based on temperature
def compute_young_modulus(temp, material):
if material == "Tango":
if temp <= 223.668:
return 1664000000
elif temp <= 232.181:
return 1386400000
elif temp <= 240.339:
return 1035300000
elif temp <= 248.14:
return 668100000
elif temp <= 254.167:
return 359200000
elif temp <= 259.128:
return 173060000
elif temp <= 263.024:
return 74740000
elif temp <= 265.855:
return 32280000
elif temp <= 268.687:
return 13441000
elif temp <= 271.163:
return 5597000
elif temp <= 274.35:
return 2417000
elif temp <= 278.246:
return 1043800
elif temp <= 283.562:
return 503000
elif temp <= 290.298:
return 290900
elif temp <= 298.455:
return 209430
elif temp <= 306.613:
return 174490
elif temp <= 315.482:
return 156390
elif temp <= 323.996:
return 156390
elif temp <= 332.865:
return 156390
elif temp <= 341.734:
return 156390
elif temp <= 350.604:
return 162210
elif temp <= 359.119:
return 168240
elif temp <= 367.278:
return 168240
else:
return 168240 # Default value for temperatures above the specified range
elif material == "Vero":
if temp <= 223.096:
return 2964000000
elif temp <= 228.196:
return 2845000000
elif temp <= 231.433:
return 2845000000
elif temp <= 237.153:
return 2742600000
elif temp <= 261.166:
return 2435800000
elif temp <= 265.844:
return 2269500000
elif temp <= 270.433:
return 2173200000
elif temp <= 275.976:
return 2043400000
elif temp <= 283.535:
return 1702400000
elif temp <= 289.721:
return 1522400000
elif temp <= 295.329:
return 1217300000
elif temp <= 301.845:
return 882500000
elif temp <= 307.939:
return 544000000
elif temp <= 314.297:
return 273140000
elif temp <= 317.531:
return 172280000
elif temp <= 321.03:
return 93470000
elif temp <= 322.469:
return 71080000
elif temp <= 324.329:
return 49120000
elif temp <= 326.675:
return 30560000
elif temp <= 328.491:
return 21261000
elif temp <= 332.011:
return 11404000
elif temp <= 334.514:
return 7917000
elif temp <= 338.391:
return 5470000
elif temp <= 345.95:
return 3947000
elif temp <= 351.072:
return 3815000
elif temp <= 358.832:
return 3780000
elif temp <= 368.278:
return 3763000
elif temp <= 373.156:
return 3815000
else:
return 3815000 # Default value for temperatures above the specified range
# Function to compute Thermal Expansion based on temperature
def compute_thermal_expansion(temp, material):
if material == "Tango":
if temp < 273:
return 0
elif temp <= 300.15:
return -0.9e-4
elif temp <= 310:
return -1.2e-4
elif temp <= 320:
return -1.2e-4
elif temp <= 343.15:
return -1.3e-4
else:
return -1.3e-4 # Default value for temperatures above the specified range
elif material == "Vero":
if temp < 273:
return 0
elif temp <= 300.15:
return -1.3e-4
elif temp <= 310:
return -1.35e-4
elif temp <= 320:
return -1.4e-4
elif temp <= 343.15:
return -1.6e-4
else:
return -1.6e-4 # Default value for temperatures above the specified range
# Function to compute thermal conductivity based on temperature
def compute_conductivity(temp, material, direction):
if material == "Tango":
return 1.4 # Replace with actual function of T if needed
elif material == "Vero":
return 0.14 # Replace with actual function of T if needed
# Function to compute effective thermal expansion coefficient
def compute_effective_thermal_expansion(temp, material):
alpha = compute_thermal_expansion(temp, material)
E = compute_young_modulus(temp, material)
nu = 0.4 # Poisson's ratio
beta = E * alpha / (1 - 2 * nu)
return beta
# Function to assign properties based on temperature
def assign_properties(properties, temp, mesh):
properties['temperature'].vector()[:] = temp
for cell in cells(mesh):
y_mid = cell.midpoint().y()
if y_mid < 0.0005: # TangoBlackPlus
nu = 0.4
alpha = compute_thermal_expansion(temp, "Tango")
E = compute_young_modulus(temp, "Tango")
beta = compute_effective_thermal_expansion(temp, "Tango")
properties['Poisson'].vector()[cell.index()] = nu
properties['conductivity_x'].vector()[cell.index()] = compute_conductivity(temp, "Tango", 'x')
properties['conductivity_y'].vector()[cell.index()] = compute_conductivity(temp, "Tango", 'y')
properties['density'].vector()[cell.index()] = 1200
properties['specific_heat'].vector()[cell.index()] = 1300
properties['young_modulus'].vector()[cell.index()] = E
properties['thermal_expansion'].vector()[cell.index()] = alpha
properties['beta'].vector()[cell.index()] = beta
else: # VeroWhitePlus
nu = 0.4
alpha = compute_thermal_expansion(temp, "Vero")
E = compute_young_modulus(temp, "Vero")
beta = compute_effective_thermal_expansion(temp, "Vero")
properties['Poisson'].vector()[cell.index()] = nu
properties['conductivity_x'].vector()[cell.index()] = compute_conductivity(temp, "Vero", 'x')
properties['conductivity_y'].vector()[cell.index()] = compute_conductivity(temp, "Vero", 'y')
properties['density'].vector()[cell.index()] = 1050
properties['specific_heat'].vector()[cell.index()] = 2300
properties['young_modulus'].vector()[cell.index()] = E
properties['thermal_expansion'].vector()[cell.index()] = alpha
properties['beta'].vector()[cell.index()] = beta
# Function to define heat generation (assumed to be 10 for both materials)
def heat_generation(material):
return 10 # Modify this if there is internal heat generation
# Define the Linear Plane Strain Elasticity Problem
def solve_elasticity(properties, mesh, V_u):
# Define trial and test functions for displacement
u = TrialFunction(V_u)
v = TestFunction(V_u)
# Extract material properties
E = properties['young_modulus']
nu = properties['Poisson']
mu = E / (2 * (1 + nu))
lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))
rho = properties['density']
# Define strain and stress
def epsilon(u):
return sym(grad(u))
def sigma(u):
return lambda_ * tr(epsilon(u)) * Identity(len(u)) + 2 * mu * epsilon(u)
# Define gravitational body force for both Tango and Vero
g = 9.81
rho_tango = 1200 # Density for Tango material in kg/m^3
rho_vero = 1050 # Density for Vero material in kg/m^3
# Define the subdomain markers
subdomains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
for cell in cells(mesh):
if cell.midpoint().y() < 0.0005:
subdomains[cell] = 1 # Tango material
else:
subdomains[cell] = 2 # Vero material
dx = Measure('dx', domain=mesh, subdomain_data=subdomains)
# Define gravitational body force for each material
b_tango = as_vector((0, -rho_tango * g))
b_vero = as_vector((0, -rho_vero * g))
# Define variational problem for displacement with gravity
a = inner(sigma(u), epsilon(v)) * dx
L = dot(b_tango, v) * dx(1) + dot(b_vero, v) * dx(2) # Gravity term for both materials
# Define boundary conditions for displacement
zero_displacement = Constant((0.0, 0.0))
bc_left = DirichletBC(V_u, zero_displacement, 'near(x[0], 0)')
bc_right = DirichletBC(V_u, zero_displacement, 'near(x[0], 0.06)')
bcs_u = [bc_left, bc_right]
# Compute solution for displacement
u_sol = Function(V_u)
solve(a == L, u_sol, bcs_u)
# Calculate volume strain
volume_strain = div(u_sol)
# Project volume strain onto a function space for visualization or further use
V_s = FunctionSpace(mesh, 'P', 1)
volume_strain_projected = project(volume_strain, V_s)
return u_sol, volume_strain_projected
# Function to calculate thermal stress
def calculate_thermal_stress(properties, T_grad, u_sol, mesh):
beta = properties['beta']
epsilon_u = sym(grad(u_sol))
a = as_vector([1.0, 1.0]) # Matching dimensions with T_grad
D = as_matrix([[properties['young_modulus'] / (1 - properties['Poisson'] ** 2), properties['Poisson'] * properties['young_modulus'] / (1 - properties['Poisson'] ** 2)],
[properties['Poisson'] * properties['young_modulus'] / (1 - properties['Poisson'] ** 2), properties['young_modulus'] / (1 - properties['Poisson'] ** 2)]])
# Calculate the thermal stress
sigma = -beta * dot(T_grad, a) * Identity(2) + D * epsilon_u
return sigma
# Function to compute von Mises stress
def compute_von_mises_stress(sigma):
s = sigma - (1.0 / 3) * tr(sigma) * Identity(2) # deviatoric stress
von_mises = sqrt(3.0 / 2 * inner(s, s))
return von_mises
# Define the time-dependent temperature function
def temperature_function(t):
return 298 + (350 - 298) / 60 * t # Linear interpolation
# Time values
t_values = np.linspace(0, 60, 13) # Adjusted to go up to 60 for heat problem
def main():
mesh, V_u, V, V_T = create_simulation_environment()
properties = initialize_properties(V)
# Define the initial condition for temperature
T = Function(V_T)
T.assign(interpolate(Constant(298.0), V_T)) # Initial temperature set to 298 K
T_previous = T.copy(deepcopy=True) # Store initial temperature
# Define the boundary condition only on the lower surface from (0, 0) to (0.06, 0)
class LowerSurfaceBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0) and (0 <= x[0] <= 0.06) and on_boundary
lower_surface_boundary = LowerSurfaceBoundary()
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
lower_surface_boundary.mark(boundary_markers, 1)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
# Initialize T_n with the initial temperature
T_n = Function(V_T)
T_n.assign(T)
# Time-stepping
t = 0
dt = t_values[1] - t_values[0]
for n in range(len(t_values)):
t = t_values[n]
temp = temperature_function(t)
# Assign properties based on the current temperature
assign_properties(properties, temp, mesh)
# Solve elasticity problem to get volume strain
u_sol, volume_strain_projected = solve_elasticity(properties, mesh, V_u)
# Define the problem for each subdomain in the heat equation
v = TestFunction(V_T)
F_tango = (
+properties['density'] * properties['specific_heat'] * (T - T_n) / dt * v * dx(1)
- properties['conductivity_x'] * dot(grad(T)[0], grad(v)[0]) * dx(1)
- properties['conductivity_y'] * dot(grad(T)[1], grad(v)[1]) * dx(1)
- heat_generation("Tango") * v * dx(1)
+ properties['beta'] * T * volume_strain_projected * v * dx(1)
)
F_vero = (
properties['density'] * properties['specific_heat'] * (T - T_n) / dt * v * dx(2)
- properties['conductivity_x'] * dot(grad(T)[0], grad(v)[0]) * dx(2)
- properties['conductivity_y'] * dot(grad(T)[1], grad(v)[1]) * dx(2)
- heat_generation("Vero") * v * dx(2)
+ properties['beta'] * T * volume_strain_projected * v * dx(2)
)
# Add heat flux boundary condition to the variational formulation
F = F_tango + F_vero + 0 * v * ds(1) # No flux but keeps the form
# Define the Jacobian
J = derivative(F, T)
# Adjust solver parameters
problem = NonlinearVariationalProblem(F, T, bcs=None, J=J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-8
prm['newton_solver']['relative_tolerance'] = 1E-9
prm['newton_solver']['maximum_iterations'] = 100 # Increased iterations
prm['newton_solver']['relaxation_parameter'] = 0.7 # Adjusted relaxation parameter
# Add more diagnostic output
prm['newton_solver']['report'] = True
prm['newton_solver']['error_on_nonconvergence'] = False
# Solve the problem
try:
solver.solve()
except RuntimeError as e:
print(f"Solver failed to converge at time step {n}, t = {t:.2f}. Error: {e}")
break
T_previous.assign(T)
T_n.assign(T) # Update T_n for the next time step
# Debugging output
print(f"Time step {n}, t = {t:.2f}, max temperature: {T.vector().max()}")
print(f"Temperature distribution at time step {n}: {T.vector().get_local()}") # Added debug info
# Update plot if desired (e.g., every few steps)
if n % 10 == 0:
plt.figure()
p = plot(T)
plt.title(f'Temperature at t={t:.1f}')
plt.colorbar(p)
plt.show()
# Final plot
plt.figure()
p = plot(T)
plt.title('Final Temperature Distribution')
plt.colorbar(p)
plt.show()
# Project the temperature gradient for the thermal stress calculation
T_grad = project(grad(T), VectorFunctionSpace(mesh, 'P', 1))
# Calculate thermal stress
sigma = calculate_thermal_stress(properties, T_grad, u_sol, mesh)
# Compute von Mises stress
von_mises_stress = compute_von_mises_stress(sigma)
# Project von Mises stress to a scalar function for plotting
V_sig = FunctionSpace(mesh, 'P', 1)
von_mises_stress_projected = project(von_mises_stress, V_sig)
# Print von Mises stress value
print(f'Von Mises stress: {von_mises_stress_projected.vector().max()} Pa')
# Plot von Mises stress
plt.figure()
p = plot(von_mises_stress_projected)
plt.title('Von Mises Stress')
plt.colorbar(p)
plt.show()
# Plot displacement components
plt.figure()
p = plot(u_sol.sub(0), title='Displacement U1 (X-axis)')
plt.colorbar(p)
plt.show()
plt.figure()
p = plot(u_sol.sub(1), title='Displacement U2 (Y-axis)')
plt.colorbar(p)
plt.show()
# Plot deformed mesh using standard FEniCS method
plt.figure()
plot(mesh, title="Original mesh")
plt.show()
plt.figure()
plot(mesh, title="Deformed mesh")
plot(u_sol, mode="displacement")
plt.show()
# Project and plot the strain
epsilon_u = sym(grad(u_sol))
V_eps = TensorFunctionSpace(mesh, 'P', 1)
epsilon_projected = project(epsilon_u, V_eps)
# Plot strain components
plt.figure()
p = plot(epsilon_projected[0, 0])
plt.title('Strain Component epsilon_xx')
plt.colorbar(p)
plt.show()
plt.figure()
p = plot(epsilon_projected[1, 1])
plt.title('Strain Component epsilon_yy')
plt.colorbar(p)
plt.show()
plt.figure()
p = plot(epsilon_projected[0, 1])
plt.title('Strain Component epsilon_xy')
plt.colorbar(p)
plt.show()
if __name__ == "__main__":
main()