一种基于有限元方法的气球膨胀动画模拟方法.pdf

本发明公开了一种基于有限元方法的气球膨胀动画模拟方法：设置气体压强、气球密度、气球初始厚度，并从obj模型中读入数据作为气球未变形的形状，包括结点位置以及三角形信息；对模型进行预计算，包括计算各个结点的质量、气球的初始体积等；基于有限元方法求出各结点弹力、气体压力，并使用牛顿-拉夫逊迭代方法求出各结点在下一时刻的位置及速度；应用保体积算法修正求出的新位置及速度，使得气球的体积与初始体积保持一致；使用POV-Ray软件渲染模型，作为下一帧的图像；将新的结果作为当前位置，进行下一帧的计算，直到气球稳定下来不再变形。本发明简单高效，能通用在保体积的薄壳模型的变形中，并可实现气球的较真实的膨胀效果及不同厚度气球在变形中的不同效果。

1.  一种基于有限元方法的气球膨胀动画模拟方法，其特征在于步骤为：
(1)设置气球内部的气体压强、气球密度、气球初始厚度参数，并从obj模型文件中读入网格数据作为气球未变形的形状，网格数据包含一系列三角形，每个三角形的数据包括三角形的三个结点位置以及三角形信息；
(2)对步骤(1)中读入的模型进行预计算，包括计算各个结点的质量、气球的初始体积、气体的初始体积；
(3)基于有限元方法，结合预计算的结果计算各基本单元的弹力、气体压力，从而求出各结点弹力、气体压力，并使用牛顿-拉夫逊迭代方法求出各结点在下一时刻的位置及速度；
(4)应用保体积算法修正步骤(3)求出的各结点位置及速度作为各结点在下一帧的位置及速度，使得气球的体积与初始体积保持一致；
(5)使用POV-Ray软件对各结点在下一帧的位置进行渲染，作为下一帧的图像；
(6)将步骤(4)中求出的各结点的位置和速度作为当前位置和速度，回到步骤(3)进行下一帧的各结点的位置的计算，直到气球稳定下来不再变形。

2.  根据权利要求1所述的基于有限元方法的气球膨胀动画模拟方法，其特征在于：所述步骤(1)中，三角形的面积及分布尽量均匀；从obj文件中读入的模型作为气球未变形的形状，适当对obj模型进行缩放，使得后续计算精确度更高。

3.  根据权利要求1所述的基于有限元方法的气球膨胀动画模拟方法，其特征在于：步骤(3)基本单元为三棱柱单元，三棱柱单元是在其他方法中常用的三角形单元的基础上增加厚度信息，即每个三棱柱单元包含3个结点1个厚度标量；每个三棱柱单元的初始厚度为气球的初始厚度，在之后的计算中，有限元方法更新3个结点的位置，保体积算法中更新结点的位置以及每个三棱柱单元的厚度，以此来保证气球的体积在变形过程中不变。

4.  根据权利要求1所述的基于有限元方法的气球膨胀动画模拟方法，其特征在于：所述步骤(2)中，各个结点的质量根据各三棱柱单元体积及密度以及质量分布权重来计算；求出各三棱柱单元的质量后，根据权重离散到3个结点上；所述步骤(2)中，气球的初始体积由各三棱柱单元体积积分得到；所述步骤(2)中，气体的初始体积由三棱柱单元底面三角形与原点形成的四面体积分得到。

5.  根据权利要求1所述的基于有限元方法的气球膨胀动画模拟方法，其特征在于：所述步骤(2)中，预计算中还包括计算未变形状态时变形梯度中的矩阵，提前计算并存储，方便后续使用。

6.  根据权利要求1所述的基于有限元方法的气球膨胀动画模拟方法，其特征在于：所述步骤(3)中，基于有限元方法求出各气球结点弹力、气体压力，并使用牛顿-拉夫逊迭代方法求出各气球结点在下一时刻的位置及速度流程为：
(31)在每个三棱柱单元中，计算结点受到的弹力和气体压力，叠加到结点的相应弹力向量和气体压力向量中；
(32)利用步骤(31)计算出的结点弹力向量和压力向量构造系统线性方程；
(33)采用共轭梯度法求解线性方程，求解过程中，根据需要求出弹力和压力的变分，作为刚度矩阵及系数矩阵与向量的乘积；
(34)根据步骤(33)求出的结果更新结点位置及速度，回到步骤(31)进入下一次迭代；迭代次数为3-5次，根据模型而定；最后一次迭代的结果即为权利要求1中步骤(3)的最终结果。

