非线性| 弧长法实例

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

python代码:

import math
def KT():
    return 1

def Fint(x):
    if x > 1 :
        return 1 - 0.15*(x-1)
    else:
        return 2*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()

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