★ FEATURED ARTICLE 去年做一批钢筋混凝土框架的推覆分析被 ABAQUS 里 CDP 模型Concrete Damaged Plasticity即混凝土塑性损伤本构模型的参数标定折腾得不轻。算出来的整体曲线看着像模像样可把中间某个积分点的应力拉出来一查明显不合理——问题就出在损伤变量 d 的定义和规范公式对不上。后来我干脆自己用 Matlab 手写了一套完整的混凝土塑性损伤本构模型代码单轴受拉受压、损伤演化、再到三维应力更新框架全部按 GB 50010-2010《混凝土结构设计规范》附录 C 的公式体系实现。从规范条文到可算函数的这段路踩了不少坑也积累了一些经验很值得单独写一篇。先说个小细节标题写的是塑形损伤严格讲这个模型应该叫塑性损伤plasticity coupled damage。塑性对应的是不可恢复的变形损伤对应的是刚度退化两个机制耦合在一起才构成 CDP 模型的全貌。下面所有内容里我都按塑性损伤来讲。1. 从规范条文到可计算函数这套代码到底在干什么1.1 CDP 模型解决什么问题混凝土在受压时一开始接近线弹性然后进入非线性强化段到达峰值后出现明显的软化段刚度逐渐减小。受拉时抗拉强度低峰值后迅速软化开裂。如果用线弹性模型来近似等于完全忽略了这个材料的后峰值行为——而抗震设计、pushover 分析、非线性时程分析里恰恰是后峰值段在决定构件的延性和耗能能力。于是需要一套描述路径依赖的本构模型。塑性部分解决加载-卸载走不同路径的问题损伤部分解决卸载刚度退化的问题。两者耦合能比较真实地复现混凝土在复杂应力状态下的拉压破坏、裂缝成核、刚度退化等现象。我在代码里对它的定位是不追求覆盖商业软件里的所有细枝末节而是把单轴拉压曲线 损伤演化 基于有效应力的塑性流动这三个核心机制用尽可能少的代码量实现出来并且保证和 GB 50010 的曲线高度吻合。1.2 规范公式在哪、怎么用GB 50010-2010 附录 C 给出了混凝土单轴受拉、受压的应力-应变曲线表达式并对受压、受拉的损伤演化参数给出了计算方法。这些公式最大的好处是所有参数都可以由混凝土强度等级推算非常适合作为工程弹塑性分析的标定源。标题里的GB5...结合混凝土本构的工程惯例指的就是这本规范。不过要注意规范附录的公式是基于名义应力-名义应变的全量形式。程序里要处理的是增量步内的应力更新所以必须把全量公式改写成由当前应变计算当前应力的函数再补充损伤变量的演化方程。这一步看着简单实际动手后会发现很多细节比如分段函数的交接点、峰值点的数值精度、损伤变量在应变很小时的奇异性。后面第 2 章和第 4 章会专门展开。1.3 代码模块划分概览我把整套逻辑拆成三层材料参数层只有一个 params 结构体承载所有标定参数方便统一修改。本构响应层单轴受拉、受压两个函数输入应变输出应力和损伤。增量更新层三维 CDP 应力更新函数输入上一增量步应力加应变增量输出更新后应力。这样拆的好处是单轴曲线可以单独用来验证参数合理性三维函数可以直接挂到后续自编有限元程序里。下面按这个层次逐一展开。2. 单轴本构的数学准备哪些公式必须吃透2.1 受拉应力-应变曲线规范把受拉曲线归一化后用 x 和 y 表示x ε/εty σ/ft其中 εt 是受拉峰值应变ft 是轴心抗拉强度。当 x ≤ 1上升段y 1.2x − 0.2x^6当 x 1下降段y x / (αt*(x − 1)^1.7 x)其中 αt 是受拉下降段参数反映曲线下降的快慢。αt 越大峰值后掉得越陡材料越脆。有个容易被忽略的点上升段公式在 x0 处的初始切线斜率是 1.2也就是说规范隐含的初始受拉模量比弹性模量 Ec 略高。这倒不是错误而是为了拟合试验数据中受拉段起始非线性明显的现象。编程时不必纠结这一点按公式写就行。2.2 受压应力-应变曲线受压曲线同样用 x ε/εcy σ/fc 归一化。上升段y n*x / ((n − 1) x^n)其中 n Ecεc/(Ecεc − fc)。这个 n 很关键当 x1 时不管 n 是多少分子分母相等y1即曲线必过峰值点。这个自洽性在编程时可以用来校验代码是否正确。下降段y x / (αc*(x − 1)^2 x)αc 是受压下降段参数和 αt 一样越大越脆。标定这些参数时工程上最常用的经验公式对应规范条文说明和表 C.2.2 的实用拟合是εc (700 172*sqrt(fc)) × 10⁻⁶fc 以 MPa 计εt 可取 ft/Ec 80e−6 ~ 100e−6作为经验值具体试验标定会更准αc 0.74 0.0015*fcαt 0.312*ft²这几个公式用来快速搭建参数库特别方便代码里我直接内置了。但有一点必须提醒标定本构参数时fc、ft 要用强度标准值或实测平均值不要用设计值——设计值里已经包含了安全系数会让曲线峰值偏矮算出来的构件承载力会整体偏低。2.3 损伤变量从刚度退化到能量等效损伤变量的本质是描述卸载刚度的退化无损状态弹性刚度是 Ec损伤后变成 (1−d)·Ec。规范给出的是名义应力下的总应变曲线而程序内部需要区分弹性应变和损伤/塑性引起的残余变形所以必须把全量应力-应变关系转化成损伤演化。工程上最常用的是能量等效法d 1 − sqrt( σ / (Ec * ε) )推导逻辑不复杂损伤材料的弹性余能等于无损材料在有效应力下的弹性余能由此得到有效应力与真实应力之间的平方根关系反推就得到上式。另一种应变等效法是 d 1 − σ/(Ec*ε)两者结果有差异。我用一个表对比方法损伤变量公式特点应变等效法d 1 − σ/(E·ε)实现简单但软化段损伤增长过快能量等效法d 1 − sqrt(σ/(E·ε))退化更平缓与构件卸载试验吻合度更好我的代码采用能量等效法和规范附录的损伤演化思路基本一致。不过规范的损伤演化公式在细节上还考虑了塑性应变与弹性应变的区分所以如果要死磕规范原文建议以纸质规范文本为准把公式先誊写一遍再写代码我就是这么干的。3. CDP 模型的理论骨架从单轴到三维3.1 屈服函数Lubliner 屈服面如果只是算单轴拉压曲线屈服函数用不上。但要把本构推广到二维、三维的积分点应力更新就必须面对屈服面。CDP 模型采用的屈服函数来自 Lubliner 等人1989在主应力空间是一个开口的锥形面。函数形式用有效应力表达F 1/(1−α) * ( αI1 sqrt(3J2) β*⟨σ̂max⟩ ) − (1−α)*c这里 I1 是有效应力张量的第一不变量反映静水压力部分说白了就是三个正应力的和形容被压成一团的程度J2 是偏应力第二不变量反映形状畸变的程度⟨σ̂max⟩ 是最大有效主应力的正部即 Macaulay 括号负值取零正值保留。系数 α 的工程标定范围通常在 0.08~0.12 之间β 需要通过单轴拉压强度比 σc/σt 来标定。屈服函数里包含 ⟨σ̂max⟩ 这一项是为了模拟混凝土在受拉时更容易破坏的机制——拉应力不大的时候屈服面就会明显内缩这是混凝土强压弱拉特性的数学表达。3.2 非关联流动法则与剪胀塑性流动的方向由塑性势决定。混凝土的剪胀角 ψ 一般取 30°~40°。如果采用关联流动即塑性势等于屈服函数会明显高估混凝土的剪胀效应所以 CDP 模型用的是 Drucker-Prager 型塑性势G sqrt( (eσttanψ)² q² ) − p*tanψ其中 p I1/3 是静水压力q sqrt(3J2) 是 Mises 等效应力e 是偏心率一般取 0.1。第二项里的 (eσttanψ)² 让塑性势在受拉区保持光滑避免数值奇异。在塑性力学上这叫非关联流动带来的直接后果是塑性应变增量的方向不等于屈服面法向整体刚度矩阵不再对称有限元求解时只能用非对称求解器——这一点在挂到自编程序时要留意。3.3 损伤耦合与刚度退化单轴标定里损伤 d 直接由应力应变关系算出来。三维情况下需要把拉、压两种损伤区分开因为混凝土在拉压反复加载比如地震作用下的框架梁时裂缝的开合对刚度退化的影响完全不同。基本表达式σ (1 − d)·D0 : (ε − εp)这里 D0 是初始弹性矩阵d 是综合损伤变量。工程上常用的拉压组合方式是d 1 − (1 − dc·sc)·(1 − dt·st)sc、st 是与应力状态有关的权重系数原则是压力主导时受压损伤权重更大拉力主导时受拉损伤权重更大。代码里我暂时用了最简单的组合方式即 d 1 − (1−dc)·(1−dt)。对于单调加载问题足够如果后面要模拟滞回行为再把权重系数补上。4. Matlab 代码实现从函数骨架到数值细节4.1 材料参数结构体先上第一段可以直接复制跑的代码。它的作用是根据混凝土强度等级快速生成标定参数% CDP_concrete_params.m % 输入 fc/ft 单位 MPaEc 单位 MPa function params CDP_concrete_params(fc, ft, Ec) params.fc fc; params.ft ft; params.Ec Ec; % 受压峰值应变 / 受拉峰值应变规范推荐经验公式 params.eps_c (700 172*sqrt(fc)) * 1e-6; params.eps_t ft / Ec 80e-6; % 受压/受拉下降段参数 params.alpha_c 0.74 0.0015 * fc; params.alpha_t 0.312 * ft^2; % 受压上升段形状系数 n params.n params.Ec * params.eps_c / ... (params.Ec * params.eps_c - params.fc); % 剪胀角默认取 35 度转弧度 params.psi 35 * pi / 180; end以某批 C30 试件的实测参数为例fc 32.4 MPaft 2.2 MPaEc 3.0e4 MPa。跑完上面的函数得到 eps_c ≈ 1.28e-3alpha_c ≈ 0.79n 约 2.7和手头几十组试验曲线的统计均值很接近。这个函数我建议每个项目跑之前先调用一次快速复核参数是否在合理区间。4.2 单轴受压/受拉函数单轴受压函数的完整实现输入压应变为负值输出压应力为负值function [sigma, d] CDP_uniaxial_comp(eps, params) fc params.fc; Ec params.Ec; eps_c params.eps_c; alpha_c params.alpha_c; n params.n; abs_eps abs(eps); % 压应变取绝对值便于处理 x abs_eps / eps_c; if x 1 y n*x / (n - 1 x^n); else y x / (alpha_c * (x-1)^2 x); end sigma -fc * y; if abs_eps eps_c * 0.1 % 进入损伤区后才计算 d d 1 - sqrt(sigma / (Ec * eps)); else d 0; end d max(0, min(1, d)); % 数值保护 end这里有个关键细节损伤变量计算式 d 1 − sqrt(σ/(Ec·ε)) 在压应变绝对值很小时ε 接近 0会得到奇异值。所以代码里加了一个阈值判断只有应变超过峰值应变的 10% 时才计算损伤早期段认为 d 0。这个阈值不是物理上的严格界限纯粹是数值保护。实际测试中只要不进入软化段太深这个近似对整体曲线几乎没有影响。受拉函数逻辑完全对称符号取正即可function [sigma, d] CDP_uniaxial_tens(eps, params) ft params.ft; Ec params.Ec; eps_t params.eps_t; alpha_t params.alpha_t; x eps / eps_t; if x 1 y 1.2*x - 0.2*x^6; else y x / (alpha_t*(x-1)^1.7 x); end sigma ft*y; if x 0.1 d 1 - sqrt(sigma / (Ec*eps)); else d 0; end d max(0, min(1, d)); end4.3 数组化调用与向量化优化单轴标定时要画完整拉压曲线通常要算几百个应变点。不少人用 for 循环一个个调用函数其实 Matlab 里可以把整个应变向量传进去。但上面函数里有 if x1 这种分支直接传向量会报错。我的做法是把内部逻辑改成向量写法用逻辑索引一次性完成分段赋值function [sigma, d] CDP_uniaxial_comp_vec(eps_vec, params) fc params.fc; Ec params.Ec; eps_c params.eps_c; alpha_c params.alpha_c; n params.n; abs_eps abs(eps_vec); x abs_eps / eps_c; y zeros(size(x)); rise x 1; fall ~rise; y(rise) n*x(rise) ./ (n - 1 x(rise).^n); y(fall) x(fall) ./ (alpha_c*(x(fall)-1).^2 x(fall)); sigma -fc*y; d zeros(size(x)); valid abs_eps 0.1*eps_c; d(valid) 1 - sqrt(sigma(valid) ./ (Ec*eps_vec(valid))); d max(0, min(1, d)); end这种写法画 500 个点的曲线基本是瞬时的。许多 Matlab 新手习惯用 for 循环遍历数组其实能用向量和逻辑索引解决的分支就不要写循环代码可读性也好很多。4.4 三维应力更新框架三维部分我给的是可以直接挂到有限元积分点的框架。核心是弹性预测-塑性修正两步根据当前应变增量算弹性预测应力σ_trial σ_old D0 : Δε。计算屈服函数 F(σ_trial)如果 F ≤ tol完全按弹性路径更新不产生塑性应变。如果 F tol进入塑性修正投影回屈服面同时更新塑性应变与损伤。骨架代码如下弹性刚度矩阵 D0 的构造大家都很熟省略% CDP_3D_stress_update.m 骨架 function [sigma_new, d_new, statev] CDP_3D_stress_update(... sigma_old, deps, params, statev) D0 elastic_stiffness(params); % 6x6 弹性矩阵 sigma_trial sigma_old D0 * deps; % 弹性预测 F yield_function(sigma_trial, params); % 有效应力屈服判断 if F 1e-8 % 弹性路径 sigma_new sigma_trial; d_new statev.d; statev.deps_p zeros(6,1); return; end % 塑性修正牛顿迭代求塑性乘子 lambda lambda 0; sigma sigma_trial; residual F; iter 0; while abs(residual) 1e-8 iter 20 % 计算塑性流动方向 nG dG/dsigma nG plastic_flow_dir(sigma, params); % 返回映射sigma sigma_trial - lambda*D0*nG sigma sigma_trial - lambda * (D0 * nG); residual yield_function(sigma, params); lambda lambda residual / (nG * D0 * nG); iter iter 1; end % 更新塑性应变和损伤 statev.eps_p statev.eps_p lambda * nG; [dc, dt] calculate_damage(statev, params); d_new 1 - (1-dc)*(1-dt); sigma_new (1 - d_new) * sigma; % 损伤耦合后的名义应力 end这个骨架保留了完整三维 CDP 的核心弹性预测、屈服面判断、返回映射、塑性应变更新、损伤耦合。但要注意真正的屈服面投影需要主应力空间的变换把 6 维应力向量转到主方向求 σ1、σ2、σ3还要处理 Macaulay 括号的解析梯度。这两块展开写够一篇单独的博客这里把框架交给你顺着填就是完整可用的版本。4.5 可视化验证脚本曲线验证是每次跑参数库前必做的一步。我习惯把受拉受压两个函数合在一起画一眼检查参数是否合理% plot_uniaxial_curves.m clear; close all; fc 32.4; ft 2.2; Ec 3.0e4; params CDP_concrete_params(fc, ft, Ec); eps_c linspace(0, -0.0035, 400); eps_t linspace(0, 0.0005, 400); [sig_c, d_c] CDP_uniaxial_comp_vec(eps_c, params); [sig_t, d_t] CDP_uniaxial_tens_vec(eps_t, params); figure; subplot(2,1,1); plot(eps_c*1e3, sig_c, b-, eps_t*1e3, sig_t, r-, LineWidth, 1.6); hold on; plot([-3.5 0.5], [0 0], k:); grid on; xlabel(应变 ε (×10^{-3})); ylabel(应力 σ (MPa)); legend(受压,受拉); subplot(2,1,2); plot(eps_c*1e3, d_c, b-, eps_t*1e3, d_t, r-, LineWidth, 1.6); xlabel(应变 ε (×10^{-3})); ylabel(损伤变量 d); legend(受压,受拉); grid on;这个脚本能迅速暴露参数问题峰值点是否落在 fc/ft 上、软化段是否太陡、损伤是否在 0~1 单调递增。我每次换一组材料参数第一件事就是跑这个图。5. 验证与调试三个刷经验的坑下面分享实际调试中遇到的三个问题比公式推导更容易帮后来者省时间。5.1 峰值点漂移罪魁祸首是 n 的有效数字程序第一次跑受压曲线时我发现峰值应力不在 εc 处而是偏了一点点。查了半天问题出在 n 的计算上。n Ecεc/(Ecεc − fc)如果圆整误差太大x1 时 y 不完全等于 1峰值点自然就漂移了。解决办法保留足够的有效数字并且在函数里用 x1 而不是 x1 做判断强制把 x1 归入上升段这样峰值点稳稳落在 (εc, fc) 上。5.2 损伤变量出现负值或超过 1软化段结束时应力降到很低但应力/应变比可能小于 0导致 sqrt 里出现负值。尤其受压曲线大应变时压应力还保留一定的残余强度σ/(E·ε) 接近 0开根号没问题但受拉软化到接近 0 时σ/(E·ε) 可能因为残余强度而出现负值。加 max/min 夹紧之后d 稳定在 0~1。严格说残余强度的存在意味着 d 永远小于 1这是材料和数值的良好衔接不要硬改成 1。5.3 单位不一致导致的差三个数量级有次从一份资料里抄了组参数算出来曲线几乎贴着坐标轴。排查半天发现对方给的 fc 单位是 kN/m²我直接当 MPa 用了。混凝土本构计算必须统一单位制应力用 MPa即 N/mm²、弹性模量用 MPa、应变无量纲。我习惯在文件头部写一行注释所有力学量单位统一为 MPa / mm 体系跑任何新模型前先核对一遍。6. 扩展方向和应用建议6.1 挂到自编有限元框架这套本构的最终目标是进积分点。如果自编有限元程序只需要把 CDP_3D_stress_update 挂到材料模型接口返回应力增量和一致切线刚度矩阵即可。如果要把成果放进商业软件通常需要改写成交互式 UMAT/VUMAT 格式本构核心完全一致。我自己的流程是先用单轴脚本验证参数再挂框架跑单点 test最后再进整体模型每一步的隔离都做得很干净。6.2 用试验数据反演参数规范经验公式给出的是平均水平实际试验测出来的曲线会有偏差。这时候可以用最小二乘法拟合 αc、αt。方法很简单把试验应变数据丢进来调整 αc 直到计算曲线和试验曲线的残差最小。整个过程十几行 Matlab 就能搞定比盲调参数快得多。反演时要注意先固定峰值前的参数由 Ec、fc 确定只调下降段参数否则多个参数一起动拟合结果没有物理意义。6.3 黏性正则化让软化段收敛软化段会导致有限元计算的强烈网格依赖这是 CDP 类模型的老问题。工程界的标准处理是给损伤变量加黏性正则化常用二阶黏性方程ḋ_vis (d − d_vis) / μμ 是黏性参数静力分析时通常取 1e-5 ~ 1e-3 量级它没有实际时间意义只是一个正则参数。把 d_vis 代入应力更新能明显改善软化段的收敛性。这一步是从能算走向算得稳的关键值得单独写一篇展开。最后分享一个小习惯我每次用这套代码算一个新构件之前都会先跑一遍 4.5 的可视化脚本把单轴拉压曲线打印出来贴在屏幕边上。看起来很简单但这个习惯真的帮我避免了好几次参数抄错和单位混乱的低级失误。模型再复杂本构的物理合理性永远是第一道门槛。
阅读完成 · 觉得有帮助？
