getfem++ 热弹性-电耦合算例 Python 实现

网上似乎没有好的译文，自己翻译了一下。原文戳这里。

问题设定

考虑一个厚度为 $\epsilon$ 的二维平板，受到外力、电势和加热作用。令 $\Omega \subset \mathrm{I}{\mathrm{R}}^2$为平板计算域，记 $\theta :\Omega \to \mathrm{I}\mathrm{R}$ 为温度场（摄氏度），$V:\Omega \to \mathrm{I}\mathrm{R}$ 为电势场，$u:\Omega \to \mathrm{I}\mathrm{R}$ 为薄板的位移场。

热问题

假设板的侧面（上下左右）是绝热的，前后表面与空气存在热交换，热传导系数为 $D$。温度方程和边界条件为
$$\{\begin{array}{l}-\text{div}(\epsilon \kappa (\nabla \theta ))+2D(\theta -T_0)-\epsilon \sigma |\nabla V{|}^2=0\text{ in }\Omega ,\\ \kappa \nabla \theta \cdot n=0\text{ on }\partial \Omega ,\end{array}$$
其中热导率由 $\kappa$ 设定，$T_0$ 是空气温度，$n$ 是外法向量。$\sigma |\nabla V{|}^2$ 是与焦耳热对应的非线性耦合项，$\sigma$ 是电导率。

电势问题

假设右侧面和左侧面之间存在一个电势差 $0.1V$。其他表面为电绝缘。电势方程为
$$\{\begin{array}{l}-\text{div}(\epsilon \sigma (\nabla V))=0\text{ in }\Omega ,\\ \sigma \nabla V\cdot n=0\text{ on the top an bottom lateral faces},\\ V=0\text{ on the right lateral face},\\ V=0.1\text{ on the left lateral face},\end{array}$$
其中电导率依赖于温度，即
$$\sigma =\frac{1}{{\rho }_0(1+\alpha (\theta -T_0))},$$
其中 ${\rho }_0$和 $\alpha$ 均为电阻温度系数。

变形问题

考虑平板在右侧面力和热的影响下发生薄壁小变形。线弹性位移方程为
$$\{\begin{array}{l}-\text{div}(\bar{\sigma }(u))=0\text{ in }\Omega ,\\ \bar{\sigma }n=0\text{ on the top an bottom lateral faces},\\ \bar{\sigma }n=F\text{ on the right lateral face},\\ u=0\text{ on the left lateral face},\end{array}$$
其中 $F$ 是右侧面上的表面力，$\bar{\sigma }$ 是柯西应力张量，定义为
$$\bar{\sigma }(u)={\lambda }^{\ast }\text{div}(u)I+2\mu \bar{\epsilon }(u)+\beta (T_0-\theta )I,$$
线性应变张量为
$$\bar{\epsilon }(u)=(\nabla u+(\nabla u{)}^T)/2$$
${\lambda }^{\ast }$ 和 $\mu$ 为拉梅系数。$\beta (T_0-\theta )I$ 对应于热膨胀，$\beta ={\alpha }_{th}E/(1-2\nu )$， 其中 ${\alpha }_{th}$ 为热膨胀系数。

弱形式

基于伽辽金法的弱形式决定了我们要在代码中添加哪些项。该问题的弱形式为
$$\begin{gathered} \text{Find }\theta ,V,u\text{ with }V=0.1,u=0\text{ on the left face},V=0\text{ on the right face}, \\ {\int }_{\Omega }\epsilon \kappa \nabla \theta \cdot \nabla {\delta }_{\theta }+2D\theta {\delta }_{\theta }d\Omega ={\int }_{\Omega }(2D{T}_0+\epsilon \sigma |\nabla V{|}^2){\delta }_{\theta }d\Omega \text{ for all }{\delta }_{\theta }, \\ {\int }_{\Omega }\epsilon \sigma \nabla V\cdot \nabla {\delta }_V=0d\Omega \text{ for all }{\delta }_V\text{ satisfying }{\delta }_V=0\text{ on the left and right faces}, \\ {\int }_{\Omega }\bar{\sigma }(u):\bar{\epsilon }({\delta }_u)d\Omega ={\int }_{{\Gamma }_N}F\cdot {\delta }_{u}d\Gamma \text{ for all }{\delta }_u\text{ satisfying }{\delta }_u=0\text{ on the left face}, \end{gathered}$$
其中 ${\delta }_{\theta },{\delta }_V,{\delta }_u$ 分别代表 $\theta ,V,u$ 的试函数。${\Gamma }_N$ 表示施加表面力的右侧边界。

Python 实现

1. 导入

   import getfem as gf
   import numpy as np
   

2. 模型参数

   epsilon = 1.; E = 21E6; nu = 0.3;
   clambda = E*nu/((1+nu)*(1-2*nu));
   cmu = E/(2*(1+nu));
   clambdastar = 2*clambda*cmu/(clambda+2*cmu);
   F = 100E2; kappa = 4.; D = 10;
   air_temp = 20; alpha_th = 16.6E-6;
   T0 = 20; rho_0 = 1.754E-8;
   alpha = 0.0039;
   
   h = 2; elements_degree = 2;
   

3. 网格生成
   getfem++ 提供较弱的网格工具，因此此处采用自带工具。一般建议借助更专业的外部软件导入网格。以下代码中 h 代表网格尺寸，网格的 degree 等于2（代表单元的插值多项式次数）。

   mo1 = gf.MesherObject('rectangle', [0., 0.], [100., 25.])
   mo2 = gf.MesherObject('ball', [25., 12.5], 8.)
   mo3 = gf.MesherObject('ball', [50., 12.5], 8.)
   mo4 = gf.MesherObject('ball', [75., 12.5], 8.)
   mo5 = gf.MesherObject('union', mo2, mo3, mo4)
   mo  = gf.MesherObject('set minus', mo1, mo5)
   
   mesh = gf.Mesh('generate', mo, h, 2)
   

4. 边界选择

   由于边界的不同部分具有不同的边界条件，我们需要将边界的不同部分进行编号（孔洞边界上假设为绝热绝缘且应力自由）。我们要选择一些网格单元的面（face）并定义为 region，编号为1,2,3,4，依次对应右、左、上、下边界。这些边界编号将用于构建模型的块 (brick)。

   fb1 = mesh.outer_faces_in_box([1., 1.], [99., 24.])
   fb2 = mesh.outer_faces_with_direction([ 1., 0.], 0.01)
   fb3 = mesh.outer_faces_with_direction([-1., 0.], 0.01)
   fb4 = mesh.outer_faces_with_direction([0.,  1.], 0.01)
   fb5 = mesh.outer_faces_with_direction([0., -1.], 0.01)
   
   RIGHT_BOUND=1; LEFT_BOUND=2; TOP_BOUND=3; BOTTOM_BOUND=4; HOLE_BOUND=5;
   
   mesh.set_region( RIGHT_BOUND, fb2)
   mesh.set_region(  LEFT_BOUND, fb3)
   mesh.set_region(   TOP_BOUND, fb4)
   mesh.set_region(BOTTOM_BOUND, fb5)
   mesh.set_region(  HOLE_BOUND, fb1)
   mesh.region_subtract( RIGHT_BOUND, HOLE_BOUND)
   mesh.region_subtract(  LEFT_BOUND, HOLE_BOUND)
   mesh.region_subtract(   TOP_BOUND, HOLE_BOUND)
   mesh.region_subtract(BOTTOM_BOUND, HOLE_BOUND)
   

   注意：上面这段代码相比于 demo 代码，少了孔洞边界的定义。但demo 代码中的 outer_faces_in_ball 函数会报错称不存在（文档里可查到，原因不明）。可以把那部分代码改为：

   rh = 8 + 0.01*h
   fb6 = mesh.outer_faces_in_box([25-rh, 12.5-rh], [25+rh, 12.5+rh]) # Left hole boundary
   fb7 = mesh.outer_faces_in_box([50-rh, 12.5-rh], [50+rh, 12.5+rh]) # Center hole boundary
   fb8 = mesh.outer_faces_in_box([75-rh, 12.5-rh], [75+rh, 12.5+rh]) # Right hole boundary

   或者干脆把孔洞边界都注释掉。

5. 网格绘制

   mesh.export_to_vtk('mesh.vtk');
   

6. 有限元和积分方法定义
   getfem++的所谓“Fem (finite element method)”其实是指单元类型，而所谓“MeshFem (mesh of finite element method)”则是指由指定类型的单元组成的离散网格。以下将采取后一种表述以防混淆。现在我们选择经典单元（即连续拉格朗日单元），且单元阶数（degree）等于2（代表二阶等参元）。（这里有一份可选择的有限元列表。最常用的就是经典有限元。）我们定义三个离散网格。第一个向量离散网格mfu用于近似位移场。第二个标量离散网格mft用于近似温度场和电势场。原则上一个离散网格允许被用于近似任意个变量。第三个是经典不连续标量离散网格（即间断标量拉格朗日离散网格），可被用于插值某个变量的导数，例如插值Mises应力。最后定义一个积分方法。getfem++没有默认积分方法，所以这一步是必须的。积分阶数必须足够高以保证积分的方便。这里令阶数等于单元阶数平方就够了。

   mfu = gf.MeshFem(mesh, 2)
   mfu.set_classical_fem(elements_degree)
   mft = gf.MeshFem(mesh, 1)
   mft.set_classical_fem(elements_degree)
   mfvm = gf.MeshFem(mesh, 1)
   mfvm.set_classical_discontinuous_fem(elements_degree)
   mim = gf.MeshIm(mesh, elements_degree*2)
   

7. 模型定义
   model 对象包含了模型的变量（unknowns）、数据 (data) 和所谓的块(model bricks)。块是弱形式中的一个项，可以是线性或非线性（注意这个线性与线弹性概念不是一回事，是指微分方程的线性。），用于一个或多个变量。real表示模型变量都是实数。

   md=gf.Model('real');
   md.add_fem_variable('u', mfu)
   md.add_fem_variable('theta', mft)
   md.add_fem_variable('V', mft)
   

   7.1 薄壁弹性变形问题

   我们使用预定义的add_isotropic_linearized_elasticity_brick构建一个块，它代表以下的项：
   $${\int }_{\Omega }({\lambda }^{\ast }\text{div}(u)I+2\mu \bar{\epsilon }(u)):\bar{\epsilon }({\delta }_u)dx,$$
   除了预定义方法，也可以用“通用弱形式语言（GWFL, generic weak form language）”自己定义每个块。 add_isotropic_linearized_elasticity_brick需要拉梅系数对应的数据（数据在 getfem++ 里称为 data，是定义在一个离散网格上的长度等于单元数量的数组，数组元素的形状与 MeshFem 变量一致。）。拉梅系数的数据可以是常数或非常数的数组。

   然后构建耦合项
   $${\int }_{\Omega }(\beta \theta I):\bar{\epsilon }({\delta }_u)dx,$$
   由于没有预定义的块，所以我们用 GWFL 构建。该项为线性项，GWFL 构建方法为add_linear_term。GWFL 定义的表达式将会在高斯点上执行，并组装到总体矩阵中。

   最后用预定义方法add_Dirichlet_condition_with_multipliers和add_source_term_brick构建了狄利克雷边界条件块。关于狄利克雷边界条件的预定义方法有很多，详见这里。add_initialized_data里的 initialized 代表该数据属于当前时间步（静力学只有一步）的起点，在该步的迭代中保持不变。

   md.add_initialized_data('cmu', [cmu])
   md.add_initialized_data('clambdastar', [clambdastar])
   md.add_initialized_data('T0', [T0])
   md.add_isotropic_linearized_elasticity_brick(mim, 'u', 'clambdastar', 'cmu')
   
   md.add_Dirichlet_condition_with_multipliers(mim, 'u', elements_degree-1, LEFT_BOUND)
   md.add_initialized_data('Fdata', [F*epsilon, 0])
   md.add_source_term_brick(mim, 'u', 'Fdata', RIGHT_BOUND)
   
   md.add_initialized_data('beta', [alpha_th*E/(1-2*nu)])
   md.add_linear_term(mim, 'beta*(T0-theta)*Div_Test_u')
   

   7.2 电势问题

   方法同上。

   sigmaeps = '(eps/(rho_0*(1+alpha*(theta-T0))))'
   md.add_initialized_data('eps', [epsilon])
   md.add_initialized_data('rho_0', [rho_0])
   md.add_initialized_data('alpha', [alpha])
   md.add_nonlinear_term(mim, sigmaeps+'*(Grad_V.Grad_Test_V)')
   md.add_Dirichlet_condition_with_multipliers(mim, 'V', elements_degree-1, RIGHT_BOUND)
   md.add_initialized_data('DdataV', [0.1])
   md.add_Dirichlet_condition_with_multipliers(mim, 'V', elements_degree-1, LEFT_BOUND, 'DdataV')
   

   7.3 热问题

   方法同上。

   md.add_initialized_data('kaeps', [kappa*epsilon])
   md.add_generic_elliptic_brick(mim, 'theta', 'kaeps')
   md.add_initialized_data('D2', [D*2])
   md.add_initialized_data('D2airt', [air_temp*D*2])
   md.add_mass_brick(mim, 'theta', 'D2')
   md.add_source_term_brick(mim, 'theta', 'D2airt')
   
   md.add_nonlinear_term(mim, '-'+sigmaeps+'*Norm_sqr(Grad_V)*Test_theta')
   

8. model 求解

   8.1. 强耦合求解

   md.solve('max_res', 1E-9, 'max_iter', 100, 'noisy')
   

   由于耦合是非线性的，所以求解器用了牛顿迭代法，迭代了4次。

   8.2. 弱耦合求解

   也可以先求解热电问题再求解变形问题，这样就使得温度和电势与变形不同步了，所以称为弱耦合。这通过 disable_variable 和 enable_variable 来实现。

   md.disable_variable('u')
   md.solve('max_res', 1E-9, 'max_iter', 100, 'noisy')
   md.enable_variable('u')
   md.disable_variable('theta')
   md.disable_variable('V')
   md.solve('max_res', 1E-9, 'max_iter', 100, 'noisy')
   

9. 导出

   数据切片和导出详见这里。

   U = md.variable('u')
   V = md.variable('V')
   THETA = md.variable('theta')
   VM = md.compute_isotropic_linearized_Von_Mises_or_Tresca('u', 'clambdastar', 'cmu', mfvm)
   CO = np.reshape(md.interpolation('-'+sigmaeps+'*Grad_V', mfvm), (2, mfvm.nbdof()), 'F')
   
   mfvm.export_to_vtk('displacement_with_von_mises.vtk', mfvm,
                     VM, 'Von Mises Stresses', mfu, U, 'Displacements')
   print ('You can view solutions with for instance:')
   print ('mayavi2 -d displacement_with_von_mises.vtk -f WarpVector -m Surface')
   mft.export_to_vtk('temperature.vtk', mft, THETA, 'Temperature')
   print ('mayavi2 -d temperature.vtk -f WarpScalar -m Surface')
   mft.export_to_vtk('electric_potential.vtk', mft, V, 'Electric potential')
   print ('mayavi2 -d electric_potential.vtk -f WarpScalar -m Surface')