Hello，I recently wanted to apply fenicsx to structural thermodynamic coupling calculations, where the elastic parameters of anisotropic materials varying with temperature, where μ=f1(T),E=f2(T),G=f3(T). I have tried its implementation below，but I don’t know much about the ufl language, and I don’t know how the global tempreture function(thetanew) is computed on the unit by C(thetanew),ufl.dot(C, strain2voigt(v)).I want to know if my implementation is correct, thanks!
import numpy as np
import ufl
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import dolfinx
from dolfinx import mesh, fem, plot, io
import petsc4py
L = 1
W = 0.2
T0 = 293.0
rho = 8.780e-9 # density, in kg/mm^3
def E_(Temp):
A1 = 138.81486
A2 = -8732.52764
A3 = 3387.36692
A4 = 476.71839
return (A2+(A1-A2)/(1+ufl.exp((Temp-A3-273.15)/A4)))*10**3
def G_(Temp):
A1 = 145.80245
A2 = -56756.35741
A3 = 4541.13659
A4 = 521.38506
return (A2+(A1-A2)/(1+ufl.exp((Temp-A3-273.15)/A4)))*10**3
def nu_(Temp):
a = 5.06136e-5
b = 0.34045
return (b+a*(Temp-273.15))
def alpha(Temp):
a = 0.00443
b = 10.62028
return (b+a*(Temp-273.15))*10**-6 # thermal diffusivity, in mm^2/s
def cp(Temp):
a = 0.34895
b = 355.7477
return (b+a*(Temp-273.15))*10**3 # heat capacity, in J/kg-K
def k(Temp):
a = 0.02109
b = 5.40548
return (b+a*(Temp-273.15))*10**-3 # thermal conductivity, in W/mm-K
def clamped_boundary(x):
return np.isclose(x[0], 0)
def boundary_D(x):
return np.isclose(x[0], 1)
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, W, W])],
[20,6,6], cell_type=mesh.CellType.hexahedron)
fdim = domain.topology.dim - 1
facets = mesh.locate_entities(domain, fdim, boundary_D)
facet_tag = mesh.meshtags(domain, fdim, facets, 1)
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
dx = ufl.dx
dt = dolfinx.fem.Constant(domain,c=0.)
ds = ufl.Measure("ds", domain=domain)
'''Thermal Problem'''
Vte = dolfinx.fem.FunctionSpace(domain, ("CG", 1))
theta_trial = ufl.TrialFunction(Vte)
theta_test = ufl.TestFunction(Vte)
thetaold = dolfinx.fem.Function(Vte)
thetaold.x.array[:] = T0
thetanew = dolfinx.fem.Function(Vte)
thetanew.name = "Tempreture"
h_robin = dolfinx.fem.Function(Vte)
t_wall = dolfinx.fem.Function(Vte)
therm_forml = ufl.inner(cp(thetaold)*rho*theta_trial / dt,theta_test)*dx + k(thetaold)*ufl.dot(ufl.grad(theta_trial), ufl.grad(theta_test))*dx
therm_formr = ufl.inner(cp(thetaold)*rho*(thetaold) / dt,theta_test)*dx - ufl.inner(h_robin*(t_wall-thetaold),theta_test)*ds(subdomain_data=facet_tag)
a_cpp_therm = dolfinx.fem.form(therm_forml)
f_cpp_therm = dolfinx.fem.form(therm_formr)
'''Mechanical problem'''
Vue = dolfinx.fem.VectorFunctionSpace(domain, ("CG", 2))
V_von_mises = dolfinx.fem.FunctionSpace(domain, ("CG", 1))
mecha_trial = ufl.TrialFunction(Vue)
mecha_test = ufl.TestFunction(Vue)
mechaold = dolfinx.fem.Function(Vue)
mechaold.x.array[:] = 0
mechanew = dolfinx.fem.Function(Vue)
mechanew.name = "displacement"
stresses = dolfinx.fem.Function(V_von_mises)
stresses.name = "VonmiseS"
def eps(v):
return ufl.sym(ufl.grad(v))
def strain2voigt(e):
"""
e is a 2nd-order tensor, returns its Voigt vectorial representation[ex,ey,ez,gxy,gyz,gzx]
"""
return ufl.as_vector([e[0,0],e[1,1],e[2,2],2*e[0,1],2*e[1,2],2*e[0,2]])
def voigt2stress(s):
"""
s is a stress-like vector (no 2 factor on last component) returns its tensorial representation;
"""
return ufl.as_tensor([[s[0], s[3], s[5]],
[s[3], s[1], s[4]],
[s[5], s[4], s[2]]])
def sigma(v,theta):
E,G,nu = E_(theta),G_(theta),nu_(theta)
C = np.array([[(1.-nu)*E/(1.-2.*nu)/(1.+nu), nu*E/(1.-2.*nu)/(1.+nu), nu*E/(1.-2.*nu)/(1.+nu), 0., 0., 0.],
[nu*E/(1.-2.*nu)/(1.+nu), (1.-nu)*E/(1.-2.*nu)/(1.+nu), nu*E/(1.-2.*nu)/(1.+nu), 0., 0., 0.],
[nu*E/(1.-2.*nu)/(1.+nu), nu*E/(1.-2.*nu)/(1.+nu), (1.-nu)*E/(1.-2.*nu)/(1.+nu), 0., 0., 0.],
[0., 0., 0., G, 0., 0.],
[0., 0., 0., 0., G, 0.],
[0., 0., 0., 0., 0., G]])
C = ufl.as_matrix(C)
return voigt2stress(ufl.dot(C, strain2voigt(eps(v))))
def thetasigma(v,theta):
E,G,nu = E_(theta),G_(theta),nu_(theta)
C = np.array([[(1.-nu)*E/(1.-2.*nu)/(1.+nu), nu*E/(1.-2.*nu)/(1.+nu), nu*E/(1.-2.*nu)/(1.+nu), 0., 0., 0.],
[nu*E/(1.-2.*nu)/(1.+nu), (1.-nu)*E/(1.-2.*nu)/(1.+nu), nu*E/(1.-2.*nu)/(1.+nu), 0., 0., 0.],
[nu*E/(1.-2.*nu)/(1.+nu), nu*E/(1.-2.*nu)/(1.+nu), (1.-nu)*E/(1.-2.*nu)/(1.+nu), 0., 0., 0.],
[0., 0., 0., G, 0., 0.],
[0., 0., 0., 0., G, 0.],
[0., 0., 0., 0., 0., G]])
C = ufl.as_matrix(C)
return voigt2stress(ufl.dot(C, strain2voigt(v)))
mecha_forml = ufl.inner(sigma(mecha_trial,thetanew), eps(mecha_test))*dx
mecha_formr = ufl.inner(thetasigma(alpha(thetanew)*(thetanew-T0)*ufl.Identity(len(mecha_trial)),thetanew), ufl.grad(mecha_test))*dx
a_cpp_mecha = dolfinx.fem.form(mecha_forml)
f_cpp_mecha = dolfinx.fem.form(mecha_formr)
'''Solve'''
Temp_results = dolfinx.io.VTKFile(MPI.COMM_WORLD, "./Result/Temp/Temp.pvd", "w")
Displacement_results = dolfinx.io.VTKFile(MPI.COMM_WORLD, "./Result/Displacement/Displacement.pvd", "w")
Sigma_results = dolfinx.io.VTKFile(MPI.COMM_WORLD, "./Result/Sigma/Sigma.pvd", "w")
dt.value = 1
dofs_D = fem.locate_dofs_geometrical(Vte, boundary_D)
h_robin.x.array[dofs_D] = -0.00056841
t_wall.x.array[dofs_D] = 1365
'''thermal solve'''
A_therm = dolfinx.fem.petsc.assemble_matrix(a_cpp_therm)
A_therm.assemble()
F_therm = dolfinx.fem.petsc.assemble_vector(f_cpp_therm)
F_therm.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)
ksp = petsc4py.PETSc.KSP()
ksp.create(domain.comm)
ksp.setOperators(A_therm)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
ksp.setFromOptions()
ksp.solve(F_therm, thetanew.vector)
thetanew.x.scatter_forward()
Temp_results.write_function(thetanew, 1)
Temp_results.close()
'''mechanic solve'''
A_mecha = dolfinx.fem.petsc.assemble_matrix(a_cpp_mecha)
A_mecha.assemble()
F_mecha = dolfinx.fem.petsc.assemble_vector(f_cpp_mecha)
F_mecha.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)
ksp = petsc4py.PETSc.KSP()
ksp.create(domain.comm)
ksp.setOperators(A_mecha)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
ksp.setFromOptions()
ksp.solve(F_mecha, mechanew.vector)
mechanew.x.scatter_forward()
Displacement_results.write_function(mechanew, 1)
Displacement_results.close()