弧长法点击这里:
对于一个非线性有限元模型,只有一个自由度 ,外荷载 ,内力为
切线刚度矩阵
假设某一荷载步迭代收敛时荷载因子, 。接下来以 开始。以下是改进弧长法手算过程
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') )
输出结果: