This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License
这项工作根据知识共享署名-相同方式共享 4.0 国际许可协议授权

The corresponding files can be obtained from:
相应的文件可以从以下地址获取：

• Jupyter Notebook: thermoelasticity_transient.ipynb
• Python script: thermoelasticity_transient.py  Python 脚本: thermoelasticity_transient.py

Thermo-elastic evolution problem (full coupling)
热弹性演化问题（完全耦合） ¶

Introduction  引言 ¶

In this tour, we will solve a transient thermoelastic evolution problem in which both thermo-mechanical fields are fully coupled, we will however assume that the evolution is quasi-static and will hence neglect inertial effects. Note that a staggered approach could also have been adopted in which one field is calculated first (say the temperature for instance) using predicted values of the other field and similarly for the other field in a second step (see for instance [FAR91]).
在这个教程中，我们将解决一个瞬态热弹性演化问题，其中热力学场完全耦合，但我们将假设演化是准静态的，因此将忽略惯性效应。请注意，也可以采用交错方法，其中先计算一个场（例如温度）的值，使用预测的另一个场的值，然后在第二步中类似地计算另一个场（例如参考文献[FAR91]）。

Static elastic computation with thermal strains is treated in the LinearThermoelasticity tour.
带热应变的静态弹性计算在 LinearThermoelasticity 教程中处理。

Problem position  问题位置 ¶

The problem consists of a quarter of a square plate perforated by a circular hole. Thermo-elastic properties are isotropic and correspond to those of aluminium (note that stress units are in MPa$\text{MPa}$, distances in m$\text{m}$ and mass in kg$\text{kg}$). Linearized thermo-elasticity will be considered around a reference temperature of T0=293 K$T_0=293\text{ K}$. A temperature increase of ΔT=+10∘C$\mathrm{\Delta }T=+{10}^{\circ }\text{C}$ will be applied on the hole boundary. Symmetry conditions are applied on the corresponding symmetry planes and stress and flux-free boundary conditions are adopted on the plate outer boundary.
该问题包括一个被圆形孔洞穿透的正方形板的四分之一。热弹性性质是各向同性的，对应于铝的性质（注意应力单位是 MPa$\text{MPa}$ ，距离单位是 m$\text{m}$ ，质量单位是 kg$\text{kg}$ ）。将在参考温度 T0=293 K$T_0=293\text{ K}$ 附近考虑线性化热弹性。将在孔洞边界上施加温度升高 ΔT=+10∘C$\mathrm{\Delta }T=+{10}^{\circ }\text{C}$ 。在相应的对称平面上应用对称条件，并在板的外边界上采用无应力和无通量的边界条件。

We first import the relevant modules and define the mesh and material parameters (see the next section for more details on the parameters).
我们首先导入相关模块，并定义网格和材料参数（有关参数的更多详细信息，请参见下一节）。

from dolfin import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
%matplotlib notebook

L = 1.
R = 0.1
N = 50  # mesh density

domain = Rectangle(Point(0., 0.), Point(L, L)) - Circle(Point(0., 0.), R, 100)
mesh = generate_mesh(domain, N)

T0 = Constant(293.)
DThole = Constant(10.)
E = 70e3
nu = 0.3
lmbda = Constant(E*nu/((1+nu)*(1-2*nu)))
mu = Constant(E/2/(1+nu))
rho = Constant(2700.)     # density
alpha = 2.31e-5  # thermal expansion coefficient
kappa = Constant(alpha*(2*mu + 3*lmbda))
cV = Constant(910e-6)*rho  # specific heat per unit volume at constant strain
k = Constant(237e-6)  # thermal conductivity

