弧長法點(diǎn)擊這里:

非線性 | 弧長法(Arc-Length Methods)

改進(jìn)弧長法點(diǎn)擊這里：

非線性|弧長法改進(jìn)

對(duì)于一個(gè)非線性有限元模型，只有一個(gè)自由度 ,外荷載 ,內(nèi)力為

切線剛度矩陣

假設(shè)某一荷載步迭代收斂時(shí)荷載因子, 。接下來以 開始。以下是改進(jìn)弧長法手算過程

Python代碼：

import math

def stiff(u):
    f_int = -4*(u-1)**2 + 4
    Kt = -8*(u-1)
    return f_int, Kt



f_ext = 1

lamd_0 = 3.84
u0 = 0.8
delta_lamd = 0.26

tol = 1e-7
conv = 2e11

max_iter = 200
niter = 0
sum_u = 0
sum_lamd = 0
print('niter      desplacement                   lamda                   conv')
while ( conv > tol and niter < max_iter ):
    niter += 1

    if niter == 1:
        f_int, Kt = stiff( u0 )
        lamd_1 = lamd_0 + delta_lamd
        f_resid = lamd_1 * f_ext - f_int

        delta_u = f_resid / Kt

        u1 = u0 + delta_u
        sum_u = sum_u + delta_u
        sum_lamd = sum_lamd + delta_lamd

        u0 = u1
        lamd_0 = lamd_1
    else:
        f_int, Kt = stiff( u0 )
        f_resid =  f_int - lamd_0 * f_ext 

        delta_u_ext = f_ext / Kt
        delta_u_resid = f_resid / Kt

        delta_lamd1 = delta_u_resid * sum_u / (delta_u_ext * sum_u + sum_lamd )

        lamd_1 = lamd_0 + delta_lamd1

        delta_u = delta_lamd1 *  delta_u_ext - delta_u_resid 

        sum_u = sum_u + delta_u
        sum_lamd = sum_lamd + delta_lamd1
        u1 = u0 + delta_u

        u0 = u1
        lamd_0 = lamd_1
        delta_lamd = delta_lamd1

    conv = math.sqrt( f_resid**2 )
    print(format(niter, '>3x'), format(u1, '>20.12f'), 
                    format(lamd_1, '>26.14f'), format(f_resid, '>28.16e') )

輸出結(jié)果：