7.  根据权利要求6所述的基于有限元方法的气球膨胀动画模拟方法，其特征在于：所述步骤(32)中的弹力向量的计算根据当前结点位置以及未变形时结点位置得到，其计算流程为：对每个三棱柱，计算变形梯度、应变张量、应力张量、变形梯度对3个结点的导数，然后计算每个三棱柱单元中3个结点受到的弹力并加到结点的弹力向量中；最终，结点的弹力向量即为该结点受到的总的弹力。

8.  根据权利要求6所述的基于有限元方法的气球膨胀动画模拟方法，其特征在于：所述步骤(32)中气体压力向量的计算根据当前气体体积、初始气体体积、初始气体压强、三棱柱单元底面面积及法向量得到，计算流程为：首先计算当前气体的体积，然后根据当前气体体积及初始气体体积和压强来得到当前气体压强，接着求出每个三棱柱单元受到的气体压力，最后将三棱柱单元气体压力离散到3个结点上并加到结点的压力向量中；最终，结点的压力向量即为该结点受到的总的气体压力。

9.  根据权利要求1所述的基于有限元方法的气球膨胀动画模拟方法，其特征在于：所述步骤(33)中弹力的变分的计算方法为：求出变形梯度的变分、应变张量的变分、应力张量的变分，然后计算每个三棱柱单元中3个结点的弹力变分并加到结点的弹力变分向量中。

10.  根据权利要求1所述的基于有限元方法的气球膨胀动画模拟方法，其特征在于：
所述步骤(4)中应用保体积算法修正步骤(3)求出的各结点位置及速度流程为：
(41)将每个三棱柱单元的厚度信息转化为虚拟结点，作为三棱柱单元的第四个结点，并将该结点加入系统的结点列表中；
(42)在每个三棱柱单元中，计算增加了第四个结点的共4个结点的新质量；
(43)计算保体积的目标函数，即当前气球体积与初始体积的差，然后计算目标函数对每个结点包括虚拟结点的导数；
(44)根据上述结果计算每个结点的位置修正量，包括虚拟结点；
(45)根据修正量，更新包括虚拟结点在内的所有结点的位置及速度，然后在每个三棱柱单元中，根据4个结点的新位置计算新的厚度。

