弧长法的Python实现_非线性_python

弧长法的Python实现

弧长法点击这里:

非线性 | 弧长法(Arc-Length Methods)

改进弧长法点击这里：

非线性|弧长法改进

对于一个非线性有限元模型，只有一个自由度 ,外荷载 ,内力为

切线刚度矩阵

假设某一荷载步迭代收敛时荷载因子, 。接下来以开始。以下是改进弧长法手算过程

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') )

输出结果：

来源：数值分析与有限元编程

材料的名义应力、应变与真实应力、应变转换公式的推导

材料的名义(Nominal)应力、应变是基于变形前的数据计算得到，其中为试件初始截面面积，为试件初始长度。名义应力、应变也叫工程(Engineering)应力、应变。CAE软件需要采用基于变形后的应力、应变，即真实的应力、应变。其中为试件当前截面面积，为试件当前长度。两种应力、应变的转化公式为：下面来推导这两个公式。一) 了解定积分的精确定义。点击这里：定积分的精确定义(重排版)二) 根据试件的体积不变的原则可得，即于是三) 假设荷载分为个增量步，且每个增量步产生相同的伸长量 ,如图所示总应变再和定积分的精确定义比较故 import mathimport matplotlib.pyplot as pltimport numpy as npplt.rcParams['font.sans-serif'] = ['SimHei'] # 正常显示中文标签plt.rcParams['axes.unicode_minus'] = False # 正常显示负号# 名义应力应变eps_N = np.array([0, 0.00060117, 0.0010815,0.0017279,0.0022288,0.0028608,0.00348034, 0.00404227,0.0045566,0.0051164,0.0058212,0.00674012,0.0078904,0.0092845, 0.010934,0.012828,0.014924,0.017178])sigma_N = np.array([490.393,503.8289, 512.22608,520.62322,529.020377,537.41752,545.81467,554.21182, 562.60897,571.00612,579.40327,587.800419,596.197567,604.594716,612.99186, 621.3890,629.78616,638.18331])n = len(sigma_N)print(n)#真实应力应变sigma_T = np.zeros((n))eps_T = np.zeros((n))for i in range(n): sigma_T[i] = sigma_N[i] *( 1 + eps_N[i] ) eps_T[i] = math.log( 1 + eps_N[i] )v1 = np.array( [ 0 ] )eps_N1 = np.hstack( (v1, eps_N) )sigma_N1 = np.hstack( (v1, sigma_N) )eps_T1 = np.hstack( (v1, eps_T) )sigma_T1 = np.hstack( (v1, sigma_T) )fig, axs = plt.subplots(1, 1, figsize=(14,6) )axs.plot(eps_N1, sigma_N1, label="M1", linewidth = 2) axs.plot(eps_T, sigma_T, label="M1", linewidth = 2, color = "dimgrey")plt.legend(["Nominal","True"]) plt.xlim(-0.0001, 0.019)plt.ylim(0, 760)axs.set_xlabel('Strain', fontsize = 18)axs.set_ylabel('Stress', fontsize = 18)fig.savefig('./f118.png', dpi = 400) #保存图片 plt.show() 这里弹性应变很小，弹性段几乎成铅锤。来源：数值分析与有限元编程