ABAQUS隐式与显式求解法对比

算法特点

隐式求解法 与显式求解法 是有限元计算 中的两种最常用的求解方法。

Abaqus/Standard 是基于隐式求解来计算，主要特点如下

• 利用迭代方式求解下个增量步的未知量；
• 无条件稳定，任何大小的增量都可以让结果在一定取值范围内；
• 存在计算收敛问题（有些模型就是无法收敛）；

Abaqus/Explicit 是基于显式求解来计算，主要特点如下

• 显式求解仅考虑前一时刻的状态量，不需要迭代；
• 有条件稳定，增量必须小于某一临界值才能让结果达到稳定状态；
• 不存在计算收敛问题（大变形可收敛）；

实例解析

下面以一个弹簧系统的例子来看两种求解的运行方式

两弹簧的刚度系数为 $k$k ，初始状态为水平30°放置，弹簧原长为10，在中间节点施加载荷 $P$P ，就会产生对应的位移 $\delta$\delta 。随着位移 $\delta$\delta 的增加，系统等效刚度会慢慢增加，趋近于并联弹簧 $2k$2k ，这是一个典型的几何非线性问题。

用静力平衡方程可以求得解析解

$$x=\sqrt{(5+\delta)^2+(5\sqrt{3})^2}, sin\theta=\frac{5+\delta}{x}\\ \tag{1-1}$$ x=\sqrt{(5+\delta)^2+(5\sqrt{3})^2}, sin\theta=\frac{5+\delta}{x}\\ \tag{1-1}

$$P=2k(x-10)sin\theta\\ \tag{1-2}$$ P=2k(x-10)sin\theta\\ \tag{1-2}
$P(\delta)=2k(\delta+10)\left( 1-\frac{10}{\sqrt{\delta^2+10\delta+100}} \right)\\ \tag{1-3}$P(\delta)=2k(\delta+10)\left( 1-\frac{10}{\sqrt{\delta^2+10\delta+100}} \right)\\ \tag{1-3}

$K_{eq}(\delta)=2ksin\theta=2k\frac{\delta+5}{\sqrt{\delta^2+10\delta+100}} \\ \tag{1-4}$K_{eq}(\delta)=2ksin\theta=2k\frac{\delta+5}{\sqrt{\delta^2+10\delta+100}} \\ \tag{1-4}

由式(1-3)可以看到，已知位移求载荷是简单的，反之已知 $P$P 从解析方法求位移 $\delta$\delta 是比较困难的。式(1-4)是受力方向等效刚度的表达式。

隐式方法过程

隐式方法将外力P分解成若干个增量加载到结构上（以4个增量为例），计算初始状态状态的系统等效刚度系数 $K_{eq}^1$K_{eq}^1 ，应用胡克定律 直接求得 $P_{1}=0.25P$P_{1}=0.25P 时的位移量 $\delta_{1}$\delta_{1}

$$\delta_{1} = \delta_{0}+\frac{P_{1}}{K_{eq}^1}\\\tag{2-1}$$ \delta_{1} = \delta_{0}+\frac{P_{1}}{K_{eq}^1}\\\tag{2-1}

根据变形后的结果计算内外力平衡，由于刚才低估了系统刚度系数而得到残余力值 $R_{1}$R_{1}

$R_{1}=P_{1}-P(\delta_1)=2k(\delta_1+10)\left( 1-\frac{10}{\sqrt{\delta_1^2+10\delta_1+100}} \right)\\ \tag{2-2}$R_{1}=P_{1}-P(\delta_1)=2k(\delta_1+10)\left( 1-\frac{10}{\sqrt{\delta_1^2+10\delta_1+100}} \right)\\ \tag{2-2}

计算变形后的等效刚度系数 $K_{eq}^2$K_{eq}^2 ，可求得残余力施加后的总位移量 $\delta_{2}$\delta_{2}

$$\delta_{2}=\delta_{1}-\frac{R_{1}}{K_{eq}^2} \\\tag{2-3}$$ \delta_{2}=\delta_{1}-\frac{R_{1}}{K_{eq}^2} \\\tag{2-3}

接下来重新计算内外力平衡，反复迭代直到满足许用误差，这样就得到第一个增量步的结果。继续施加下一个增量外力，重复上述过程，直到左右外力施加完成。

通过python代码完成上述计算过程

import numpy as np

def P(D):
    return 2*(D+10)*(1-10/(D**2 + 10*D + 100)**0.5)

def K(D):
    return 2*(D+5)/(D**2 + 10*D + 100)**0.5


def analytical_solve(dmax,interval):
    dd = dmax / interval
    res_x = np.arange(0,dmax+dd,dd)
    res_y = P(res_x)
    return res_x, res_y

def implicit_solve(f,df,f0,d0,interval,itermax=50,errmax=1e-6):
    res_x = [d0]
    res_y = [f0]
    for i in range(1,interval+1):
        f1 = f0 + df
        dfin = df
        for j in range(1,itermax+1):
            k1 = K(d0)
            d1 = d0 + dfin / k1
            res_x.append(d1)
            res_y.append(f1)

            r1 = f1 - P(d1)
            dfin = r1
            err = abs( r1 / P(d1))
            if err<errmax:
                f0 = f1
                print("Interval-{:d} Interation-{:d}".format(i,j))
                break
            else:
                res_x.append(d1)
                res_y.append(P(d1))
                d0 = d1
    return np.array(res_x), np.array(res_y)


