CAE有那些重要的算法？？大家讨论讨论

关键词：自由面；渗流有限元分析；丢单元；算法

　　多孔介质的渗流分析问题一直是水利水电工程的一项热门研究课题。如何高效简捷地处理并确定出渗流自由面及渗流逸出点位置，并保证迭代过程收敛的稳定性是该问题的关键所在。多年来，从变网格有限元数值分析到固定网格的各种算法中，反映出该问题至今仍未得到完全地解决。

1 各种算法及存在的问题

　　传统的变网格算法最大的优点就是，渗流自由面和逸出点可以随着求解渗流场的迭代过程逐步稳定而自行形成，迭代过程是收敛的。但该算法对有复杂夹层和复杂排水系统的水工结构处理起来太困难，几乎不可能实现；另外对初始渗流自由面位置的假定要求也较高，如果初始位置与最终自由面位置相距甚远，则极易造成单元严重畸变，影响计算的精度；由于网格的变形，渗流与应力耦合分析也无法在同一网格下进行，其应用受到一定的局限。

　　目前应用较广的是固定网格法，即保持有限元网格不变的前提下求解渗流场。固定网格法在克服变网格法不足的同时，又面临了新的困难。较有代表性的当属剩余流量法和传导矩阵调整法。剩余流量法通过不断求解流过自由面的法向流量(称为剩余流量)建立求解水头增量的线性代数方程组，达到修正全场水头和调整新的自由面位置的目的。迭代过程中只需一次形成总传导矩阵，但需要判断自由面被单元分割的各种情形，要求算出穿过单元的自由面被单元切割的面积及流过自由面的法向流速，计算工作量很大，难以推广到三维问题中。初流量法在此基础上作了重大改进，利用高斯点的水头求出节点的初流量作为求解水头增量的右端项，避免了求自由面被切割的面积，同时避免了每次迭代中确定自由面的位置的做法，大大简化了剩余流量法的计算工作量，但是收敛稳定性较差，而且由于两种算法的整个迭代过程依赖于第一次有限元计算的结果，精度受到一定的影响。

　　单元传导矩阵调整法利用对渗流场有限元计算的结果，根据单元8节点水头与节点位置势的比较，把渗流厨行分区，各区的渗透系数给不同的值，通过不断调整单元传导矩阵，模拟渗流不饱和区的作用，来确定出真实的渗流饱和区及渗流场。该算法实际上是把边界不确定的非线性问题转化成了材料非线性问题来考虑。但是大量算例表明，该算法对含自由面单元组成的区域的处理往往使迭代过程的收敛性和解的稳定性难尽人意，结果出现振荡的现象。针对这一问题，有了许多改进的算法，如加密高斯点法，虚单元法，复合材料单元法等。加密高斯点法旨在对含自由面单元通过加密高斯点的做法，试图更精确地模拟该类单元对总传导矩阵的贡献；复合材料单元法把渗透系数在不饱和区取实际渗透系数的千分之一，在饱和区取实际渗透系数，在介于饱和区与不饱和区的界面渗透系数出现跳跃的现象改为在一个不大的邻域(-ε，ε)内，用连接该邻域端点的渗透系数值作一线性函数过渡，避免由于含自由面单元的复合缝面处渗透系数发生突变影响收敛的稳定性；事实上用一个渗透系数的小值模拟不饱和区仍然会对真实的渗流场起作用，而且单元内任一点的水头和坐标都是节点水头和坐标的插值函数，是一个连续可微函数，而上述做法并不能真实反映这一区域的透水特性，计算精度和收敛稳定性都受到影响。因此虚单元法通过变动局部网格，调整自由面穿过单元的某些节点，使之落在该次迭代形成的自由面上，计算区域被严格地分成自由面以上区域(称为虚区)和自由面以下区域(称为实区)，只求解实区的渗流场，虚区不参与计算。这种算法的不足是显而易见的，每次迭代都要通过移动单元的某些节点形成自由面，尤其在迭代过程中自由面若出现上下反复振荡的时间，调整局部单元网格时容易使单元发生畸变，对三维复杂问题不适用；另外，由于单元网格的变形，高斯点位置也随之发生了变化，对单元传导矩阵的影响是必然存在的，易产生结果收敛不稳定的现象。