We now define the relevant FunctionSpace for the considered problem. Since we will adopt a monolithic approach i.e. in which both fields are coupled and solved at the same time, we will need to resort to a Mixed FunctionSpace for both the displacement u$\mathit{u}$ and the temperature variation Θ=T−T0$\mathrm{\Theta }=T-T_0$. For an introduction on the use of Mixed FunctionSpace, check out the FEniCS tutorials on the mixed Poisson equation or the Stokes problem. Let us just point out that the constructor using MixedFunctionSpace has been deprecated since version 2016.2 (see this post). In the following code, MixedElement([Vue, Vte]) could also be replaced by Vue * Vte.
我们现在定义所考虑问题的相关 FunctionSpace。由于我们将采用整体方法，即同时耦合并求解两个场，因此我们需要为位移 u$\mathit{u}$ 和温度变化 Θ=T−T0$\mathrm{\Theta }=T-T_0$ 采用混合 FunctionSpace。有关混合 FunctionSpace 使用的介绍，请查看 FEniCS 教程中的混合泊松方程或 Stokes 问题。我们只需指出，使用 MixedFunctionSpace 的构造函数自 2016.2 版本起已被弃用（请参阅此帖子）。在以下代码中， MixedElement([Vue, Vte]) 也可以被替换为 Vue * Vte 。

Vue = VectorElement('CG', mesh.ufl_cell(), 2) # displacement finite element
Vte = FiniteElement('CG', mesh.ufl_cell(), 1) # temperature finite element
V = FunctionSpace(mesh, MixedElement([Vue, Vte]))

def inner_boundary(x, on_boundary):
    return near(x[0]**2+x[1]**2, R**2, 1e-3) and on_boundary
def bottom(x, on_boundary):
    return near(x[1], 0) and on_boundary
def left(x, on_boundary):
    return near(x[0], 0) and on_boundary

bc1 = DirichletBC(V.sub(0).sub(1), Constant(0.), bottom)
bc2 = DirichletBC(V.sub(0).sub(0), Constant(0.), left)
bc3 = DirichletBC(V.sub(1), DThole, inner_boundary)
bcs = [bc1, bc2, bc3]

Dirichlet boundary conditions must be defined from the full FunctionSpace V using the appropriate subspaces that is V.sub(0) for the displacement (and .sub(0) or .sub(1) for the corresponding x/y component) and V.sub(1) for the temperature. Note also that in the following, we will in fact work with the temperature variation Θ=T−T0$\mathrm{\Theta }=T-T_0$ as a field unkwown instead of the total temperature. Hence, the boundary condition on the hole boundary reads indeed as ΔT=+10∘C$\mathrm{\Delta }T=+{10}^{\circ }\text{C}$ .
Dirichlet 边界条件必须使用完整的 FunctionSpace V 并通过适当的子空间定义，即 V.sub(0) 用于位移（以及 .sub(0) 或 .sub(1) 用于相应的 x/y 分量）和 V.sub(1) 用于温度。请注意，在以下内容中，我们将实际上使用温度变化 Θ=T−T0$\mathrm{\Theta }=T-T_0$ 作为未知场，而不是总温度。因此，孔边界上的边界条件确实读作 ΔT=+10∘C$\mathrm{\Delta }T=+{10}^{\circ }\text{C}$ 。

Note: The definition of the inner_boundary region used the syntax near(expr, value, tolerance) to define the region where expr = value up to the specified tolerance. The given tolerance is quite large given that mshr generates a faceted approximation of a Circle (here using 100 segments). Mesh nodes may therefore not lie exactly on the true circle.
注意： inner_boundary 区域的定义使用了 near(expr, value, tolerance) 语法来定义从 expr = value 到指定 tolerance 的区域。考虑到 mshr 生成了一个由多边形近似表示的圆（这里使用了 100 个分段），给定的容差相当大。因此，网格节点可能不会完全位于真正的圆上。

Variational formulation and time discretization
变分格式和时间离散化 ¶

The linearized thermoelastic constitutive equations are given by:
线性化热弹性本构方程为：

