# Inreasing of temperature during time

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):

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
beta = properties['beta']

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)
- 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)
- 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)

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

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

# 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
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()

``````

Please try to reduce the code to a minimal reproducible example. The code you have produced introduces so many other concepts, that is very time consuming for anyone to go through in detail.
It seems like you would like to have a time dependent `T`, and in your code you seem to have a temporal loop in your code, which would make it time dependent. It is to me unclear what you mean with

Do you want to prescribe this data to your problem?

If not, I would consider starting with the simple heat equation tutorial to teach you time dependency, see chapter 3.1 of:
https://fenicsproject.org/pub/tutorial/pdf/fenics-tutorial-vol1.pdf