2 丢单元法

2.1 算法的基本方程 将三维计算域剖分为n个单元，取八节点等参元的插值函数，根据变分原理，求解无源汇、无流量交换边界的渗流场的有限元方程为：

而

式中：K为总传导矩阵；H为渗流场各节点水头组成的列向量；B_i为单元i的3×8阶几何矩阵；M_i为单元i的3×3阶材料矩阵，各向同性材料M_i=k_iI.

　　如果引入第一类边界条件，即已知边界水头值，则方程(1)可变成：

2.2 算法的基本思路　丢单元法不同于虚单元法，它不仅吸取了虚单元法在求解水头时虚区不参与计算的优点，还克服了虚单元法每次迭代需调整自由面节点位置容易造成单元出现畸异的缺陷，具体做法如下：由于真实的渗流域大小未知，自由面位置不确定，迭代初始先假定一个与上游水面同一高程的水平面为自由面边界，比最终的自由面位置显然要高。随着迭代过程的进行，第r步迭代时，利用对含自由面单元8节点水头与节点位置势的比较结果，可对第r-1步确定的计算域进行分区，记KJ为单元8节点的节点水头小于节点位置势的节点个数，KI0为单元8节点的节点水头等于节点位置势的节点个数，当KI=KJ+KI0=8时，说明该单元完全处于自由面之上，当KI<8，KJ>0时，说明自由面穿过该单元，当KI<8，KJ=0时，说明该单元完全处于自由面之下。保留含自由面及自由面以下的单元，将8节点完全处于自由面以上的单元抛掉，消除不饱和区对渗流场的影响，单元号、单元信息号、材料信息等与单元有关的一切信息都要进行修改。丢单元的同时，单元上的节点也随单元一起被丢掉，值得注意的是与含自由面单元相邻的节点仍保留。每次迭代还需要对含自由面单元调整单元传导矩阵，为了避免在含自由面单元的复合缝面处材料发生突变，拟采用分段线性函数过渡的办法，如图1所示，ε取值范围的含义是自由面附近的一个带状区域，在此区域内渗透系数取比实际渗透系数略小的值，如果有限元网格较均匀，则ε取一个小值对结果的收敛不会造成影响，如果网格形状不好，各边长的比例相差较大，

则ε取值不能太小。当穿过自由面的单元各节点水头与位置势差的总和在该领域内时，即

.还可以采

　　在处理逸出段时，对自由面以下的逸出段上的节点预先赋边界值h_i＝z_i，称这些节点为固定点。与自由面相交的逸出段上的一点不带任何约束，定义为活动点。如果第r步迭代时，逸出点水头刚好等于该点的z坐标，说明其上的单元刚刚被丢掉，该逸出点还是有约束的固定点，这时应将该点的固定边界值释放，使该点成为逸出段上的活动点；如果第r步迭代时，逸出点水头不等于该点的z坐标，说明该点仍是活动点，继续参与下一步迭代。直到最终收敛。结果收敛后，将含自由面的单元挑选出来，利用插值技术确定自由面的具体位置，而逸出点的z坐标则取活动点的水头值，x坐标再利用插值技术确定。

2.3 算法实例　考察一各向同性且均质介质的河间地块，长100m，高80m，无源汇，无垂向补给，上游河流稳定水位为h1＝60m，下游水位为h₂＝20m.

　　文献[1]给出了该例自由面位置的解析解与高斯点法形成的数值解。在该算例中ε取0.01，w取1，丢单元算法的计算结果与之比较于表1所示。从表1可见，丢单元法精度是较高的，其结果值得信赖。