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

最新推荐文章于 2021-07-18 23:12:38 发布

鼠标君

最新推荐文章于 2021-07-18 23:12:38 发布

阅读量812

点赞数

分类专栏：有限元编程文章标签：有限元

本文链接：https://blog.csdn.net/iamzhtr/article/details/109213821

版权

编程同时被 2 个专栏收录

10 篇文章 0 订阅

订阅专栏

有限元

3 篇文章 3 订阅

订阅专栏

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

问题设定

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

在这里插入图片描述

热问题

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

电势问题

假设右侧面和左侧面之间存在一个电势差 $0.1 V$ 。其他表面为电绝缘。电势方程为
$\left\{\begin{array}{l} -\text{div}(\varepsilon\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} \right.$
其中电导率依赖于温度，即
$\sigma = \dfrac{1}{\rho_0(1+\alpha(\theta - T_0))},$
其中 $\rho_0$ 和 $\alpha$ 均为电阻温度系数。

变形问题

考虑平板在右侧面力和热的影响下发生薄壁小变形。线弹性位移方程为
$\left\{\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} \right.$
其中 $F$ 是右侧面上的表面力， $\bar{\sigma}$ 是柯西应力张量，定义为
$\bar{\sigma}(u) = \lambda^* \text{div}(u) I + 2\mu \bar{\varepsilon}(u) + \beta(T_0-\theta) I,$
线性应变张量为
$\bar{\varepsilon}(u) = (\nabla u + (\nabla u)^T)/2$
$\lambda^*$ 和 $\mu$ 为拉梅系数。 $\beta(T_0-\theta) I$ 对应于热膨胀， $\beta = \alpha_{th} E/(1-2\nu)$ ，其中 $\alpha_{th}$ 为热膨胀系数。

弱形式

基于伽辽金法的弱形式决定了我们要在代码中添加哪些项。该问题的弱形式为
$\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} \varepsilon\kappa\nabla\theta\cdot\nabla\delta_{\theta} + 2D\theta\delta_{\theta}d\Omega = \int_{\Omega} (2DT_0 + \varepsilon\sigma|\nabla V|^2)\delta_{\theta} d\Omega ~~~\text{ for all } \delta_{\theta}, \\ \int_{\Omega} \varepsilon\sigma\nabla V\cdot\nabla\delta_V = 0 d\Omega ~~~ \text{ for all } \delta_V \text{ satisfying } \delta_V = 0 \text{ on the left and right faces}, \\ \int_{\Omega} \bar{\sigma}(u):\bar{\varepsilon}(\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},$
其中 $\delta_{\theta}, \delta_V, \delta_u$ 分别代表 $\theta, V, u$ 的试函数。 $\Gamma_N$ 表示施加表面力的右侧边界。

Python 实现

导入
```
import getfem as gf
import numpy as np
```

模型参数

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;

网格生成
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)

边界选择

由于边界的不同部分具有不同的边界条件，我们需要将边界的不同部分进行编号（孔洞边界上假设为绝热绝缘且应力自由）。我们要选择一些网格单元的面（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

或者干脆把孔洞边界都注释掉。
网格绘制
```
mesh.export_to_vtk('mesh.vtk');
```
有限元和积分方法定义
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)
```
模型定义
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^* \text{div}(u) I + 2\mu \bar{\varepsilon}(u)):\bar{\varepsilon}(\delta_u)dx,$
除了预定义方法，也可以用“通用弱形式语言（GWFL, generic weak form language）”自己定义每个块。 add_isotropic_linearized_elasticity_brick需要拉梅系数对应的数据（数据在 getfem++ 里称为 data，是定义在一个离散网格上的长度等于单元数量的数组，数组元素的形状与 MeshFem 变量一致。）。拉梅系数的数据可以是常数或非常数的数组。

然后构建耦合项
$\int_{\Omega} (\beta\theta I) :\bar{\varepsilon}(\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')
```
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')
```

导出

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

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')

在这里插入图片描述

鼠标君

关注

0
点赞
踩
3

收藏

觉得还不错? 一键收藏
0
评论
getfem++ 热弹性-电耦合算例 Python 实现

原文戳这里。问题设定考虑一个厚度为 ε\varepsilonε 的二维平板，受到外力、电势和加热作用。令 Ω⊂IR2\Omega \subset \rm I\hspace{-0.15em}R^2Ω⊂IR2为平板计算域，记 θ:Ω→IR\theta : \Omega \rightarrow \rm I\hspace{-0.15em}Rθ:Ω→IR 为温度场（摄氏度），V:Ω→IRV : \Omega \rightarrow \rm I\hspace{-0.15em}RV:Ω→IR 为电势场，u:Ω→IRu
复制链接

扫一扫