有关弧长法的内容已经有很多了,程序也有。即便这样,离实用还有很远的路要走。这里再发一个例子。《混凝土结构有限元分析第二版》第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()
对于荷载因子的理解有些混乱,还得继续查文献。