CN108491636B - 基于几何约束的弹性体网格变形方法 - Google Patents

Abstract

本发明提供一种基于几何约束的弹性体网格变形方法，包括以下步骤：步骤S100：根据待变形的弹性体网格受外力变形后，弹性体网格内部网格单元产生的弹性应变ε和网格单元的应力状态σ，构建含网格单元的弹性模量E的待变形弹性体网格变形模型；步骤S200：建立网格单元的几何约束E_G为E_G＝E_V+E_θ，得到约束弹性模量E'：E'＝E_G·E＝(E_V+E_θ)·E；步骤S300：将约束弹性模量E'代入待变形弹性体网格变形模型中得到约束方程，求解约束方程，得到具有几何约束的变形后的网格节点坐标，根据网格节点坐标得到变形后的弹性体网格。能充分利用网格单元的几何属性约束变形效果较差的网格单元，阻止整体变形过程中畸变的扩大，提高变形方法的鲁棒性。

Description

基于几何约束的弹性体网格变形方法

技术领域

本发明涉及工程设计技术领域，具体的涉及一种基于几何约束的弹性体网格变

形方法。

背景技术

在日益复杂的工程问题仿真中，经常需要求解包含运动边界的非定常流场问题，尤其随着计算流体力学的快速发展，网格变形方法已经成为支撑非定常流场仿真的关键技术。网格变形可以基于现有的网格节点信息和变形算法，根据动边界的变化计算网格内部节点的位移，不用改变网格拓扑关系就可得到变形网格。在工程设计领域，如自由表面振动、机翼气动弹性、流固耦合、飞行器气动外形优化设计等，用户可以利用网格变形技术对研究对象的计算网格重复利用，既能提高设计效率，又能减少更替网格带来的计算误差。

目前，研究较多的网格变形方法有两个类别：物理模型法和代数法(也称数学插值法)。其中物理模型法是基于物理模型的网格变形方法，通过求解模型建立的控制方程得到变形后的网格坐标，变形效果较好但是建模过程复杂，主要有：弹簧比拟法和弹性体法；而代数法是根据网格节点的坐标信息实现网格变形，效率较高但是稳定性不佳，主要有：径向基函数法和Delaunay背景网格插值法。综上各种网格变形方法，弹性体法把网格变形当成弹性介质力学问题来求解，其网格变形能力、网格质量和局部控制效果都要优于弹簧比拟法、径向基函数法和Delaunay背景网格插值法等其它方法。但是在运动边界的大变形问题中，弹性体网格变形方法的鲁棒性较差，特别是将其运用于运动边界附近的网格单元时，常会导致非法的“负体积”单元出现，即网格的拓扑结构发生破坏，导致网格变形过早而失败。此外，现有弹性体法虽然能将网格变形进行到较大程度，但是个别单元仍会因过度挤压或拉伸而发生严重畸变，导致变形网格无法满足数值计算的稳定性要求。

发明内容

本发明的目的在于提供一种基于几何约束的弹性体网格变形方法，该发明解决了现有弹性体网格变形方法鲁棒性差、大变形时易导致个别单元严重畸变，影响计算稳定性的技术问题。

参见图1，本发明提供的基于几何约束的弹性体网格变形方法，包括以下步骤：

步骤S100：根据待变形的弹性体网格受外力变形后，弹性体网格内部网格单元产生的弹性应变ε和网格单元的应力状态σ，构建含网格单元的弹性模量E的待变形弹性体网格变形模型；

本文中弹性体网格包括多个相互连接的网格单元。各网格单元的端点称为网格节点。

步骤S200：建立网格单元的尺寸约束

其中，V为二维网格单元的面积或三维网格单元的体积，建立网格单元的形状约束E_θ，

其中，θ_min为二维三角形网格单元的最小内角或三维四面体网格单元的最小二面角，取网格单元的几何约束E_G为E_G＝E_V+E_θ，得到约束弹性模量E'：E'＝E_G·E＝(E_V+E_θ)·E；

