✨ 长期致力于地基、圆砾、本构模型、剪切屈服条件、随动硬化律、等向硬化律、体积硬化律研究工作擅长数据搜集与处理、建模仿真、程序编写、仿真设计。✅ 专业定制毕设、代码✅如需沟通交流点击《获取方式》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 MPaH_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, pa101.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
融合圆砾非线性压硬与剪缩突变的循环本构模型与数值实现【附仿真】
✨ 长期致力于地基、圆砾、本构模型、剪切屈服条件、随动硬化律、等向硬化律、体积硬化律研究工作擅长数据搜集与处理、建模仿真、程序编写、仿真设计。✅ 专业定制毕设、代码✅如需沟通交流点击《获取方式》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 MPaH_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, pa101.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