if __name__=="__main__":
    p = 1.0
    dp = 0.25
    itermax = 50
    interval = int(p/dp)
    d0,p0 = 0.0, 0.0

    x_analytical, y_analytical = analytical_solve(1, 100)
    x_implicit, y_implicit = implicit_solve(p,dp,p0,d0,interval,itermax)

    res_analytical = np.vstack((x_analytical,y_analytical)).T
    np.savetxt("analytical-100.csv", res_analytical, fmt='%.6f', delimiter=",")

    res_implicit = np.vstack((x_implicit,y_implicit)).T
    np.savetxt("implicit-{:.1f}-{:.3f}.csv".format(p,dp), res_implicit, fmt='%.6f', delimiter=",")

每迭代一次就记录下当前的结果和更新一次后的位置，追踪迭代过程。将计算结果数组保存到csv文件中，作图用excel就好了。

这是在分成2步加载的时候得到的结果，可以看到很快就会得到一个稳定平衡的状态，计算结果向解析解靠近。放大来看的话每个增量也是经历了4~5次迭代才达到收敛结果，再继续往下加载。

测试更多的分步加载，也都很快达到稳定状态（无条件稳定，求解系统可以自动调整增量大小，更快得到结果）。

显式方法过程

显式方法同样将外力 $P$P 分解成若干个增量加载到结构上（4个为例），计算初始状态状态的系统等效刚度系数 $K_{eq}^1$K_{eq}^1 ，应用胡克定律直接求得 $P_{1}=0.25P$P_{1}=0.25P 时的位移量 $\delta_{1}$\delta_{1}，见式（2-1）。接下来直接计算变形后的等效刚度系数 $K_{eq}^2$K_{eq}^2 ，

$$K_{eq}^2=K_{eq}(\delta_{1})\\\tag{3-1}$$ K_{eq}^2=K_{eq}(\delta_{1})\\\tag{3-1}

可求得下个增量的位移量 $\delta_{2}$\delta_{2}，

$$\delta_2=\delta_1 +\frac{0.25P}{ K_{eq}^2}\\\tag{3-2}$$ \delta_2=\delta_1 +\frac{0.25P}{ K_{eq}^2}\\\tag{3-2}

直到所有增量加载完毕，计算结束。每个增量都是直接求解的，无需迭代，整个过程如下

def explicit_solve(f,df,f0,d0,interval,errmax=1e-6):
    res_x = [d0]
    res_y = [f0]
    for i in range(1,interval+1):
        f1 = f0 + df
        k1 = K(d0)
        d1 = d0 + df / k1
        print("Interval-{:d}".format(i))
        res_x.append(d1)
        res_y.append(f1)
        # update increment
        d0 = d1
        f0 = f1

    return np.array(res_x), np.array(res_y)

可以看到，显式方法在增量较大时累积误差也较大，增量必须小于某一临界值才能让结果达到稳定状态，这个临界值就是稳定时间增量 ，它与系统的最高自然频率有关

$$\Delta t \leq \frac{2}{\omega_{max}},without-damping\\\tag{4-1}$$ \Delta t \leq \frac{2}{\omega_{max}},without-damping\\\tag{4-1}

$\Delta t \leq \frac{2}{\omega_{max}}(\sqrt{1+\xi^2}-\xi),with-damping\\\tag{4-2}$\Delta t \leq \frac{2}{\omega_{max}}(\sqrt{1+\xi^2}-\xi),with-damping\\\tag{4-2}

而系统自然频率 $\omega_{max}$\omega_{max} 与材料中疏密波传播速度 $c_{d}$c_{d} 有关

$$\omega_{max}=\frac{2c_{d}}{L}=\frac{2}{L}\sqrt{\frac{E}{\rho}}\\\tag{4-3}$$ \omega_{max}=\frac{2c_{d}}{L}=\frac{2}{L}\sqrt{\frac{E}{\rho}}\\\tag{4-3}

$$\Delta t \leq \frac{2}{\omega_{max}}=L\sqrt{\frac{\rho}{E}}\\\tag{4-4}$$ \Delta t \leq \frac{2}{\omega_{max}}=L\sqrt{\frac{\rho}{E}}\\\tag{4-4}

$L$L 表示单元在形变方向的长度或特征尺寸， $E$E 是材料弹性模量， $\rho$\rho 是单元密度或质量。

Abaqus 的 Explicit 求解步在设置时选择 Automatic 就会自动计算稳定时间增量，自动调整增量步的大小。

通常 Abaqus/ Explicit 将问题拆解成大量的小增量求解，强烈建议开启双精度，避免计算过程中的误差累积。另外，如果求解时的增量超过2000万步，也需要开启双精度，否则系统会强制计算终止，因为系统估计单精度下的累积误差已经超出许可范围。

根据式(4-4)，如果想要增加稳定时间增量，一方面是增加单元的特征尺寸，也就是用较大的网格单元，避免出现细长的或扁平的畸形单元；一方面可以增加单元密度，也就是Mass Scaling 技术。

Abaqus/Standard 与 Abaqus/Explicit 比较

Abaqus/Standard 时间增量大小没有限制，系统会根据平衡状态进行调整，逐渐增大增量，遇到不收敛时则会减小增量。适用于轻度非线性和连续平滑的非线性问题。

Abaqus/Explicit 时间增量大小受限，必须小于一定大小才可进行，与单元尺寸、质量等高度相关。不需要判断单步的误差，在可靠的增量下不存在收敛问题，除非单元被破坏。非常适用于大变形、不连续、高度非线性的问题，如碰撞、爆炸、断裂和破坏等。