Nn?x?1j?,x3?,x2jj???nj （2.64）
?nj为克罗内克符号，通过定义质心，这项?x?jdV?0，将公式2.62和2.63带入
V2.60中得：
4Eb??bi??vic0Vn?1nn （2.65）
由式2.64和形心的性质可知c0n4?14，带入式2.65中，得：
（2.66）
E?b?n?1?vin?biV4将式2.62带入2.61中，得：
4E????vin?1In?V?NndvidtdV （2.67）
最后将式2.66和2.67带入2.59中可得：
4E??n?1?vi?fi??n?n?biV4??V?Nndvi?dV? dt? （2.68）
在四面体的静力平衡方程中，在任何虚速度下，外力虚功率等于内力虚功率：由式2.58和2.68可知：
?fi?nTin3??biV4??V?NndvidtdV （2.69）
假定在四面体内，加速度是不变的，则式2.69中最后一项可改写为：
?V?Nndvi?dv?dV??i?dt?dt?n?V?NdV （2.70）
n又因为?为常数，?x?jdV?0，c0n?14，则式2.7为：
V?V?NndvidtdV??V?dvi???4?dt?n （2.71）
定义?V4为虚拟的节点质量mn，节点质量可由式2.72确定，节点质量可以保证在平衡的过程中数值计算的稳定。式2.71改写为：
?V?Nndvi?dv?dV?m?i?dt?dt?nn （2.72）
因此式2.69改写为：
?fi?nTin3??biV4n?dv??m?i??dt?n （2.73）
至此等价系统的平衡方程已建立，如式2.73，在每个节点，总的静态平衡力??f?，包括所有四面体的贡献、节点贡献荷载?p?以及集中力为0。为了描述以上规律，我们约定如下记号：如果一变量带有
献之和。采用以上记号，节点牛顿运动定律可以写成如下形式：
Fi?l??M?l??dvi????dt??l? l?1,nn
?l? （2.74）
nn为介质的节点总数，节点质量M定义如下：
M?l?????m????l? （2.75）
不平衡力?F??l?定义如下：
Fi?l???T?biV?????i??4?????3?l??Pi?l? （2.76）
当系统达到平衡时，不平衡力渐渐趋于0。 （3） 时间导数的显示有限差分
考虑本构方程2.43和变形率和节点速度之间的关系2.51，公式2.74可以写成如下形式：
dvidt?l??1M?l?Fi?l??t,?vi?1?,vi?2?,vi?3?,?,vi?p???l?,? l?1,nn
? （2.77）
记号{}?l?表示计算节点l速度的速度子集。在FLAC3D中求解2.77一般采用对时间t的显示有限差分。在这个方法中假定节点速度在时间间隔?t内线性变化。公式2.77中左边的导数采用中心差分来代替。中心差分时，当计算位移、力时只存取一半时间步的速度。节点速度通过以下迭代公式计算：
vi?l?(t??t2)?vi?l?(t??t2)??tM?l?Fi?l??t,?vi?1?,vi?2?,vi?3?,?,vi?p???l?,?? （2.78）
同理节点坐标更新也是通过中心差分的形式得到：
xi?l?(t??t)?xi?l?(t)??tvi?l?(t??t2) （2.79）
从式2.78和2.79中可知，一阶误差项不存在，因此上述迭代方法是二阶精度的节点位移迭代公式如下：
ui?l??l?(t??t)?ui?l?(t)??tvi?l?(t??t2) （2.80）
ui(0)?0
（4） 本构方程的增量形式
在FLAC3D里面假定在时间间隔?t，速度为常量，因此本构方程2.43可以写成如下形式：
?ij?Hij??ij,?ij?t? ??* （2.81）
?ij??*为共轭应力增量，Hij为一已知函数。
假定位移在时间?t内，线性变化：
?ij?t???ij （2.82）
??ij为时间t时的应变改变量。
应力增量通过???ij计算得到，如下：
?ij???ij （2.83） ??ij?????ij为一应力修正，定义如下：
??ij???ik?kj??ik?kj??tCCC （2.84）
转动率张量由公式2.40或是2.51计算
（5） 大小应变模式
以上描述的微分方程都是描述大变形，涉及到大位移，以及位移的线性变化和转动。在FLAC3D里面称之为大变形模式。
在实际应用中，转动足够小，以至于分量?ij??ij与整体相比很小，可以忽略不计。张量???可以用?I?来代替，公式2.83中的应力修正项可以忽略。当然在小位移和位移变化中，应变率张量中涉及的空间导数，可以通过初始值估计。节点坐标也不再由2.79式更新。
在FLAC3D里面小变形假定小位移、位移梯度和转动，节点坐标不再更新，以及不考虑应力转动修正。 （6） 保持数值稳定的力学时间步
公式2.78能求解出正确解答的前提是数值计算的稳定。通过一个理想的质点-弹簧系统，可以得出一些观点。质点-弹簧系统的平衡方程如下（矩阵符号表示）：
?P???K??u???M??*?dv?? dt?? （2.85）
这里大括号表示节点矢量，?P*?表示外力，?K?表示弹簧系统的刚度矩阵，?M?表示分布质量矩阵。如果我们能解释公式2.74中的不平衡力作为外力施加到弹簧系统上，以及公式2.85中弹簧的反力，那么理想材料的近似就是可取的。
在采用有限差分来分析质点-弹簧振荡系统时，时间步不能超过临界时间步，而这临界时间步与整个系统的最小固有周期有关。因此在有限差分中必须提供时间的一个上限，以此来保证数值计算的稳定。
推导临界时间步的关系，需要一些固有周期方面的知识。然而在实际中，全局的特征值分析是不可行的，也不经常采用这种方法来得到固有周期。在FLAC3D里面采用在稳定系统中的局部扰动方法来实现一个目标。在整个数值计算过程中，整个系统都是采用一不变的单位时间步?t?1。
假定公式2.74中右边节点质量为变量，并且参与整个局部平衡过程。
图3 质点弹簧系统1
首先考虑一维自由度弹簧系统，见图3.质点在给定初始位移情况下的运动方程如下：
?kx?mdxdt22 （2.86）
k为弹簧的劲度，m节点质量。这个方程的二阶有限差分临界时间步如下：
?t?T? （2.87）
这里T为系统周期。
T?2?mk （2.88）
图4 质点弹簧系统2
现在假定无穷个质点和弹簧见图4，由于结构对称，这个组合系统的力学行为可以采用图4所示简化方法分析。图4中（a）所示系统可以简化为一刚度矩阵为4k的单一弹簧系统。从式2.87和2.88中可推出极限稳定临界公式如下：
m?k??t?2 （2.89）
假定?t?1，如果节点质量大于或是等于弹簧劲度，则系统式稳定的。在局部分析中，公式2.89的有效性，通过解释m作为节点l的节点质量贡献ml，以及k和相关的节点刚度贡献kl扩展到四面体。节点质量贡献通过无穷的准则来提供一个时间步的上限来导出。节点刚度通过局部劲度矩阵的简单对角化得到。 四面体在节点l的内力贡献为：Til3，这个力可以假定弹簧的恢复力?kijuj（式
ll2.85）。假定在时间间隔dt内，可得如下公式：
dTi3l??kijvjdtll （2.90）
将公式2.57带入2.90中可得：
d?ij3njS(l)(l)??kijvjdtll （2.91）
假定在节点l的q方向上有一单位速度，其他方向上的速度为0，则从2.91式可得：
kqqdt??d?qj3njS(l)(l) （2.92）