首页/文章/ 详情

非线性| 弧长法实例

7月前浏览6475

有关弧长法的内容已经有很多了,程序也有。即便这样,离实用还有很远的路要走。这里再发一个例子。《混凝土结构有限元分析第二版》第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.51.25]
iter1_0 = [0.751.5]

iter2 = [0.51.519]
iter2_0 = [0.751.231]

iter3 = [0.51.618]
iter3_0 = [0.751.023]
iter4 = [0.51.64]
iter4_0 = [0.750.929]
iter5 = [0.51.643]
iter5_0 = [0.750.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()

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

弧长法的Python实现

非线性| 弧长法算例

非线性|弧长法改进

来源:数值分析与有限元编程
非线性pythonUM
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-04-03
最近编辑:7月前
太白金星
本科 慢慢来
获赞 5粉丝 10文章 322课程 0
点赞
收藏
作者推荐

理想塑性材料的残余应力

如图1所示,圆杆 为理想塑性材料, , 作用在 点,然后撤去,求杆的残余应力。已知杆的半径为 。▲图1荷载 作用在杆 处,可能会有四种情况: 都处于弹性状态; 塑性而 还是弹性; 塑性而 还是弹性; 都进入塑性状态。按照弹性分析,得到 , 段已经屈服,而 段还处于弹性状态。实际上, 时,就开始屈服。此时 段内力为两段的应力分别为: 由于 段仍然处于弹性,伸长量为 屈服应变为: ▲图2▲图3当 作用在 点时, 段的应力应变行为由 移动到 , 段的应力应变行为由O移动到 。如图3所示,当撤去 ,相当于反方向施加一个大小相等的 。此时有完全的弹性变形发生, 段分别有反方向的 ,或者 此时,两部分的残余应力分别为: 移除外部荷载载将导致支反座力对弹性恢复作出响应。由于这些力会限制构件完全恢复,因此会在构件中产生残余应力。为了解决这类问题,可以将构件加载和卸载的整个循环视为正荷载(加载)与负荷载(卸载)的叠加。从O到C的荷载导致塑性应力分布,而沿CD的卸载仅导致弹性应力分布。叠加需要抵消这些荷载;然而,应力分布不会取消,因此残余应力将保留在构件中。▲图4来源:数值分析与有限元编程

未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