σ=C:(ε−α(T−T0)1)=λtr(ε)1+2με−κ(T−T0)1$$\mathit{\sigma }=\mathbb{C}:(\mathit{\epsilon }-\alpha (T-T_0)\mathbf{1})=\lambda \text{tr}(\mathit{\epsilon })\mathbf{1}+2\mu \mathit{\epsilon }-\kappa (T-T_0)\mathbf{1}$$

ρs=ρs0+ρCεT0(T−T0)+κtr(ε)$$\rho s=\rho s_0+{\displaystyle \frac{\rho C_{\epsilon }}{T_0}}(T-T_0)+\kappa \text{tr}(\mathit{\epsilon })$$

• λ,μ$\lambda ,\mu$ the Lamé coefficients
  λ,μ$\lambda ,\mu$ Lamé系数
• ρ$\rho$ material density   ρ$\rho$ 材料密度
• α$\alpha$ thermal expansion coefficient
  α$\alpha$ 热膨胀系数
• κ=α(3λ+2μ)$\kappa =\alpha (3\lambda +2\mu )$
• Cε$C_{\epsilon }$ the specific heat at constant strain (per unit of mass).
  Cε$C_{\epsilon }$ 恒定应变下的比热容（单位质量）。
• s$s$ (resp. s0$s_0$) the entropy per unit of mass in the current (resp. initial configuration)
  s$s$ （或 s0$s_0$ ）当前（或初始构型）的单位质量熵

These equations are completed by the equilibrium equation which will later be expressed in its weak form (virtual work principle) and the linearized heat equation (without source terms):
这些方程由平衡方程补充，该方程稍后将以其弱形式（虚功原理）和线性化热方程（无源项）表示：

ρT0s˙+divq=0$$\rho T_0\dot{s}+\text{div}\mathit{q}=0$$

where the heat flux is related to the temperature gradient through the isotropic Fourier law: q=−k∇T$\mathit{q}=-k\mathrm{\nabla }T$ with k$k$ being the thermal conductivity. Using the entropy constitutive relation, the weak form of the heat equation reads as:
其中热通量通过各向同性傅里叶定律与温度梯度相关： q=−k∇T$\mathit{q}=-k\mathrm{\nabla }T$ ，其中 k$k$ 为热导率。使用熵本构关系，热方程的弱形式为：

∫ΩρT0s˙TˆdΩ−∫Ωq⋅∇TˆdΩ=−∫∂Ωq⋅nTˆdS∀Tˆ∈VT$${\int }_{\mathrm{\Omega }}\rho T_0\dot{s}\hat{T}d\mathrm{\Omega }-{\int }_{\mathrm{\Omega }}\mathit{q}\cdot \mathrm{\nabla }\hat{T}d\mathrm{\Omega }=-{\int }_{\mathrm{\partial }\mathrm{\Omega }}\mathit{q}\cdot \mathit{n}\hat{T}dS\mathrm{\forall }\hat{T}\in V_T$$

∫Ω(ρCεT˙+κT0tr(ε˙))TˆdΩ+∫Ωk∇T⋅∇TˆdΩ=∫∂Ωk∂nTTˆdS∀Tˆ∈VT$${\int }_{\mathrm{\Omega }}(\rho C_{\epsilon }\dot{T}+\kappa T_0\text{tr}(\dot{\mathit{\epsilon }}))\hat{T}d\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}k\mathrm{\nabla }T\cdot \mathrm{\nabla }\hat{T}d\mathrm{\Omega }={\int }_{\mathrm{\partial }\mathrm{\Omega }}k{\mathrm{\partial }}_{n}T\hat{T}dS\mathrm{\forall }\hat{T}\in V_T$$

with VT$V_T$ being the FunctionSpace for the temperature field.
VT$V_T$ 是温度场的 FunctionSpace。

The time derivatives are now replaced by an implicit Euler scheme, so that the previous weak form at the time increment n+1$n+1$ is now:
现在时间导数被隐式欧拉格式替代，因此时间增量 n+1$n+1$ 之前的弱形式现在是：

