python有限元_sfepy: python有限元分析模块介绍-【3】交互式方式进行有限元分析

本文是简单的搬运加意译，原文：Tutorial - SfePy 2019.3+git.77e93336de469766c0aa904f07e28b0b17d8a467 documentation​sfepy.org

对于sfepy来说，上一节说的用sfepy-run命令行方法运行program-file的方法，可以说是非常简单，但是控制性不足，下面介绍用ipython交互式环境，一步步的进行分析。

【1】 问题描述

我们要解决线弹性(linear elasticity)问题，其二维描述方程如下：

其中 应力(stress)：

和

称为拉梅常数(Lamé’s constants)

其中

是拉梅常数的一阶常数，表示材料的压缩性，等价与 体弹性模量(modulus of volume elasticity) 或者 杨氏模量(Young's modulus)

是拉梅常数的二阶常数，在线弹性力学中，表示剪切模量，shear modulus

是应变(strain)，

是 体积力 (volume forces)

令:

于是式(1)可以改写为 通用式(general form)：

这里

是测试函数(test function)，

、

都属于 suitable function space。

其中

是刚度矩阵。

【2】用sfepy进行分析：

执行ipython，进入交互式界面：

【2.1】导入相关的模块：

In [1]: import numpy as nm

In [2]: from sfepy.discrete.fem import Mesh, FEDomain, Field

【2.2】导入mesh文件，这定义了分析的主区域( main domain)

:

In [3]: from sfepy import data_dir #sfepy模块的目录

In [4]: mesh = Mesh.from_file( data_dir + '/meshes/2d/rectangle_tri.mesh')

创建主区域：

In [5]: domain = FEDomain('domain', mesh)

【2.3】创建区域region，用于加载初始条件，载荷，边界条件等：

In [6]: min_x, max_x = domain.get_mesh_bounding_box()[:, 0]

In [7]: eps = 1e-8 * (max_x - min_x)

In [8]: omega = domain.create_region('Omega', 'all')

In [9]: gamma1 = domain.create_region('Gamma1',

...: 'vertices in x

...: 'facet')

In [10]:gamma2 = domain.create_region('Gamma2',

...: 'vertices in x >%.10f' % (max_x - eps),

...: 'facet')

上面min_x和max_x是几何模型的最小最大x坐标，接着创建了3个区域： (标记：Omega)：针对所有单元

(标记：Gamma1)：所有节点的 x

(标记：Gamma2)：所有节点的 x>max_x-eps的面片(facet)单元

这里

(标记为：eps)是一个非常小的范围：1e-8*(max_x-min_x)，是加载边界条件的区域宽度，如下图显示：

【2.4】定义实际有限元近似(actual finite element approximation)：

In [11]: field = Field.from_args('fu', nm.float64, 'vector', omega, approx_order=2)

两个变量：未知变量unknown variable:

测试变量test variable：

可以用下面的语句来创建：

In [12]: from sfepy.discrete import (FieldVariable, Material, Integral, Function,

....: Equation, Equations, Problem)

In [13]: u = FieldVariable('u', 'unknown', field)

In [14]: v = FieldVariable('v', 'test', field, primary_var_name='u')

注意要先定义u，再定义v。

【2.5】定义材料参数：

In [15]: from sfepy.mechanics.matcoefs import stiffness_from_lame

In [16]: m = Material('m', D=stiffness_from_lame(dim=2, lam=1.0, mu=1.0))

In [17]: f = Material('f', val=[[0.02], [0.01]])

这里是通过拉梅常数来定义材料的刚度系数

(标记为: D)

此外，把 volume forces也定义为材料参数了，可以直接叠加到每一个单元上，其值为[0.02, 0.01]的转置。

【2.6】 定义各个积分项：

In [18]: integral = Integral('i', order=3)

In [19]: from sfepy.terms import Term

In [20]: t1 = Term.new('dw_lin_elastic(m.D, v, u)',

....: integral, omega, m=m, v=v, u=u)

In [21]: t2 = Term.new('dw_volume_lvf(f.val, v)',

....: integral, omega, f=f, v=v)

In [22]: eq = Equation('balance', t1 + t2)

In [23]: eqs = Equations([eq])

这里：

，

，最后建立方程，组成方程组(这里只有一个方程)。

【2.7】设置边界条件：

In [24]: from sfepy.discrete.conditions import Conditions, EssentialBC

In [25]: fix_u = EssentialBC('fix_u', gamma1, {'u.all' : 0.0})

In [26]: def shift_u_fun(ts, coors, bc=None, problem=None, shift=0.0):

....: val = shift * coors[:,1]**2

....: return val

In [27]: bc_fun = Function('shift_u_fun', shift_u_fun,

....: extra_args={'shift' : 0.01})

In [28]: shift_u = EssentialBC('shift_u', gamma2, {'u.0' : bc_fun})

注意，在Gamma1的区域设置不能移动，在Gamma2的区域则设置一个分布力。

【2.8】设置求解器(solvers):

In [29]: from sfepy.base.base import IndexedStruct

In [30]: from sfepy.solvers.ls import ScipyDirect

In [31]: from sfepy.solvers.nls import Newton

In [32]: ls = ScipyDirect({})

In [33]: nls_status = IndexedStruct()

In [34]: nls = Newton({}, lin_solver=ls, status=nls_status)

线性求解器ls，和非线性求解器nls，后者用了牛顿迭代的算法。

【2.9】 创建problem：

In [35]: pb = Problem('elasticity', equations=eqs)

In [36]: pb.save_regions_as_groups('regions')

设置求解的范围是弹性(elasticity)，并把相关的regions的信息提取出来输出到regions.vtk文件中方便查看。

【2.10】查看regions信息：

In [37]: from sfepy.postprocess.viewer import Viewer

In [38]: view = Viewer('regions.vtk')

In [39]: view()

如下图：

上图中左边显示了Gamma1区域，中间显示了Gamma2区域，最右边显示Omega区域。

【2.11】求解

In [39]: pb.set_bcs(ebcs=Conditions([fix_u, shift_u]))

In [40]: pb.set_solver(nls)

In [41]: status = IndexedStruct()

In [42]: vec = pb.solve(status=status)

In [43]: print('Nonlinear solver status:\n', nls_status)

In [44]: print('Stationary solver status:\n', status)

In [45]: pb.save_state('linear_elasticity.vtk', vec)

In [46]: view = Viewer('linear_elasticity.vtk')

In [47]: view()

上面设置了边界条件，然后开始求解，结果保存在linear_elasticity.vtk文件中，显示如下：

让我们以更直观的对比显示来比较：

In [48]: view(vector_mode='warp_norm', rel_scaling=2,

...: is_scalar_bar=True, is_wireframe=True)

下一篇，我们介绍gmsh的使用方法。