Elastic parameters of anisotropic materials varying with temperature

dolfinx
1

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!

image
image782×313 6.23 KB

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

I have a need for something similar - did this implementation work out for you? It makes sense to me. (It seems odd to have the thermal problem as a function of dt without a time-stepping loop at the end, but you may have left that out for brevity.)