非线性| 弧长法实例

有关弧长法的内容已经有很多了，程序也有。即便这样，离实用还有很远的路要走。这里再发一个例子。《混凝土结构有限元分析第二版》第245页的例题，书中是手算演示，我这里用python实现。

python代码:

import math
defKT():
return1
defFint(x):
if x > 1 :
return1 - 0.15*(x-1)
else:
return2*x - x**2
f_ext = 1# 适应书中计算，假设外荷载为1
list_lamd = [0.75]  # 荷载步只有一步
tol  = 2e-2# 收敛参数，这个值主要是适应书中 5 步迭代
conv = 2e11
max_iter = 200# 单个荷载步最大迭代步数
niter = 0
last_sum_lamd = 0.75# 上一个收敛点，以下过程跟他没有关系？
last_fe  = last_sum_lamd * f_ext # 上一个收敛点对应的外荷载   1*0.75=0.75
last_sum_x = 0.5# 上一个收敛点
print('niter          x             lamda           conv')
for delta_lamd in list_lamd:
    delta_f = f_ext * delta_lamd   # 1*0.75 = 0.75
    lamd = 1
    Kt = KT()            # 常刚度迭代，同一荷载步下，切线刚度矩阵相同
    sum_x = 0
while ( conv > tol and niter < max_iter ):
        niter += 1
if niter == 1:            
            x1 = delta_f / Kt  
            sum_x = last_sum_x + x1  
            f_int = Fint(sum_x)
            sum_fe = last_fe + lamd * delta_f  #  
            f_resid = sum_fe - f_int
else:
            delta_x1 = delta_f / Kt
            delta_x2 = f_resid / Kt
            delta_lamd = (x1 * delta_x2)/(x1 * delta_x1 + 
                lamd * delta_f**2)
            delta_lamd = -1 * delta_lamd
            lamd = lamd + delta_lamd
            sum_fe = last_fe + lamd * delta_f  #
            dx = delta_lamd * delta_x1 + delta_x2
            x2 = x1 + dx
            sum_x = last_sum_x + x2
            f_int = Fint(sum_x)
            sum_fe = last_fe + lamd * delta_f  #  
            f_resid = sum_fe - f_int
            x1 = x2
        conv = math.sqrt( f_resid**2 )
        print(format(niter, '>3d'), format(sum_x, '>15.5f'), 
                    format(lamd, '>16.6f'), 
                    format(f_resid, '>18.6e') ) 
    last_sum_lamd = sum_fe
    last_sum_x =  sum_x
print(last_sum_lamd)
print(last_sum_x)

迭代过程如图所示

代码：

import numpy as np
import matplotlib.pyplot as plt
x1 = np.linspace(0,1,10)
y1 = 2*x1-x1**2
x2 = np.linspace(1,3,10)
y2 =  1 - 0.15*(x2-1)
# Figure 并指定大小
plt.figure(num=3,figsize=(8,5))
# 绘制F图像，设置 color 为 red，线宽度是 2，
plt.plot(x1, y1, color='red',linewidth=2.0)
plt.plot(x2, y2, color='red',linewidth=2.0)
# 四次迭代的数据
iter1 = [0.5, 1.25]
iter1_0 = [0.75, 1.5]
iter2 = [0.5, 1.519]
iter2_0 = [0.75, 1.231]
iter3 = [0.5, 1.618]
iter3_0 = [0.75, 1.023]
iter4 = [0.5, 1.64]
iter4_0 = [0.75, 0.929]
iter5 = [0.5, 1.643]
iter5_0 = [0.75, 0.909]
plt.plot(iter1, iter1_0, label="M1", linewidth = 2, color = "green",marker='o', markersize=7)
plt.plot(iter2, iter2_0, label="M1", linewidth = 2, color = "green",marker='o', markersize=7)
plt.plot(iter3, iter3_0, label="M1", linewidth = 2, color = "green",marker='o', markersize=7)
plt.plot(iter4, iter4_0, label="M1", linewidth = 2, color = "green",marker='o', markersize=7)
plt.plot(iter5, iter5_0, label="M1", linewidth = 2, color = "green",marker='o', markersize=7)
# 设置 x，y 轴的范围以及 label 标注
plt.xlim(0,3)
plt.ylim(0,1.7)
plt.xlabel('x',fontsize = 18)
plt.ylabel('F',fontsize = 18)
# 显示图像
plt.show()

对于荷载因子的理解有些混乱，还得继续查文献。