步骤S300：将约束弹性模量E'代入待变形弹性体网格变形模型中得到约束方程，求解约束方程，得到具有几何约束的变形后的网格节点坐标，根据网格节点坐标得到变形后的弹性体网格。

本发明提供的方法通过对待变形弹性体网格变形模型中的弹性模量E重新定义，从而有效避免了网格的拓扑结构在运动边界附近发生破坏，提高了网格变形的鲁棒性。避免网格畸变的发生。以上未详述的步骤按现有方法进行即可。

进一步地，步骤S100包括以下步骤：

步骤S110：解析待变形的弹性体网格中所有网格节点的坐标和所有网格单元的拓扑关系，定义弹性体网格内任意网格节点(x,y,z)在外力作用下产生的位移矢量U(x,y,z)＝(u,v,w)，其中u、v、w分别为网格节点(x,y,z)在

三个方向的位移分量，弹性体网格内部网格单元的弹性应力张量满足

步骤S120：网格单元满足线性运动定律：

其中，ε为网格单元的应变状态。本文中(·)^T均表示转秩；

步骤S130：根据广义胡克定律，可得：

σ＝λTr(ε)I+2με (2)

其中，λ和μ为代表弹性网格材料属性的拉梅常量，

E为弹性模量，ν为泊松比，Tr(ε)＝ε_x+ε_y+ε_z，σ为网格单元的应力状态，I表示单位向量；

步骤S140：根据公式(1)(2)得到弹性体网格变形模型：

其中，

E为弹性模量和ν为泊松比。

进一步的，步骤S300中求解约束方程的步骤，包括以下步骤：

步骤S310：将待变形的弹性体网格作为有限元网格，建立有限元方程求解约束方程；

步骤S320：建立弹性体网格变形的边界条件，边界条件包括受力边界条件和位移边界条件；

步骤S330：求解网格单元变形后的线性代数方程组，得到具有几何约束的变形网格节点坐标。

进一步地，求解线性代数方程组的方法为高斯-赛德尔迭代法。

参见图1具体的，本发明提供的基于几何约束的弹性体网格变形方法，包括以下步骤：

步骤一，构建待变形的弹性体网格变形方程：

(1)导入待变形的网格，解析待变形的网格所有节点的坐标信息和所有单元的拓扑关系。

(2)定义弹性体内任意网格节点(x,y,z)的位移矢量U(x,y,z)＝(u,v,w)，其中u、v、w分别为网格节点(x,y,z)在

三个方向的位移分量。

根据弹性力学的基本原理，如果对弹性体网格施加一个外力，网格内部的所有节点会产生弹性应力，进而发生弹性应变，导致网格节点发生位移变化，从而网格系统再次达到平衡状态。

(3)定义在外力的作用下，弹性体网格产生的弹性应力张量满足：

其中，σ为网格单元的应力状态，σ_x，σ_y和σ_z分别为正应力沿笛卡尔坐标系的分量，τ_xy、τ_zy和τ_xz分别为剪应力沿笛卡尔坐标系的分量。本文中

均表示散度。

(4)在弹性应力的作用下，弹性体网格产生的弹性应变ε的状态表示为：

其中，ε为网格单元的应变状态，ε_x，ε_y和ε_z分别为正应变沿笛卡尔坐标系的分量，γ_xy、γ_zy和γ_xz分别为剪应变沿笛卡尔坐标系的分量。

在弹性应力的作用下，弹性体网格会发生弹性应变，从而得到关于应变ε的上述表述。

(5)结合位移矢量U和弹性应变ε的表达式，网格单元满足线性运动定律：

(6)在弹性体变形中，应力和应变满足广义胡克定律，整理后有：

σ＝λTr(ε)I+2με (7)

其中，λ和μ为代表弹性网格材料属性的拉梅常量，Tr(ε)＝ε_x+ε_y+ε_z，σ为网格单元的应力状态，I表示单位向量。

(7)综合上述弹性体网格变形建模过程，可得网格节点位移矢量U的表达式：

其中，拉梅常量λ和μ可用弹性模量E和泊松比ν来表示，即

在弹性体网格变形方法中，弹性模量E和泊松比ν扮演着重要角色，它们分别控制单元的刚度和网格可压缩程度，与弹簧比拟法中的刚度系数K类似。具体而言，E为胡克定律中比例常数，E值越大材料刚性越强，网格单元在变形过程中的变化程度就越小，而ν通常取[-1,0.5]中的常数。

步骤二，引入网格单元的几何约束。

(1)首先针对运动边界附近的网格单元容易出现“负体积”单元的问题，根据网格单元的面积或体积信息，建立网格单元的尺寸约束E_V，防止网格变形过早失败，表达式为

其中，V为二维网格单元的面积或三维网格单元的体积。

(2)然后针对单元因过度挤压或拉伸而发生严重畸变的问题，根据网格单元的最小内角或最小二面角信息，建立网格单元的形状约束E_θ，防止网格单元质量过差而失效，其表达式为

其中，θ_min为二维三角形网格单元的最小内角或三维四面体网格单元的最小二面角。

(3)整合上述建立的尺寸约束E_V和形状约束E_θ为网格单元的几何约束E_G，其表达式为

E_G＝E_V+E_θ (11)

(4)最后在弹性体网格变形模型中引入几何约束E_G，通过约束每个网格单元的弹性模量E，得到基于几何约束的约束弹性模量E'，其表达式为

E'＝E_G·E＝(E_V+E_θ)·E (12)

当网格单元的尺寸(面积、体积)V→0时，尺寸约束E_V→∞，进而弹性模量E'变大，使得网格单元不易变形，保证了该单元不会过早成为“负体积”的非法单元。

当网格单元的形状(最小内角、最小二面角)θ_min→0°时，形状约束E_θ→∞，进而弹性模量E'变大使得网格单元不易变形，避免了该单元过度遭受挤压而畸变严重，在一定程度上增强了弹性体网格变形的鲁棒性。

步骤三，弹性体网格变形方程求解。

(1)将待变形的弹性体网格看作有限元网格，建立有限元方程求解网格节点位移变化的方程(即约束方程，即

其中，

E'＝E_G·E＝(E_V+E_θ)·E。)

(2)建立弹性体网格变形的边界条件，包括受力边界条件和位移边界条件。

(3)将边界条件引入建立的有限元方程，即联立网格变形的有限元方程和边界条件，求解网格变形的线性代数方程组，得到更新后的网格节点坐标。

用求解线性代数方程组的方法求解即可，比如高斯-赛德尔迭代法。

参见图2，本发明的另一方面还提供了一种基于几何约束的弹性体网格变形装置，包括：

应变应力模块100，用于根据待变形的弹性体网格受外力变形后，弹性体网格内部网格单元产生的弹性应变ε和网格单元的应力状态σ，构建含网格单元的弹性模量E的待变形弹性体网格变形模型；

几何约束模块200，用于立网格单元的尺寸约束

其中，V为二维网格单元的面积或三维网格单元的体积，建立网格单元的形状约束E_θ，

其中，θ_min为二维三角形网格单元的最小内角或三维四面体网格单元的最小二面角，取网格单元的几何约束E_G为E_G＝E_V+E_θ，得到约束弹性模量E'：E'＝E_G·E＝(E_V+E_θ)·E；

求解模块300，用于将约束弹性模量E'代入待变形弹性体网格变形模型中得到约束方程，求解约束方程，得到具有几何约束的变形后的网格节点坐标，根据网格节点坐标得到变形后的弹性体网格。

进一步地，应变应力模块100，包括：

拓扑模块，用于解析待变形的弹性体网格中所有网格节点的坐标和所有网格单元的拓扑关系，定义弹性体网格内任意网格节点(x,y,z)在外力作用下产生的位移矢量U(x,y,z)＝(u,v,w)，其中u、v、w分别为网格节点(x,y,z)在

三个方向的位移分量，弹性体网格内部网格单元的弹性应力张量满足

线性运动模块，用于网格单元满足线性运动定律：

其中，ε为网格单元的应变状态，

ε_x，ε_y和ε_z分别为正应变沿笛卡尔坐标系的分量，γ_xy、γ_zy和γ_xz分别为剪应变沿笛卡尔坐标系的分量；

应力状态模块，用于根据广义胡克定律，可得：

σ＝λTr(ε)I+2με (2)

其中，λ和μ为代表弹性网格材料属性的拉梅常量，

E为弹性模量，ν为泊松比，Tr(ε)＝ε_x+ε_y+ε_z，σ为网格单元的应力状态，I表示单位向量；

变形模型模块，用于根据公式(1)和(2)得到弹性体网格变形模型：

其中，

E为弹性模量和ν为泊松比。

进一步地，求解模块300，包括：

有限元模块，用于将待变形的弹性体网格作为有限元网格，建立有限元方程求解约束方程；

边界条件模块，用于建立弹性体网格变形的边界条件，边界条件包括受力边界条件和位移边界条件；

代数求解模块，用于求解网格单元变形后的线性代数方程组，得到具有几何约束的变形网格节点坐标。

进一步地，求解线性代数方程组的方法为高斯-赛德尔迭代法。

相对现有技术具有的优点：

本发明提供的基于几何约束的弹性体网格变形方法，通过在所建立的待变形弹性体网格变形模型中的弹性模量E的定义进行改进，提高了所得弹性体网格变形方程的鲁棒性，使其具有广泛的适用性，能有效提高大变形时网格的变形能力和网格质量，为求解包含运动边界的非定常流场问题提高更加有力的支撑。

本发明提供的基于几何约束的弹性体网格变形方法，针对弹性体网格变形在运动边界的大变形问题中，1)运动边界附近的网格单元容易出现非法的“负体积”单元破坏网格拓扑结构的缺点；2)个别网格单元容易过度挤压或拉伸而发生严重畸变的缺点，分别引入了网格单元的尺寸约束和形状约束，并整合为几何约束，既增强了网格变形的最大变形能力，也保证了变形过程中较好的稳定性，较好解决了上述问题。

本发明提供的基于几何约束的弹性体网格变形方法，利用网格单元的几何属性(如二维问题中的单元面积、最小内角，三维问题中的单元体积、最小二面角等)来约束较差的网格单元，阻止其在整体变形过程中向畸变的方向继续变形，从而提高该方法鲁棒性。

本发明提供的基于几何约束的弹性体网格变形装置，能充分利用网格单元的几何属性约束变形效果较差的网格单元，阻止整体变形过程中畸变的扩大，提高变形方法的鲁棒性。

具体请参考根据本发明的基于几何约束的弹性体网格变形方法提出的各种实施例的如下描述，将使得本发明的上述和其他方面显而易见。

附图说明

图1为本发明提供的基于几何约束的弹性体网格变形方法流程示意框图；

图2为本发明提供的基于几何约束的弹性体网格变形装置结构示意图；

图3为本发明提供的优选实施例中所处理的初始网格模型示意图，其中(a)为网格模型的全局网格，(b)为网格模型的运动边界表面网格；

图4为本发明提供的优选实施例在平移运动中的网格质量随平移时间变化曲线示意图，其中(a)为平均网格质量的变化曲线，(b)为最小网格质量的变化曲线；

图5为本发明提供的优选实施例在旋转运动中的网格质量随旋转角度变化曲线示意图，其中(a)为平均网格质量的变化曲线，(b)为最小网格质量的变化曲线；

图6为本发明提供的优选实施例在伸展变形中的网格质量随伸展长度变化曲线示意图，其中(a)为平均网格质量的变化曲线，(b)为最小网格质量的变化曲线；

图7为在最大伸展长度时，本发明提供的优选实施例中经典方法和改进方法所得伸展变形网格结果对比示意图，其中(a)为现有弹性体网格变形方法所得结果，(b)为本发明优选实施例所得结果。

具体实施方式

构成本发明的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。

下面将结合一个具体的实施例，对本发明的技术方案进行清楚、完整的描述，显然所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的其他实施例，都属于本发明保护的范围。

以下是发明人给出的一个具体实施例，采用一般意义的三维非结构四面体网格单元，包括运动边界的平移、旋转和伸展变形三种运动形式。该实施例具体通过下列步骤来进行：

