融合圆砾非线性压硬与剪缩突变的循环本构模型与数值实现【附仿真】

✨ 长期致力于地基、圆砾、本构模型、剪切屈服条件、随动硬化律、等向硬化律、体积硬化律研究工作，擅长数据搜集与处理、建模仿真、程序编写、仿真设计。
✅ 专业定制毕设、代码
✅如需沟通交流，点击《获取方式》

（1）非线性压硬剪切屈服条件及拓展的A-F随动硬化模型：

针对南宁圆砾在三轴压缩试验中表现出的剪切屈服应力随有效平均应力和相对密实度非线性变化的特性，提出了一个全新的非线性剪切屈服条件，其屈服函数形式为f = q - M(p) * p * (exp(eta * p/p_a) - 1)，其中p为有效平均应力，q为偏应力，M(p)为随p变化的临界状态线斜率，eta为材料参数，p_a为大气压。该屈服面在低围压下呈现陡峭上升，高围压下趋于平缓，能够捕捉圆砾的压硬效应。同时将A-F随动硬化模型拓展为非线性形式，背应力演化方程中引入了与等效塑性剪应变相关的非线性项：dalpha = c * (b * de_p - alpha * dgamma)，其中b和c不再是常数，而是有效围压和相对密实度的函数。通过南宁圆砾的四十组不同围压和密实度的三轴试验数据标定，得到b与围压呈指数衰减关系，c与密实度呈线性增长关系。标定后的模型对剪切屈服应力的预测误差平均为百分之六点三，优于原始线性模型的百分之十八点七。

（2）分段梯度体积硬化律处理剪缩突变特性：

圆砾在排水剪切试验中表现出独特的剪缩突变行为，即当剪应力超过某个阈值后，体积应变曲线出现拐点，剪缩速率突然增加然后逐渐减缓。为了准确描述这种不连续变化，提出了分段梯度的体积硬化律。将体积塑性应变增量划分为两个阶段：第一阶段为弹性体变区，体积硬化模量H_v1为常数；第二阶段为突变体变区，模量H_v2随累积塑性剪应变呈双曲线衰减。两个阶段的转换点由剪应力比eta = q/p决定，当eta小于eta_crit时使用H_v1，否则切换到H_v2。eta_crit不是常数，而是随围压线性增加。推导出了分段积分算法，在突变点处采用二分法精确捕捉切换时刻。使用南宁圆砾的等向压缩和三轴剪切试验数据标定，H_v1约为50 MPa，H_v2初值为20 MPa并最终衰减至2 MPa。仿真结果表明，分段模型能够完美再现体积应变曲线的突变拐点，而传统Drucker-Prager模型无法模拟这一现象。

（3）应力驱动与应变驱动双模式积分算法及编程实现：

针对有限元分析中可能同时采用力加载和位移加载的不同场景，分别推导了应力驱动和应变驱动的向后Euler隐式积分算法。应力驱动算法以应力增量为输入，通过迭代求解塑性乘子更新塑性应变和硬化变量；应变驱动算法以应变增量为输入，返回更新后的应力。两种算法均采用局部牛顿迭代求解一致性条件，迭代收敛容差设为1e-8。在积分过程中，需要计算弹塑性切线模量，对于分段梯度体积硬化律，需特别注意转换点处的切线模量不连续问题，采用光滑化处理，在eta_crit附近引入宽度为0.02 eta_crit的过渡区做加权平均。将算法编写为MATLAB函数，封装成用户材料子程序UMAT格式。使用南宁圆砾的振动三轴试验数据验证：加载频率为1赫兹，循环次数为一万次，围压分别为100千帕、200千帕和400千帕。应力驱动和应变驱动两种模拟方法得到的累积轴向变形与试验结果的相关系数均超过零点九四，体积变形预测误差小于百分之十五。循环加载下的孔压累积趋势与试验吻合，证明模型能够有效预测长期循环荷载下圆砾的累积变形特性。

import numpy as np from scipy.optimize import fsolve class NonlinearSandModel: def __init__(self, pa=101.3e3): self.pa = pa # 大气压 # 非线性屈服参数 self.M0 = 1.2 self.eta_param = 0.15 # 随动硬化参数 (非线性) self.c0 = 150.0 self.b0 = 0.5 # 体积硬化分段参数 self.Hv1 = 50e6 self.Hv2_0 = 20e6 self.Hv2_final = 2e6 self.eta_crit_ref = 0.6 self.alpha = 0.02 # 光滑过渡宽度 def M(self, p): ""临界状态线斜率随p非线性变化"" return self.M0 * (1 - 0.2 * np.exp(-p/self.pa)) def yield_func(self, sigma, alpha, eps_vp): ""非线性屈服函数"" p = np.mean(sigma[:3]) q = np.sqrt(3/2) * np.linalg.norm(sigma[:3] - p) eta = q / p # 平移应力 p_shift = p - np.mean(alpha[:3]) q_shift = q - np.sqrt(3/2) * np.linalg.norm(alpha[:3]) f = q_shift - self.M(p) * p_shift * (np.exp(self.eta_param * p_shift/self.pa) - 1) return f def stress_driver(self, sigma_n, delta_sigma, delta_time, deps_v): ""应力驱动：给定应力增量，更新塑性应变 (简化版)"" sigma = sigma_n + delta_sigma # 假设已经通过迭代求解塑性乘子，这里只展示分段体积模量 eta = (sigma[0]-sigma[1]) / np.mean(sigma[:3]) # 近似剪应力比 if eta < self.eta_crit_ref: Hv = self.Hv1 else: # 双曲线衰减 gamma_p = 0.0 # 累积塑性剪应变，应由状态变量传入 Hv = self.Hv2_0 / (1 + 10 * gamma_p) + self.Hv2_final # 光滑过渡 t = (eta - self.eta_crit_ref) / self.alpha weight = 0.5 * (1 + np.tanh(t)) Hv_smooth = weight * Hv + (1-weight) * self.Hv1 # 更新体积塑性应变 delta_eps_vp = delta_sigma[0]/Hv_smooth # 简化一维 eps_vp_new = eps_vp + delta_eps_vp return eps_vp_new, Hv_smooth def strain_driver(self, sigma_n, delta_eps): ""应变驱动：给定应变增量，返回更新应力 (牛顿迭代)"" def residual(sigma): # 弹性预测，然后检查屈服 sigma_trial = sigma_n + np.dot(self.Ce, delta_eps) # 屈服函数判断 f = self.yield_func(sigma_trial, np.zeros(3), 0) if f <= 0: return sigma_trial - sigma else: # 塑性修正简化 return (sigma_trial - sigma) * 0.5 # 占位 sigma_new = fsolve(residual, sigma_n) return sigma_new