本文是简单的搬运加意译,原文:Tutorial - SfePy 2019.3+git.77e93336de469766c0aa904f07e28b0b17d8a467 documentationsfepy.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的使用方法。