步骤一，建立弹性体网格变形模型。

(1)导入待变形的网格，解析网格所有节点的坐标信息和所有单元的拓扑关系，本发明实施例的初始网格模型如图3所示；

(2)定义任意网格节点(x,y,z)的位移矢量为U(x,y,z)＝(u,v,w)，其中u、v、w分别为网格节点(x,y,z)在

三个方向的位移分量。

(3)假设弹性体网格因为外力的作用发生变形，其内部网格单元的弹性应力张量满足守恒关系，即

(4)在弹性应力的作用下，弹性体网格单元发生弹性应变ε；

(5)弹性体网格单元满足线性运动定律，

(6)弹性体网格单元的应力和应变满足广义胡克定律，σ＝λTr(ε)I+2με (2)；

(7)建立待变形弹性体网格的变形模型，

(8)将拉梅常量λ和μ用弹性模量E和泊松比ν来表示，即

(9)本算例中，弹性模量E>0，泊松比ν＝0.3。

步骤二，引入网格单元的几何约束。

(1)根据网格单元的面积或体积信息，建立网格单元的尺寸约束

(2)根据网格单元的最小内角或最小二面角信息，建立网格单元的形状约束

(3)整合尺寸约束E_V和形状约束E_θ为几何约束E_G，E_G＝E_V+E_θ；

(4)在弹性模量E中引入网格单元的几何约束E_G，经过改进后的弹性模量为E'＝E_G·E＝(E_V+E_θ)·E，至此建立了基于几何约束的弹性体网格变形模型。

步骤三，弹性体网格变形方程求解。

(1)将要变形的弹性体网格看作有限元网格，建立有限元方法求解网格节点位移变化的方程。

(2)建立弹性体网格变形的位移边界条件，包括运动边界的平移、旋转和伸展变形三种运动形式。

(3)将位移边界条件引入建立的有限元方程，采用高斯-赛德尔迭代法求解网格变形的线性代数方程组，得到更新后的网格节点坐标，完成弹性体网格变形。

本发明包括以上三个步骤：建立弹性体网格变形模型、引入网格单元的几何约束和弹性体网格变形方程的求解。此外，为方便对比引入几何约束前后的变形网格鲁棒性，额外建立了网格质量参数，通过对比说明引入几何约束后的改进效果相对现有技术的技术效果。

步骤四，建立网格质量参数。

(1)针对本实施例的三维非结构四面体网格单元，建立如下所示的网格质量参数

其中，λ＝V/ξ,V为四面体单元体积，ξ为以该单元平均表面积构建的正三角形组成的等边四面体体积，l_i为四面体的任意边长，S_i为四面体的任意平面面积。

由上式可知，f_size-shape∈(0,1]，当网格单元质量较好时f_size-shape→1，网格单元质量较差时f_size-shape→0。

(2)进而，以每个单元的综合参数f_size-shape为基础，从网格质量的平均水平和最小极值两个层次，建立平均网格质量f_mean和最小网格质量f_min，其表达式为

其中，Ne为网格单元数目。

步骤五，对比引入几何约束前后的变形网格鲁棒性。

(1)由步骤四建立的网格质量参数可知，f_mean和f_min的数值越大，网格的鲁棒性越强，可用于量化评比变形网格的鲁棒性。

(2)本文中为简化表述，将现有的弹性体网格变形方法称为经典方法，将本发明提供的基于几何约束的弹性体网格变形方法称为改进方法。

(3)针对经典方法和改进方法，分别完成实施例由初始网格模型在三种不同方式(平移运动、旋转运动和伸展变形)下的网格变形，基于网格质量参数对比变形网格的鲁棒性。

图4为本发明实施例中，网格模型由初始状态沿坐标轴

三个方向，以2v_x＝v_y＝v_z＝10m/s的速度作0.5s的匀速平移运动，其网格质量随平移运动时间t的变化结果，其中(a)为平均网格质量，(b)为最小网格质量。基本可以看出，改进方法的平均网格质量和最小网格质量的下降速率小于经典方法，尤其对最小网格质量改进效果明显。