∫Ω(ρCεT−TnΔt+κT0tr(ε−εnΔt))TˆdΩ+∫Ωk∇T⋅∇TˆdΩ=∫∂Ωk∂nTTˆdS∀Tˆ∈VT(1)$${\int }_{\mathrm{\Omega }}(\rho C_{\epsilon }{\displaystyle \frac{T-T_n}{\mathrm{\Delta }t}}+\kappa T_0\text{tr}\left({\displaystyle \frac{\mathit{\epsilon }-{\mathit{\epsilon }}_n}{\mathrm{\Delta }t}}\right))\hat{T}d\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}k\mathrm{\nabla }T\cdot \mathrm{\nabla }\hat{T}d\mathrm{\Omega }={\int }_{\mathrm{\partial }\mathrm{\Omega }}k{\mathrm{\partial }}_{n}T\hat{T}dS\mathrm{\forall }\hat{T}\in V_T(1)$$

where T$T$ and ε$\mathit{\epsilon }$ correspond to the unknown fields at the time increment n+1$n+1$. For more details on the time discretization of the heat equation, see also the Heat equation FEniCS tutorial.
其中 T$T$ 和 ε$\mathit{\epsilon }$ 对应于时间增量 n+1$n+1$ 时的未知场。有关热方程时间离散的更多细节，请参阅热方程 FEniCS 教程。

In addition to the previous thermal weak form, the mechanical weak form reads as:
除了之前的热弱形式，力学弱形式读作：

∫Ω(λtr(ε)1+2με−κ(T−T0)1):∇svˆ dΩ=Wext(vˆ)∀vˆ∈VU(2)$${\int }_{\mathrm{\Omega }}(\lambda \text{tr}(\mathit{\epsilon })\mathbf{1}+2\mu \mathit{\epsilon }-\kappa (T-T_0)\mathbf{1}):{\mathrm{\nabla }}^s\hat{\mathit{v}}\text{ d}\mathrm{\Omega }=W_{ext}(\hat{\mathit{v}})\mathrm{\forall }\hat{\mathit{v}}\in V_U(2)$$

where VU$V_U$ is the displacement FunctionSpace and Wext$W_{ext}$ the linear functional corresponding to the work of external forces.
其中 VU$V_U$ 是位移 FunctionSpace， Wext$W_{ext}$ 是对应外力功的线性泛函。

The solution of the coupled problem at t=tn+1$t=t_{n+1}$ is now (un+1,Tn+1)=(u,T)∈VU×VT$({\mathit{u}}_{n+1},T_{n+1})=(\mathit{u},T)\in V_U\times V_T$ verifying (1)$(1)$ and (2)$(2)$. These two forms are implemented below with zero right-hand sides (zero Neumann BCs for both problems here). One slight modification is that the temperature unknown T$T$ is replaced by the temperature variation Θ=T−T0$\mathrm{\Theta }=T-T_0$ which appears naturally in the stress constitutive relation.
在 t=tn+1$t=t_{n+1}$ 处的耦合问题的解现在是 (un+1,Tn+1)=(u,T)∈VU×VT$({\mathit{u}}_{n+1},T_{n+1})=(\mathit{u},T)\in V_U\times V_T$ ，它满足 (1)$(1)$ 和 (2)$(2)$ 。这两种形式在下面实现，并且右侧为零（这里两个问题的 Neumann 边界条件都为零）。一个微小的修改是温度未知数 T$T$ 被温度变化 Θ=T−T0$\mathrm{\Theta }=T-T_0$ 所替代，后者在应力本构关系中自然出现。

Note: We will later make use of the lhs and rhs functions to extract the corresponding bilinear and linear forms.
注意：稍后我们将使用 lhs 和 rhs 函数来提取相应的双线性形式和线性形式。

U_ = TestFunction(V)
(u_, Theta_) = split(U_)
dU = TrialFunction(V)
(du, dTheta) = split(dU)
Uold = Function(V)
(uold, Thetaold) = split(Uold)


def eps(v):
    return sym(grad(v))


