混凝土塑形损伤本构模型Matlab代码混凝土塑形损伤本构模型Matlab代码，参照规范为GB50010-2010（2015版）。 下图两算例为C40混凝土材料的受压应力应变曲线及损伤因子。 具体混凝土材料可根据代码中注释提示修改！
混凝土结构设计里最头疼的就是材料本构模型的选择，尤其是做非线性分析的时候。最近在搞一个剪力墙滞回分析的项目，发现直接用理想弹塑性模型根本没法模拟混凝土的压碎破坏，于是翻出了GB50010附录C里的混凝土应力-应变关系公式，自己捣鼓了个损伤本构的Matlab实现。
先上干货，核心算法部分长这样：
function [stress, damage] = concrete_constitutive(strain, fc, ft, Ec, alpha_c) % 应变符号判断（受压为负） if strain <= 0 % 受压区计算（规范公式C.2.3-1） x = -strain / (0.002 + 0.0338*(fc/10)^0.5); % 无量纲化应变 if x <= 1 stress = fc*(1.2*x - 0.2*x.^6); else stress = fc*(x/(alpha_c*(x-1).^1.7 + x)); end % 受压损伤因子（经验公式） damage = 1 - exp(-22*(x^2)); else % 受拉区计算（规范公式C.3.4） x = strain / 0.00015; stress = ft*(1.2*x - 0.2*x.^6).*(x<=1) + ft*(x/(0.31*(x-1).^1.7 + x)).*(x>1); % 受拉损伤简化计算 damage = 0.9*(1 - exp(-10*x)); end end这个函数最妙的地方在于用分段函数处理压拉状态。注意看第4行的应变归一化处理，分母里的0.002对应混凝土峰值压应变，后面的0.0338(fc/10)^0.5这个系数是规范里考虑不同强度等级混凝土的调整项。比如C40混凝土的fc=26.8MPa，算出来的分母就是0.002+0.0338(26.8/10)^0.5≈0.0035，和规范表格里的数值能对上。
损伤因子的计算是个经验活，这里用了指数函数来模拟损伤发展。受压损伤在第12行的1 - exp(-22*(x^2))，这个22的系数控制着损伤发展速度。做过标定的同学应该知道，系数越大损伤发展越陡峭。之前试过用幂函数，但发现指数函数在应变较大时更稳定。
调用的时候需要注意材料参数的输入顺序，举个栗子：
% C40混凝土参数 fc = 26.8; % 轴心抗压强度(MPa) ft = 2.39; % 轴心抗拉强度 Ec = 3.25e4; % 弹性模量 alpha_c = 2.0; % 下降段调整系数 strain = linspace(-0.0035, 0.002, 500); % 应变范围 [stress, damage] = arrayfun(@(e) concrete_constitutive(e, fc, ft, Ec, alpha_c), strain);这里有个坑要注意：arrayfun虽然写起来优雅，但实际计算大规模数据时会明显变慢。建议改成向量化运算，把代码里的所有运算符都加上点运算（比如x.^6改成x.^6），这样计算速度能提升5倍以上。
混凝土塑形损伤本构模型Matlab代码混凝土塑形损伤本构模型Matlab代码，参照规范为GB50010-2010（2015版）。 下图两算例为C40混凝土材料的受压应力应变曲线及损伤因子。 具体混凝土材料可根据代码中注释提示修改！
损伤因子曲线有个很有意思的现象——受拉损伤比受压损伤发展得快得多。从下面的典型输出可以看到，当拉应变达到0.001时损伤因子已经到0.7了，而压应变要到-0.002左右才会出现明显损伤。这其实和混凝土微裂缝发展规律一致，受拉区一旦开裂就迅速失去刚度。
!应力-应变曲线与损伤因子
（假装这里有图）
参数alphac需要根据试验数据调整。做过钢管混凝土柱的朋友应该知道，约束混凝土的下降段更平缓，这时候可以把alphac调到3.0以上。最近帮某设计院调试剪力墙模型时，他们拿出的试验数据显示下降段比规范公式更陡，把alpha_c调到1.5才吻合。
最后提醒改参数的朋友：规范里的公式都是基于棱柱体试件，如果要做梁柱节点这种存在复杂应力状态的区域，建议在损伤因子里叠加剪切损伤项。我们项目里就加了这么一段：
% 在原有损伤基础上叠加剪切损伤 tau = ... % 剪切应力计算 damage = max(damage, 0.7*(1 - exp(-5*(tau/3.5).^3))); % 3.5MPa为剪切强度当然这属于魔改范畴了，发论文的时候记得在方法里写清楚修正项。