图5为本发明实施例中，网格模型由初始状态绕OY轴，以ω＝0.01745rad/s的角速度作最大角度φ_max＝60°的旋转运动，其网格质量随旋转角度φ的变化结果，其中(a)为平均网格质量，(b)为最小网格质量。基本可以看出，改进方法的平均网格质量和最小网格质量的下降速率小于经典方法，尤其对最小网格质量改进效果明显。

图6为本发明实施例中，网格模型由初始状态沿

方向，以ux＝10m/s的速度作最大伸展量为Lmax＝6.0m的匀速伸展变形，其网格质量随远端伸展长度L的变化结果，其中(a)为平均网格质量，(b)为最小网格质量。基本可以看出，改进方法的最小网格质量的下降速率小于经典方法，但是平均网格质量的下降速率略高于经典方法。

图7为本发明实施例在伸展变形运动中，经典方法和改进方法在最大伸展长度时的网格对比图，其中(a)为经典方法结果，(b)为改进方法结果。基本可以看出，经典方法出现了“负体积”的非法单元，网格拓扑结构已经破坏而变形失败，但是改进方法保持了完好的网格结构，具有更强的网格变形鲁棒性。

本领域技术人员将清楚本发明的范围不限制于以上讨论的示例，有可能对其进行若干改变和修改，而不脱离所附权利要求书限定的本发明的范围。尽管己经在附图和说明书中详细图示和描述了本发明，但这样的说明和描述仅是说明或示意性的，而非限制性的。本发明并不限于所公开的实施例。

通过对附图，说明书和权利要求书的研究，在实施本发明时本领域技术人员可以理解和实现所公开的实施例的变形。在权利要求书中，术语“包括”不排除其他步骤或元素，而不定冠词“一个”或“一种”不排除多个。在彼此不同的从属权利要求中引用的某些措施的事实不意味着这些措施的组合不能被有利地使用。权利要求书中的任何参考标记不构成对本发明的范围的限制。。

Claims (8)

1.一种基于几何约束的弹性体网格变形方法，其特征在于，包括以下步骤：

步骤S100：根据待变形的弹性体网格受外力变形后，所述弹性体网格内部网格单元产生的弹性应变ε和网格单元的应力状态σ，构建含网格单元的弹性模量E的待变形弹性体网格变形模型；所述待变形的弹性体网格是从研究对象的运动边界提取的，用于解决研究对象的运动边界的非定常流场问题；所述研究对象包括：自由表面、机翼以及飞行器气动外形；

步骤S200：建立所述网格单元的尺寸约束

其中，V为二维网格单元的面积或三维网格单元的体积，建立所述网格单元的形状约束E_θ，

其中，θ_min为二维三角形网格单元的最小内角或三维四面体网格单元的最小二面角，取网格单元的几何约束E_G为E_G＝E_V+E_θ，得到约束弹性模量E'：E'＝E_G·E＝(E_V+E_θ)·E；

步骤S300：将约束弹性模量E'代入待变形弹性体网格变形模型中得到约束方程，求解约束方程，得到具有几何约束的变形后的网格节点坐标，根据所述网格节点坐标得到变形后的所述弹性体网格。

2.根据权利要求1所述的基于几何约束的弹性体网格变形方法，其特征在于，所述步骤S100包括以下步骤：

步骤S110：解析待变形的所述弹性体网格中所有网格节点的坐标和所有网格单元的拓扑关系，定义所述弹性体网格内任意网格节点(x,y,z)在外力作用下产生的位移矢量U(x,y,z)＝(u,v,w)，其中u、v、w分别为网格节点(x,y,z)在

三个方向的位移分量，所述弹性体网格内部网格单元的弹性应力张量满足▽·σ＝0；

步骤S120：所述网格单元满足线性运动定律：

其中，ε为网格单元的应变状态，

ε_x，ε_y和ε_z分别为正应变沿笛卡尔坐标系的分量；

步骤S130：根据广义胡克定律，可得：

σ＝λTr(ε)I+2με (2)

