弧长法的Python实现
弧长法点击这里:非线性 | 弧长法(Arc-Length Methods)改进弧长法点击这里:非线性|弧长法改进对于一个非线性有限元模型,只有一个自由度 ,外荷载 ,内力为 切线刚度矩阵 假设某一荷载步迭代收敛时荷载因子, 。接下来以 开始。以下是改进弧长法手算过程Python代码:import mathdef stiff(u): f_int = -4*(u-1)**2 + 4 Kt = -8*(u-1) return f_int, Ktf_ext = 1lamd_0 = 3.84u0 = 0.8delta_lamd = 0.26tol = 1e-7conv = 2e11max_iter = 200niter = 0sum_u = 0sum_lamd = 0print('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') ) 输出结果:来源:数值分析与有限元编程