def sigma(v, Theta):
    return (lmbda*tr(eps(v)) - kappa*Theta)*Identity(2) + 2*mu*eps(v)


dt = Constant(0.)
mech_form = inner(sigma(du, dTheta), eps(u_))*dx
therm_form = (cV*(dTheta-Thetaold)/dt*Theta_ +
              kappa*T0*tr(eps(du-uold))/dt*Theta_ +
              dot(k*grad(dTheta), grad(Theta_)))*dx
form = mech_form + therm_form

Resolution  分辨率 ¶

The problem is now solved by looping over time increments. Because of the typical exponential time variation of temperature evolution of the heat equation, time steps are discretized on a non-uniform (logarithmic) scale. Δt$\mathrm{\Delta }t$ is therefore updated at each time step. Note that since we work in terms of temperature variation and not absolute temperature all fields can be initialized to zero, otherwise T$T$ would have needed to be initialized to the reference temperature T0$T_0$.
现在通过循环时间增量来求解该问题。由于热传导方程中温度演化的典型指数时间变化，时间步长在非均匀（对数）尺度上进行离散化。因此， Δt$\mathrm{\Delta }t$ 在每一步时间更新。请注意，由于我们基于温度变化而非绝对温度进行计算，所有场都可以初始化为零，否则 T$T$ 将需要初始化为参考温度 T0$T_0$ 。

Nincr = 100
t = np.logspace(1, 4, Nincr+1)
Nx = 100
x = np.linspace(R, L, Nx)
T_res = np.zeros((Nx, Nincr+1))
U = Function(V)
for (i, dti) in enumerate(np.diff(t)):
    print("Increment " + str(i+1))
    dt.assign(dti)
    solve(lhs(form) == rhs(form), U, bcs)
    Uold.assign(U)
    T_res[:, i+1] = [U(xi, 0.)[2] for xi in x]

At each time increment, the variation of the temperature increase Θ$\mathrm{\Theta }$ along a line (x,y=0)$(x,y=0)$ is saved in the T_res array. This evolution is plotted below. As expected, the temperature gradually increases over time, reaching eventually a uniform value of +10∘C$+{10}^{\circ }\text{C}$ over infinitely long waiting time.
在每个时间增量中，沿线 (x,y=0)$(x,y=0)$ 的温度增量 Θ$\mathrm{\Theta }$ 的变化被保存在 T_res 数组中。这一演变过程绘制在下图中。正如预期的那样，温度随着时间的推移逐渐升高，最终在无限长的等待时间内达到均匀值 +10∘C$+{10}^{\circ }\text{C}$ 。

plt.figure()
plt.plot(x, T_res[:, 1::Nincr//10])
plt.xlabel("$x$-coordinate along $y=0$")
plt.ylabel("Temperature variation $\Theta$")
plt.legend(["$t={:.0f}$".format(ti) for ti in t[1::Nincr//10]], ncol=2)
plt.show()

The final stresses and temperature variation are plotted using matplotlib. Note that we can add colorbar and change the plot axis limit.
使用 matplotlib 绘制最终应力和温度变化。请注意，我们可以添加色条并更改绘图轴的限制。

u, Theta = split(U)

plt.figure()
p = plot(sigma(u, Theta)[1, 1], title="$\sigma_{yy}$ stress near the hole")
plt.xlim((0, 3*R))
plt.ylim((0, 3*R))
plt.colorbar(p)
plt.show()

plt.figure()
p = plot(Theta, title="Temperature variation")
plt.xlim((0, 3*R))
plt.ylim((0, 3*R))
plt.colorbar(p)
plt.show()

References

[FAR91]: Farhat, C., Park, K. C., & Dubois-Pelerin, Y. (1991). An unconditionally stable staggered algorithm for transient finite element analysis of coupled thermoelastic problems. Computer Methods in Applied Mechanics and Engineering, 85(3), 349-365. https://doi.org/10.1016/0045-7825(91)90102-C