一种基于有限元方法的气球膨胀动画模拟方法
技术领域
本发明属于计算机图形学领域，具体涉及一种基于有限元方法的气球膨胀动画模拟方法。
背景技术
随着计算机的高速发展，基于物理的模拟大量用在游戏、动画、影视特效等方面，比如刚体的运动、水流等流体的运动、肌肉皮肤等变形体的变形，布料、头发等的动态模拟。在变形体的模拟中，薄壳模型的模拟是一大难点，气球即为典型的薄壳模型，且是形变较大的薄壳模型。
薄壳模型与薄板模型非常相似，他们都是厚度都远远小于其他两个维度的物体。薄壳模型与薄板模型的区别在于未变形状态的形状。薄壳模型在未变形时是曲面，比如树叶、帽子、气球等。而薄板模型在未变形时是平面，比如布料。薄壳模型与薄板模型的厚度维度远远小于其他维度，所以无法将其当作空间中的一般三维物体来模拟，否则容易出现病态问题，大大降低模拟的精确度。由于薄壳模型与薄板模型的相似性，很多人使用同样的方法对两种模型进行模拟，忽略厚度，将模型近似为空间中的平面或曲面来计算。但是对于薄壳模型，特别是非常不规则的薄壳模型来说，很难计算其二维参数坐标，因此上面的近似方法在实现时有较大的困难。另一方面，忽略厚度对于某些物体来说失掉了很重要的信息，比如本方法中的气球模型。不同的厚度，气球各处产生的拉伸不同，显示效果也会不同，忽略厚度信息，明显会造成绘制时的不真实。
另外，气球与其他薄壳模型相比，具有独特的特性：(1)气球在变形中起主要作用的是拉伸产生的力，弯曲产生的力相对而言是非常小的。(2)气球在变形中，横切面上基本没有明显的错切，即在厚度方向，没有明显的错切。(3)气球在变形过程中，保持体积不变。这里的体积指的是气球本身的体积，而不包括气球内气体的体积。随着气球的膨胀，气球的厚度变薄，反之亦然。这些独特的特性进一步增加了气球变形模拟的难度，目前的方法无法实现这些特性。
发明内容
本发明的目的在于提供一种基于有限元方法的气球膨胀动画模拟方法，应用到气球 膨胀模拟动画中，模拟效果真实、高效，可以实现不同厚度下的不同变形效果，是现有技术无法实现的。
为实现上述目的，本发明采用如下的技术方案：
一种基于有限元方法的气球膨胀动画模拟方法，其步骤为：
(1)设置气球内部的气体压强、气球密度、气球初始厚度参数，并从obj模型文件中读入网格数据作为气球未变形的形状，网格数据包含一系列三角形，每个三角形的数据包括三角形的三个结点位置以及三角形信息；
(2)对步骤(1)中读入的模型进行预计算，包括计算各个结点的质量、气球的初始体积、气体的初始体积；
(3)基于有限元方法，结合预计算的结果计算各基本单元的弹力、气体压力，从而求出各结点弹力、气体压力，并使用牛顿-拉夫逊迭代方法求出各结点在下一时刻的位置及速度；
(4)应用保体积算法修正步骤(3)求出的各结点位置及速度作为各结点在下一帧的位置及速度，使得气球的体积与初始体积保持一致；
(5)使用POV-Ray软件对各结点在下一帧的位置进行渲染，作为下一帧的图像；
(6)将步骤(4)中求出的各结点的位置和速度作为当前位置和速度，回到步骤(3)进行下一帧的各结点的位置的计算，直到气球稳定下来不再变形。
所述步骤(1)中，三角形的面积及分布尽量均匀；从obj文件中读入的模型作为气球未变形的形状，适当对obj模型进行缩放，使得后续计算精确度更高。
所述步骤(3)基本单元为三棱柱单元，三棱柱单元是在其他方法中常用的三角形单元的基础上增加厚度信息，即每个三棱柱单元包含3个结点1个厚度标量；每个三棱柱单元的初始厚度为气球的初始厚度，在之后的计算中，有限元方法更新3个结点的位置，保体积算法中更新结点的位置以及每个三棱柱单元的厚度，以此来保证气球的体积在变形过程中不变。
所述步骤(2)中，各个结点的质量根据各三棱柱单元体积及密度以及质量分布权重来计算；求出各三棱柱单元的质量后，根据权重离散到3个结点上。
所述步骤(2)中，气球的初始体积由各三棱柱单元体积积分得到。
所述步骤(2)中，气体的初始体积由三棱柱单元底面三角形与原点形成的四面体积分得到。
所述步骤(2)中，预计算中还包括计算未变形状态时变形梯度中的矩阵，提前计算并存储，方便后续使用。
所述步骤(3)中，基于有限元方法求出各气球结点弹力、气体压力，并使用牛顿-拉夫逊迭代方法求出各气球结点在下一时刻的位置及速度流程为：
(31)在每个三棱柱单元中，计算结点受到的弹力和气体压力，叠加到结点的相应弹力向量和气体压力向量中；
(32)利用步骤(31)计算出的结点弹力向量和压力向量构造系统线性方程；
(33)采用共轭梯度法求解线性方程，求解过程中，根据需要求出弹力和压力的变分，作为刚度矩阵及系数矩阵与向量的乘积；
(34)根据步骤(33)求出的结果更新结点位置及速度，回到步骤(31)进入下一次迭代；迭代次数为3-5次，根据模型而定；最后一次迭代的结果即为权利要求1中步骤(3)的最终结果。
所述步骤(32)中的弹力向量的计算根据当前结点位置以及未变形时结点位置得到，其计算流程为：对每个三棱柱，计算变形梯度、应变张量、应力张量、变形梯度对3个结点的导数，然后计算每个三棱柱单元中3个结点受到的弹力并加到结点的弹力向量中；最终，结点的弹力向量即为该结点受到的总的弹力。
所述步骤(32)中气体压力向量的计算根据当前气体体积、初始气体体积、初始气体压强、三棱柱单元底面面积及法向量得到，计算流程为：首先计算当前气体的体积，然后根据当前气体体积及初始气体体积和压强来得到当前气体压强，接着求出每个三棱柱单元受到的气体压力，最后将三棱柱单元气体压力离散到3个结点上并加到结点的压力向量中；最终，结点的压力向量即为该结点受到的总的气体压力。
所述步骤(33)中弹力的变分的计算方法为：求出变形梯度的变分、应变张量的变分、应力张量的变分，然后计算每个三棱柱单元中3个结点的弹力变分并加到结点的弹力变分向量中。
所述步骤(4)中应用保体积算法修正步骤(3)求出的各结点(谁的)位置及速度流程为：
(41)将每个三棱柱单元的厚度信息转化为虚拟结点，作为三棱柱单元的第四个结点，并将该结点加入系统的结点列表中；
(42)在每个三棱柱单元中，计算增加了第四个结点的共4个结点的新质量；
(43)计算保体积的目标函数，即当前气球体积与初始体积的差，然后计算目标函数对每个结点包括虚拟结点的导数；
(44)根据上述结果计算每个结点的位置修正量，包括虚拟结点；
(45)根据修正量，更新包括虚拟结点在内的所有结点的位置及速度，然后在每个 三棱柱单元中，根据4个结点的新位置计算新的厚度。
本发明与现有技术相比具有的优点和积极效果如下：
(1)本发明基于有限元方法，提出了简单的薄壳模型的变形算法，与现有技术相比简单高效，而又不会丢失厚度信息。
(2)本发明在基于有限元的变形算法基础上增加了保体积算法，实现保体积的薄壳模型模拟算法，现有技术没有实现真正的保体积效果。
(3)本发明将算法应用到气球膨胀模拟动画中，模拟效果真实、高效，可以实现不同厚度下的不同变形效果，是现有技术无法实现的。
附图说明
图1基于有限元方法的气球膨胀动画模拟方法流程图；
图2有限元方法求解流程图；
图3保体积算法流程图。
具体实施方式
为了使本技术领域的一般技术人员更好的理解本发明，以下结合附图进一步详细描述本发明所提供的一种基于有限元方法的气球膨胀动画模拟方法，但不构成对本发明的限制。
本发明以简单三棱柱单元作为有限元方法的基本单元，将厚度信息引入，模拟气球这种橡胶薄壳模型的变形，实现厚度在气球变形中的重要作用，并用简单的均匀气体模型来模拟气体，实现吹气球的动态过程，模拟不同厚度的气球在充气时的不同变形过程。该方法速度快，基本实现气球的特性，气球的体积在变形过程中基本保持不变，变形过程基本符合真实气球的变形效果，且不同厚度的气球变形不同，是其他方法所不能达到的。
图1是基于有限元方法的气球膨胀动画模拟方法流程图。方法的主要步骤如下：
(1)设置参数，如气体压强、气球密度等。
(2)从obj模型中读取数据，要求模型是三角形网格，尽量均匀。
(3)预计算。包括计算结点质量、气球初始体积等。
(4)使用有限元方法求解下一帧各个结点的位置及速度。
(5)使用保体积算法修正上一步求出的结点位置及速度，使气球体积与初始体积相等。
(6)使用POV-Ray软件渲染模型，作为下一帧的图像。
(7)回到步骤(4)，将计算出的结果作为当前位置，继续计算下一帧的图像，直 到气球稳定下来不再变形。
步骤(4)是本发明的一大重要部分，有限元方法求解部分。该部分是算法中耗时较长的一部分，该部分的速度决定了整体算法的运算速度，求出的结点的初步位置决定了本发明的模拟效果。具体求解过程如图2所示，方法如下：
(1)在每个单元中，计算结点受到的弹力和气体压力，叠加到结点的相应向量中。
(2)使用上一步计算出的结点弹力向量和压力向量构造系统线性方程。
(3)使用共轭梯度法求解线性方程，求解过程中，根据需要求出弹力和压力的变分，作为刚度矩阵及系数矩阵与向量的乘积。
(4)根据上一步求出的结果更新结点位置及速度，回到步骤(1)进入下一次迭代。迭代次数一般为3-5次，根据模型而定。最后一次迭代的结果即为图1流程图中步骤(4)的最终结果。
图1流程图中步骤(5)保体积算法是本发明另一大主要部分，该部分在有限元方法的基础上进一步修正结点的位置，使其符合气球体积不变的特性。该部分决定了本发明是否能模拟出不同厚度的效果。具体求解过程如图3所示，方法如下：
(1)将每个单元的厚度信息转化为虚拟的结点，并将该节点加入系统的结点列表中。
(2)在每个单元中，计算4个结点的新质量。
(3)计算保体积的目标函数，即当前气球体积与初始体积的差。然后计算目标函数对每个结点的导数，包括虚拟结点。
(4)根据上述结果计算每个结点的位置修正量，包括虚拟结点。
(5)根据修正量，更新包括虚拟结点在内的所有结点的位置及速度。然后在每个单元中，根据4个结点的新位置计算新的厚度。
下面具体详细说明，其步骤为：
1)模型导入、参数设置及预计算
要求模型使用三角形网格，模型数据文件为.obj格式，三角形分布尽量均匀。读入.obj格式的文件，读取结点的信息及三角形的信息，设置结点初始位置即点i在0时刻的位置也就是文件中点i的坐标，设置初始速度为0，设置参数初始气体压强p₀气球密度ρ和初始气球厚度h。考虑到气球以及薄壳模型的特性，每个单元都是由三角形的三个结点与厚度属性组成的三棱柱，引入了厚度信息的同时，而又不过大增加方法的复杂度。根据单元的体积及密度计算结点质量：
mi=Σe∈Tiwiρ·Ve]]>
其中Τⁱ表示包含结点i的单元集合，V^e是三棱柱单元的体积，是底面积与厚度的乘积。w_i是单元e将质量加载到结点i上的权重，这里使用1/3，后续的计算中会使用不同的权重计算新的质量。
计算初始气体体积，气体的体积即为气球内部的体积，可由三角形网格求出，每个三角形与坐标系的原点可以组成一个四面体，将所有四面体的体积相加，即为气球内部的体积也就是气体的体积：
V0=Σe13·13(x‾1e+x‾2e+x‾3e)·ne·Ae]]>
其中，A^e是单元底面三角形的体积，n^e是单元底面三角形的法向量，为三棱柱单元的三个结点。三个结点的中心在法向量方向的投影即为四面体的高，与底面三角形的面积的乘积的三分之一即为四面体的体积。
为方便后续计算，本步骤对每个单元作如下预计算：
Me=e‾12e‾13d‾h‾e-1]]>
d‾=e‾12×e‾13|e‾12×e‾13|]]>
其中单元中底面三角形的3个结点为厚度为边为单元底面三角形的单位法向量。
2)计算结点受到的弹力
从该步骤开始进入迭代，每次迭代都根据当前结点的速度和位置来计算下个时刻每个结点的速度和位置，即动画中下一帧的状态。每次迭代开始时，使用上次迭代中得到的下一时刻的结点速度和位置作为当前速度和位置。如果是第一次迭代，则使用步骤1)中的结点初始位置和初始速度作为当前位置和速度。
计算每个单元的变形梯度。设单元e中3个结点的当前位置为单元厚度为h^e∈R，边d为单元底面三角形的单位法向量，则通过下式计算单元e的变形梯度：
F^e＝[e₁₂ e₁₃ dh^e]·M^e
d=e12×e13|e12×e13|]]>
其中，M^e是步骤1)中的预计算结果，代表未变形状态单元的信息。变形梯度代表 物体的形变程度，与当前结点的位置和未变形位置有关，当前位置与未变形位置差距越大，则说明物体的形变越大。该变形梯度的形式从三棱柱单元推导而来，根据三棱柱中任意一点在变形前和变形后的坐标得到。不同的单元有不同的变形梯度表达形式，该形式的繁简对后续的计算有重要作用。
根据变形梯度求出单元的应变张量，应变张量代表物体的变形，更能反映物体的特性，本方法使用格林张量：
Ee=12(FeTFe-I)]]>
根据应变张量计算每个单元的1st Piola-kirchhoff应力张量。物体发生形变后，物体内部会产生应力，正如弹簧拉伸后会产生抵抗拉伸的弹力一样。应力张量即为应力的张量形式。不同的材料，应力张量的形式不同，与材料的特性以及材料的参数有关。本方法通用于不同的材料，但由于气球的橡胶材料的大形变的特殊性，本方法主要针对非线性材料。St.Venant-Kirchhoff材料的应变张量由下式计算：
P^e(F^e)＝F^e[2μE^e+λtr(E^e)I^e]
其中，μ,λ是材料的拉梅系数。
最后，根据应力张量计算出各个结点受到的弹力：
fielastic=Σefie=Σe-VeP(Fe)∂Fe∂xi]]>
在这一步中，首先根据之前求出的单元e的变形梯度对三个结点的导数，然后根据上式求出该单元中三个结点受到的弹力，加到3个结点的弹力向量中。对所有单元均进行上述计算，最终可以得到每个结点受到的弹力。总得来说，结点受到的弹力就是它在周围所有单元中受到的弹力的总和。
总结本步骤的计算过程：对每个单元，计算变形梯度、应变张量、应力张量、变形梯度对3个结点的导数，然后计算每个单元中3个结点受到的弹力并加到结点的弹力向量中。最终，结点的弹力向量即为该结点受到的总的弹力。
3)计算结点受到的气体压力
在气球的充气过程中，气体的流动不明显，气体与气球很快达到平衡状态，因此本方法采用简单的均匀气体模型，即气球内各个点的压强相同。
封闭的气球中，气体的压强与体积的乘积是一个常数。首先，与步骤1)一样，计算气体的当前体积：
V=Σe13·13(x1e+x2e+x3e)·ne·Ae]]>
气体的当前压强可由当前体积、初始体积、初始压强得到：
p=p0V0V]]>
气体内部各个点的压强都是相等的，因此根据压强以及单元底面三角形的面积，可以得到单元所受气体压力的大小，气体压力的方向即为底面三角形的法向量的方向。结点在该单元内受到的气体压力为单元受到的气体压力的三分之一。计算所有单元的气体压力，并将其离散到单元结点上，最终可得到结点的气体压力，即结点的气体压力为其周围所有单元受到的气体压力在结点上的叠加。
fip=Σe∈Ti13pAene=Σe∈Ti13p·12(e12e×e13e)]]>
4)构造系统线性方程并求解
如果是第一次迭代，则设初值定义如下变量：
Δx(k)=x(k+1)n+1-x(k)n+1,Δv(k)=v(k+1)n+1-v(k)n+1]]>
然后根据下面公式构造系统线性方程，使用牛顿-拉夫逊方法迭代求解，求出一系列近似值以及当k→∞，limk→∞x(k)n+1=xn+1,limk→∞v(k)n+1=vn+1.]]>该过程重复3-5次即可得到较准确的结果，根据具体模型而定，每次迭代中使用共轭梯度方法求出线性方程的解，求出解之后，回到步骤2)进入下一次迭代，迭代结束后进入步骤5)。
A△v_(k)＝b
A=M+ΔtC+Δt2(Kelastic(x(k)n+1)+Kp(x(k)n+1))]]>
b=M(vn-v(k)n+1)+Δt(felastic(x(k)n+1)+fext+fp(x(k)n+1)-Cv(k)n+1)]]>
v(k+1)n+1=Δv(k)+v(k)n+1]]>
x(k+1)n+1=Δt·Δv(k)+x(k)n+1]]>
其中，M是质量矩阵，是对角矩阵，对角线是各个结点的质量，其他元素为0。△t是时间步长，即当前时刻与下一时刻的时间差。vⁿ是时刻n的速度矩阵，C是阻尼矩阵，一般使用瑞利阻尼，C＝αM+βK^elastic，f^ext是除弹力、压力外其他的力，比如重力。K^elastic和K^p分别是弹力的刚度矩阵以及气体压力的系数矩阵，它们是力对结点坐标的导数，即：
Kelastic=-∂felastic∂x,Kp=-∂fp∂x]]>
之前的方法中，大多直接求出这两个矩阵。本方法中，为了加快运算速度以及降低运算的复杂度，并且考虑到这两个矩阵在求解过程中并不是单独出现的，而是与向量乘积的形式出现，即K^elastic·w,K^p·w的形式出现，w是任意向量，因此我们只需求出乘积的值即可。
对于弹力的刚度矩阵，对任意向量w，令δx＝-w，则求解目标为δf^elastic＝-K^elasticδx
具体的求解步骤如下：
对每个单元：
●求出变形梯度的变分δF^e＝[δe₁₂ δe₁₃ δdh^e]·M^e
●求出相应的应变张量
●求出相应的应力张量δP(F；δF)＝δF[2μE+λtr(E)I]+F[2μδE+λtr(δE)I]
●求出单元内各结点的弹力微分并叠加到结点的弹力微分向量中。
所有单元求解结束之后，结点的弹力微分向量即为所求。
气体压力的系数矩阵求解步骤如下：
●对每个结点i，求出结点的法向量
●将所有结点的法向量与对应δx_i＝-w_i计算内积并相加，得到标量
●对每个结点i，求出结点的压力微分该公式的具体推导过程这里不赘述。
以上是弹力微分以及压力微分的求解步骤。在系统线性方程的求解过程中，只要遇到K^elastic·w,K^p·w的形式就用上述步骤求解。
5)保体积算法
上一步求出结点在下个时刻的位置及速度后，可能导致气球的体积变化，因此这一步使用保体积算法修正结点的位置及速度，使得气球的体积不变。
为了简化算法，首先将每个单元的厚度量h转化为虚拟的点，作为单元的第四个结点。该结点位于三棱柱上底面的中心，其他三个结点在三棱柱的下底面。
x4e=13(x1e+x2e+x3e)+dhe]]>
根据新的结点重新计算各结点的质量，原有的三个结点质量计算权重为1/6，新增加的结点质量权重为1/2。然后计算保体积约束的目标函数，即当前气球体积与初始体积的差体积的计算依然由各三棱柱单元的体积积分得到。接着计算目标函数对各结点的导数：
∂C(x)∂xi=Σe∂Ce(x)∂xi=Σe:t1e=i12e24×e23+Σe:t2e=i12e13×e14+Σe:t3e=i12e14×e12+Σe:t4e=i12e12×e13]]>
该导数由每个三棱柱体积对结点的导数积分得到，其中是单元e中4个结点的索引。最后，根据以上结果得到每个结点位置的修正量：
Δxi=-(N+N1)wiΣjwjwiC(x)Σj||∂C(x)∂xj||2·∂C(x)∂xi]]>
其中w_i是体积约束中结点i的权重，这里使用质量的倒数，即w_i＝1/m_i'。N为原来的结点数，N₁是厚度h转换而来的结点数。最后根据以上修正量更新结点位置xin+1=xin+1+Δxi,]]>速度vin+1=(xin+1-xin)/Δt,]]>并将虚拟结点转化回厚度h。
6)渲染
求出结点位置后，使用POV-Ray软件渲染模型，显示结果，然后回到步骤2)进行下一帧的计算。
以上对本发明所述的一种基于有限元方法的气球膨胀动画模拟方法进行了详细的说明，但显然本发明的具体实现形式并不局限于此。对于本技术领域的一般技术人员来说，在不背离本发明所述方法的精神和权利要求范围的情况下对它进行的各种显而易见的改变都在本发明的保护范围之内。