H everyone,
I am currently trying to simulate a 3D heat equation problem, my heat source within my volume domain was initially described by:
Laser = Expression('2*A*Iden*exp(-pow((x[0] - 0), 2)/(w*w)-pow((x[1]-0), 2)/(w*w)-pow((x[2]-0), 2)/(w2*w2))',degree=3, A=0.25, Iden=Iden,w= 0.56, w2= 1)
When this is used, my code runs smoothly and iterates consistently with each time step. However, when I try and constrain this heat source by reducing w2 to 1/alpha and multiply the entire expression by alpha (alpha is the absorption coefficient of the material), my source term becomes this:
Laser = Expression('2*A*alpha*Iden*exp(-pow((x[0] - 0), 2)/(w*w)-pow((x[1]-0), 2)/(w*w)-pow((x[2]-0), 2)/(w2*w2))',degree=3, A=0.25, alpha=alpha, Iden=Iden,w= 0.56, w2= absorb_depth)
When I run the code here, initially compiles well with two or three iterations required each time step, but after a while, the number of iterations drops to zero:
Time 25.0
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 6.154e-10 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-03)
Newton iteration 1: r (abs) = 6.522e-13 (tol = 1.000e-10) r (rel) = 1.060e-03 (tol = 1.000e-03)
Newton solver finished in 1 iterations and 1 linear solver iterations.
Time 26.0
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 2.387e-10 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-03)
Newton iteration 1: r (abs) = 6.468e-13 (tol = 1.000e-10) r (rel) = 2.710e-03 (tol = 1.000e-03)
Newton solver finished in 1 iterations and 1 linear solver iterations.
Time 27.0
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 9.255e-11 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-03)
Newton solver finished in 0 iterations and 0 linear solver iterations.
Time 28.0
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 9.255e-11 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-03)
Newton solver finished in 0 iterations and 0 linear solver iterations.
When the solution is checked in Paraview, I find that no evolution of the system has taken place and has remained unchanged across all the timesteps. Why does computation of the solution after each timestep become very difficult when the heat source becomes constrained? Is there any way in which this system can be altered such the the problem can fully evolve with time?
My full code is shown below
from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
from ufl import as_tensor
from ufl import Index
import math
import mshr
parameters["allow_extrapolation"] = True
parameters["form_compiler"]["cpp_optimize"] = True
set_log_level(20)
# SOLVE THE HEAT TRANSFER PROBLEM FOR A 3D TARGET
#----------------
t_start = 0.0
t_end = 60
nstep = 60
dtime = (t_end - t_start)/ nstep #time step
#-----------------
#Configure for dimensions of system
#----------THERMAL--PROPERTIES--[Element: W]--------
kappa = 0.0172 #conductivity [W/mm K]
c = 0.132 #Specific Heat Capacity [J/gK]
rho = 0.01925 #Density [g/mm^3]
const = kappa /(c * rho)
tau_T = 0.0 #Temperature Gradient Lag
tau_q = 0.0 #Heat Flux Lag
alpha = 4.3943e+4 #absorption coefficient (mm^-1)
#---------------------------------------
pi = 3.141592653589793
T_am = 298 #ambient vacuum temperature (K)
T_a4 = T_am**4
epsilon = 0.25 # material emissivity
sigma = 5.67E-14 # W/(mm**2.K**4)
es = epsilon*sigma
Pmax = 500 #peak power of laser in W
w = 0.5 #mm
R = 1.5 #mm
area = pi*R*R #area of target (mm^2)
depth = 8.0 #thickness of target (mm)
volume = area * depth #volume of target
Iden = Pmax / (pi*w*w*volume*rho) #Intensity per unit mass
absorb_depth = 1/alpha
Laser = Expression('2*A*alpha*Iden*exp(-pow((x[0] - 0), 2)/(w*w)-pow((x[1]-0), 2)/(w*w)-pow((x[2]-0), 2)/(w2*w2))',degree=3, A=0.25, alpha=alpha, Iden=Iden,w= 0.56, w2= absorb_depth) #power (w2 localises the z-coordinates)
#A =! emissivity for system in thermal equillibrium
geometry = mshr.Cylinder(Point(0,0,0),Point(0,0,-8),R,R) #8mm long target
# Create mesh
mesh = generate_mesh(geometry, 50) # generate a mesh from the given geometry
Space = FunctionSpace((mesh), 'P', 1) #define finite element function space, defined via basis functions
VectorSpace = VectorFunctionSpace(mesh, 'P', 1)
cells = MeshFunction('size_t',mesh,mesh.topology().dim()) #codimension of 0
facets = MeshFunction('size_t',mesh,mesh.topology().dim()-1) #codimension of 1
da = Measure('ds', domain=mesh, subdomain_data = facets, metadata = {'quadrature_degree': 2}) #area element
dv = Measure('dx', domain=mesh, subdomain_data = cells, metadata = {'quadrature_degree': 2}) #volume element
initial_T = Expression("Tini", Tini=T_am, degree=3) # extrapolate an expression for the temperature before heating
T0 = interpolate(initial_T, Space)
T = Function(Space) # Temperature
V = TestFunction(Space) # Test Function used for FEA
dT = TrialFunction(Space) # Temperature Derivative
q0 = Function(VectorSpace) # heat flux at previous time step
i = Index()
G = as_tensor(T.dx(i), (i)) #gradient of T
G0 = as_tensor(T0.dx(i), (i)) # gradient of T at previous time step
q = as_tensor(dtime/(dtime + tau_q) * (tau_q/dtime*q0[i] - kappa*(1+tau_T/dtime)*G[i] + kappa*tau_T/dtime*G0[i]),(i)) #heat
F = (rho*c/dtime*(T-T0)*V - q[i]*V.dx(i) - rho*Laser*V ) * dv + es*(T**4 - T_a4)*V*da #final form to solve
Gain = derivative(F, T, dT) # Gain will be usedas the Jacobian required to determine the evolution of a dynamic system
file_T = File('target3D/solution.pvd')
for t in np.arange(t_start,t_end,dtime):
print( "Time", t)
solve(F==0, T, [], J = Gain, solver_parameters={"newton_solver":{"linear_solver": "mumps", "relative_tolerance": 1e-3} }, form_compiler_parameters={"cpp_optimize": True, "representation": "uflacs"} )
file_T << (T,t)
q_tmp = project(q, VectorSpace)
q0.assign(q_tmp)
T0.assign(T) #change so that T0 is equal to T for the next time step