其中，λ和μ为代表弹性网格材料属性的拉梅常量，

E为弹性模量，ν为泊松比，Tr(ε)＝ε_x+ε_y+ε_z，σ为网格单元的应力状态；

步骤S140：根据公式(1)和(2)得到弹性体网格变形模型：

▽·μ▽U+▽λ(▽·U)+▽·(μ▽U)^T＝0 (3)

其中，

E为弹性模量和ν为泊松比。

3.根据权利要求1所述的基于几何约束的弹性体网格变形方法，其特征在于，所述步骤S300中求解约束方程的步骤，包括以下步骤：

步骤S310：将待变形的所述弹性体网格作为有限元网格，建立有限元方程求解约束方程；

步骤S320：建立所述弹性体网格变形的边界条件，边界条件包括受力边界条件和位移边界条件；

步骤S330：求解所述网格单元变形后的线性代数方程组，得到所述具有几何约束的变形网格节点坐标。

4.根据权利要求3所述的基于几何约束的弹性体网格变形方法，其特征在于，所述求解线性代数方程组的方法为高斯-赛德尔迭代法。

5.一种基于几何约束的弹性体网格变形装置，其特征在于，包括：

应变应力模块，用于根据待变形的弹性体网格受外力变形后，所述弹性体网格内部网格单元产生的弹性应变ε和网格单元的应力状态σ，构建含网格单元的弹性模量E的待变形弹性体网格变形模型；

几何约束模块，用于立所述网格单元的尺寸约束

其中，V为二维网格单元的面积或三维网格单元的体积，建立所述网格单元的形状约束E_θ，

其中，θ_min为二维三角形网格单元的最小内角或三维四面体网格单元的最小二面角，取网格单元的几何约束E_G为E_G＝E_V+E_θ，得到约束弹性模量E'：E'＝E_G·E＝(E_V+E_θ)·E；

求解模块，用于将约束弹性模量E'代入待变形弹性体网格变形模型中得到约束方程，求解约束方程，得到具有几何约束的变形后的网格节点坐标，根据所述网格节点坐标得到变形后的所述弹性体网格。

6.根据权利要求5所述的基于几何约束的弹性体网格变形装置，其特征在于，所述应变应力模块，包括：

拓扑模块，用于解析待变形的所述弹性体网格中所有网格节点的坐标和所有网格单元的拓扑关系，定义所述弹性体网格内任意网格节点(x,y,z)在外力作用下产生的位移矢量U(x,y,z)＝(u,v,w)，其中u、v、w分别为网格节点(x,y,z)在

三个方向的位移分量，所述弹性体网格内部网格单元的弹性应力张量满足▽·σ＝0；

线性运动模块，用于所述网格单元满足线性运动定律：

其中，ε为网格单元的应变状态，

ε_x，ε_y和ε_z分别为正应变沿笛卡尔坐标系的分量；

应力状态模块，用于根据广义胡克定律，可得：

σ＝λTr(ε)I+2με (2)

其中，λ和μ为代表弹性网格材料属性的拉梅常量，

E为弹性模量，ν为泊松比，Tr(ε)＝ε_x+ε_y+ε_z，σ为网格单元的应力状态；

变形模型模块，用于根据公式(1)和(2)得到弹性体网格变形模型：

▽·μ▽U+▽λ(▽·U)+▽·(μ▽U)^T＝0 (3)

其中，

E为弹性模量和ν为泊松比。

7.根据权利要求5所述的基于几何约束的弹性体网格变形装置，其特征在于，所述求解模块，包括：

有限元模块，用于将待变形的所述弹性体网格作为有限元网格，建立有限元方程求解约束方程；

边界条件模块，用于建立所述弹性体网格变形的边界条件，边界条件包括受力边界条件和位移边界条件；

代数求解模块，用于求解所述网格单元变形后的线性代数方程组，得到所述具有几何约束的变形网格节点坐标。

8.根据权利要求5所述的基于几何约束的弹性体网格变形装置，其特征在于，所述求解线性代数方程组的方法为高斯-赛德尔迭代法。

Images