基于Matlab的预应力混凝土梁弯矩-曲率全过程分析与延性计算

简介预应力混凝土梁全过程分析Matlab源码包面向结构工程与土木工程专业学生、课程设计者及科研人员帮助完成预应力混凝土梁从受弯到破坏的全过程数值模拟并绘制弯矩-曲率关系图为理解非线性受力特性提供直观依据。压缩包体积小巧仅24KB含9个文件包括5个.m源码文件、3个.asv自动备份文件及1张运行结果示意图.m文件承载主程序与关键计算函数.asv可作为编辑过程留痕图片则展示典型输出效果。源码围绕混凝土受压应力计算、预应力筋与非预应力筋应力求解、截面弯矩曲率迭代分析等环节组织模块划分清楚适合对照学习或在此基础上扩展参数研究。该资源目前已有158人学习下载对于正在开展预应力混凝土相关课题、希望快速上手Matlab仿真分析的读者是一份简洁实用的参考资料。1. 预应力混凝土梁全过程分析为什么要画弯矩-曲率图把一根预应力混凝土梁从零加载到受压区混凝土压碎记录每个状态的弯矩和曲率连成一条完整的 M-φ 曲线这是全过程分析最核心的产出。它不止给出极限弯矩斜率代表刚度退化拐点标记开裂、钢绞线屈服、峰值和下降段积分后还能换算延性系数——这是 pushover 和抗震性能化设计里截面本构的输入。对画施工图的人这条曲线用来校核超筋对做分析的人它是把截面弹塑性行为送进 OpenSees、Abaqus 前的最后提炼。Matlab 做这事很顺手条带离散用数组求根用 fzero积分用 trapz几十行就是一个可复现的求解器。拿到这类 zip 里的源码先别急着解压运行重点对照三件事材料本构怎么写、预应力进平衡方程的方式、曲率在主循环里怎么推进。2. 弯矩-曲率分析的理论基础应变协调与预应力初始应变普通钢筋混凝土梁的 M-φ 常被简化成三段描述预应力梁在开裂前后多出一个“消压”过渡工程上按四个阶段看更清楚。每个阶段对应一种刚度状态和一套控制变量理解这张表比记公式重要。2.1 四个受力阶段与 M-φ 曲线的形态特征阶段截面状态M-φ 形态控制因素未开裂全截面受压或拉应力很小接近直线斜率约 E_c·I_g有效预应力、混凝土弹性模量开裂至消压拉区裂缝开展中和轴上移斜率明显变缓刚度折减开裂弯矩 M_cr 与配筋量屈服至峰值筋材应变快速增长受压区压应变升高曲线趋向平坦抵达峰值 M_max配筋率、混凝土强度下降段受压区混凝土软化并压碎承载力回落直到 0.85M_max 以下ε_cu、箍筋约束程度开裂和屈服要分开看正常使用状态验算裂缝宽度和挠度用的是未裂段与开裂段的刚度承载力、延性和塑性转动能力由屈服点和下降段决定。全过程分析的价值就是把两类问题统一到同一张图上取数据而不是分别套公式。传统弹性分析只给一个使用荷载下的应力结果无法回答“截面进入塑性后还能转多少”这类延性问题这正是 M-φ 分析补上的空缺。2.2 钢绞线初始应变怎么进应变协调ε_p ε_pe φ(d_p − c)预应力筋的总应变不是从零开始的。张拉阶段先建立了预拉应变 ε_pe σ_pe / E_p外荷载产生曲率 φ 后钢绞线重心高度又获得一个附加应变。把中和轴到梁顶的距离记为 cd_p 为钢绞线重心到梁顶的距离应变协调写成ε_p ε_pe φ·(d_p − c) 钢绞线总应变拉为正 ε_s φ·(d_s − c) 普通钢筋应变拉为正 ε_c(y) φ·(c − y) 混凝土条带应变压为正三个式子必须用同一套坐标。程序里常见的做法是混凝土本构按“压为正”写钢筋本构按“拉为正”写汇总轴力时统一符号N F_c − F_p − F_s 0。这里的负号是新手最容易写反的地方写错的表现是 fzero 报区间两端函数值同号或者曲线在原点附近出现不合理的“负弯矩”。预拉应变里用的是有效预应力而不是张拉控制应力。张拉控制应力是千斤顶读数锚具回缩、孔道摩擦、混凝土收缩徐变等损失已经消耗掉一部分全过程分析关心的是使用状态下继续加载的行为自然要从 σ_pe 起算。没有实测损失数据时σ_pe 取 0.600.75 f_pk 是常用做法0.70 对多数后张构件是合理中值。严格讲 φ0 时截面已有预应力引起的初始压缩求解这个状态等价于自动做一次“平衡截面”修正后续曲线都叠在它上面。多数 Matlab 源码直接用二分法把初始压应变一起解出来物理上等于把自重下的受压状态当作原点工程结论层面误差很小不必单独写一步。2.3 材料本构混凝土抛物线、钢筋双折线与钢绞线名义屈服混凝土受压采用“上升段抛物线 下降段直线”的两段式本构即可不需要 Mander 约束模型那样的复杂度除非研究箍筋对延性的定量贡献。上升段 σ_c f_c[2(ε/ε_0) − (ε/ε_0)²]下降到 0.85 f_c 处取 ε_cu向量化小函数function sig conc_stress(eps, fc, e0, ecu) % 混凝土受压本构: 压应变为正, 受拉一律输出 0 sig zeros(size(eps)); comp eps 0; sig(comp) fc * (2*eps(comp)/e0 - (eps(comp)/e0).^2); tail comp (eps e0); sig(tail) fc * (1 - 0.15*(eps(tail)-e0)./(ecu-e0)); ende0 取 0.002无约束混凝土 ecu 取 0.00330.0038截面有箍筋加密时 ecu 放宽到 0.006 以上下降段变缓延性系数明显改善。高强混凝土上升段更瘦高e0 调到 0.00220.0025 更贴近试验。普通钢筋用弹性段加屈服平台的双折线钢绞线没有屈服平台按 0.2% 残余应变对应的名义屈服强度 f_p0.2 画双折线强化段刚度取 0.010.02 E_p对峰值和下降段影响很小。拉区混凝土从开裂那一刻起不再参与受力条带法里直接让负应变输出零应力。受拉刚化效应主要影响裂缝宽度和挠度的精细估计对极限曲率和延性系数影响很小截面级计算一般忽略。3. 用 Matlab 实现弯矩-曲率求解器条带离散加二分迭代把矩形截面沿高度切成等厚条带每条带用中心点的应变代表整条带平均应变就是条带法。条带数太少应力突变在求和里产生明显锯齿条带数太多每次 fzero 要对上百个点算本构速度线性下降。120 条带对 400×800 梁是性价比不错的折中参数扫描时降到 60 也不会有数量级误差。3.1 参数表与截面离散120 条带怎么切参数含义示例值单位b, h梁宽、梁高400, 800mmn混凝土条带数120个fc, eps0, epscu抗压强度、峰值应变、极限应变38.4, 0.002, 0.0038MPa, —dp, ds预应力筋、普通钢筋重心到梁顶距离720, 760mmAps, As预应力筋、普通钢筋面积1387, 1884mm²fp02, fsy钢绞线名义屈服、钢筋屈服强度1720, 400MPaEp, Es两种筋材弹性模量1.95e5, 2.0e5MPasigma_pe有效预应力已扣损失0.70 × fp02MPa对应代码b 400; h 800; n 120; dy h / n; yc (0.5:1:n-0.5) * dy; % 每条带中心到梁顶的深度(mm) fc 38.4; eps0 0.002; epscu 0.0038; dp 720; Aps 1387; fp02 1720; Ep 1.95e5; ds 760; As 1884; fsy 400; Es 2.0e5; sigma_pe 0.70 * fp02; % 有效预应力(MPa) eps_pe sigma_pe / Ep; % 预拉应变 p struct(b,b,dy,dy,yc,yc,fc,fc,eps0,eps0,epscu,epscu, ... dp,dp,Aps,Aps,fp02,fp02,Ep,Ep, ... ds,ds,As,As,fsy,fsy,Es,Es,eps_pe,eps_pe);yc 用列向量后续应变、应力、力都按列向量运算比 for 循环快一个数量级参数打包进 struct 传递函数签名不会越写越长。注意 dp、ds 要扣除保护层和孔道灌浆的影响预应力波纹管是圆的实际等效重心按孔道中心取别直接用“主筋到边缘”的图纸标注值。3.2 应力-应变函数与截面力汇总符号约定是关键钢筋和钢绞线共用一个双折线函数钢绞线传强化刚度比 0.01钢筋传 0.005function sig steel_stress(eps, fy, E, hard) % 双折线本构: 拉应变为正, 输出拉应力为正 epsy fy / E; sig E * eps; k abs(eps) epsy; sig(k) sign(eps(k)) .* (fy hard*E*(abs(eps(k)) - epsy)); end截面力汇总写成独立函数输入顶部应变和曲率输出轴力残差和弯矩。fzero 只处理轴力平衡弯矩在主循环顺手取一次function [N, M] section_FM(et, phi, p) % et: 梁顶混凝土应变(压为正), phi: 曲率(1/mm) eps_c et - phi * p.yc; % 混凝土条带应变 sig_c conc_stress(eps_c, p.fc, p.eps0, p.epscu); % 应力(MPa) F_c sum(sig_c * p.dy) * p.b; % 混凝土受压合力(N) M_c sum(sig_c .* p.yc * p.dy) * p.b; % 绕梁顶取矩 eps_p p.eps_pe phi * p.dp - et; % 钢绞线总应变(拉为正) eps_s phi * p.ds - et; % 普通钢筋应变(拉为正) F_p steel_stress(eps_p, p.fp02, p.Ep, 0.01) * p.Aps; F_s steel_stress(eps_s, p.fsy, p.Es, 0.005) * p.As; N F_c - F_p - F_s; % 轴力平衡残差, 压为正 M M_c - F_p * p.dp - F_s * p.ds; % 截面弯矩 end这里未知量取梁顶应变 et 而不是中和轴高度原因有两个φ0 时截面应变为均匀分布中和轴不存在对 c 做区间搜索会直接失败用顶部应变后预压状态被自动计入且 et 随曲率单调上升fzero 搜索区间全程有效。符号规则在函数头写清楚混凝土压为正钢筋拉为正汇总时折成受压为正。绕梁顶取矩的写法对对称配筋没问题非对称配筋时绕塑性中心取矩更合理但对 M-φ 形态没有影响。3.3 主循环推进曲率fzero 反解顶部应变phi_max 3e-5; % 初估曲率上限, 量级约 epscu/(0.2h) phi_list linspace(0, phi_max, 600); M_list zeros(size(phi_list)); for k 1:numel(phi_list) et fzero((e) section_FM(e, phi_list(k), p), [0, 0.01]); [~, M_list(k)] section_FM(et, phi_list(k), p); end plot(phi_list*1e3, M_list/1e6, LineWidth, 1.2); xlabel(曲率 \phi (×10^{-3}/mm)); ylabel(弯矩 M (kN·m)); grid on;顶部应变搜索区间 [0, 0.01] 的含义下界 0 是零应变上界 1% 压应变已经超过无约束混凝土极限区间的两个端点必然一个使 N0钢筋拉力主导、一个使 N0混凝土压力主导符号相反保证 fzero 稳定。600 个曲率步是经验值下降段要够密才能准确找到 0.85M_max 位置太密则 fzero 调用次数线性上涨单截面 600 步计算量是秒级放心跑。跑完全曲线后如果还想要中和轴高度做应变校核直接用 c et / φ 换算φ0 时 c 趋于无穷对应均匀受压这是符合物理的。顺带一提现在有人拿 codex 这类工具直接执行 Matlab 任务模型抄循环和求根问题不大但符号约定、顶部应变换中和轴这类物理判断它看不出来所以拿到现成源码也建议按这个顺序逐行过一遍再跑。4. 弯矩-曲率图的参数敏感性分析与收敛错误排查曲线跑通只是第一步设计判断靠扫参。常见的扫参维度就三个有效预应力、普通钢筋面积、混凝土强度正好对应设计里最常调整的三个输入。4.1 有效预应力扫参开裂弯矩与屈服点的变化规律用 3.3 的求解器做扫参只需要改一个参数 p.eps_pe把 σ_pe 从 0.50 fp02 扫到 0.75 fp02曲线叠在一张图上hold on; for r 0.50:0.05:0.75 p.eps_pe r * fp02 / Ep; % 注意循环结束后要恢复原值 for k 1:numel(phi_list) et fzero((e) section_FM(e, phi_list(k), p), [0, 0.01]); [~, M_list(k)] section_FM(et, phi_list(k), p); end plot(phi_list*1e3, M_list/1e6, ... DisplayName, sprintf(σ_pe%.2f fp02, r)); end legend show; hold off;扫完会看到三类变化开裂弯矩随预压力基本线性抬升因为开裂条件近似是“预压应力 截面模量×混凝土抗拉强度”屈服点略向左下方移动预拉应变占掉了钢绞线一部分应变容量外荷载能贡献的应变增量变小峰值弯矩几乎不变因为破坏由受压区混凝土控制。极限曲率基本不动延性系数因此随预应力增大略有下降这是“预应力越大越脆”在截面层次上的直接表现。扫参项开裂弯矩峰值弯矩极限曲率延性系数σ_pe 提高明显上升基本不变基本不变略降受拉筋 As 增大不变上升下降明显下降fc 提高略升上升下降下降方向性结论用于设计判断具体数值因截面而异不要直接背比例。这里的表格是趋势参考自己扫参时记录同一曲率步长下的对比才有意义。4.2 配筋率和混凝土强度对延性的影响普通钢筋面积 As 是控制延性的第一杠杆。As 增大后受拉钢筋拉力大中和轴下移受压区高度变高峰值后受压区更容易压碎极限曲率下降同时屈服弯矩上升等效屈服曲率也变大两个方向一起把延性系数压下来。曲线上的直观特征是屈服平台变短、峰值变尖。从 M-φ 形状可以直接判断超筋倾向适筋梁的曲线在峰值附近有平缓过渡超筋梁几乎直线冲到峰值后急转直下。工程上抗震梁的受拉配筋率一般控制在平衡配筋率一定比例以内这条经验在截面曲线里对应的就是屈服平台长度。混凝土强度 fc 提高的效果类似受压区应力块更饱满极限压应变更早到达曲线变“脆”。高性能混凝土配高强钢绞线时延性系数最容易跌破预期设计里要靠箍筋约束把 ecu 提上来而不是靠加大截面这一点从 2.3 的本构参数可以直观看出来——ecu 变大下降段变缓延性直接受益。4.3 fzero 不收敛与曲线锯齿的三种常见原因第一种fzero 报区间端点函数值同号。N(et) 在 [0, 0.01] 上没跨越零点。条带数太少会让求和震荡先把 n 提到 80 以上预压力特别大时et 的解非常接近 0下限从 0 改成1e-6避免数值截断。第二种曲线在屈服点附近呈锯齿。多半是钢筋双折线强化段 hard 设成 0应力在 ε_y 处导数突变条带法对这种突变很敏感。给强化刚度 0.0050.02 E 后曲线立即光滑。第三种下降段中途跳变或回折。截面级 M-φ 的下降段本身有负刚度fzero 可能跳到另一支解。处理办法不是换求根算法而是限制曲率范围跑到第一次出现 dM/dφ 0 且 M 低于 0.85M_max 就截断后续延性计算只用峰值点索引之前的数据。判断下降段是否到底用phi_list(k) * et 0.01之类顶部应变上限做保险避免深度压碎段的无意义计算。提示收敛出问题先怀疑符号再怀疑区间最后才怀疑条带数。F_c − F_p − F_s 里任何一个符号写反现象都是曲线整体异常而不是局部波动这个特征能帮你快速定位。5. 用等能量法把弯矩-曲率曲线换算成截面延性系数5.1 提取极限曲率与等效屈服曲率的 Matlab 脚本弯矩-曲率数据本身不直接给设计者“延性”这个数需要先双线性化。常用做法是等能量法用“直线上升 水平段”的等效双线性曲线代替实际曲线要求两条曲线与横轴围出的面积相等等效屈服弯矩取峰值弯矩 Mu水平段延伸到极限曲率 φu。极限曲率取下降段第一次降到 0.85Mu 的位置[~, imax] max(M_list); Mu M_list(imax); idx find(M_list(imax:end) 0.85*Mu, 1, first) imax - 1; phi_u phi_list(idx); Area trapz(phi_list(1:idx), M_list(1:idx)); % 原曲线实际面积 phi_y 2 * (phi_u - Area / Mu); % 面积相等反解等效屈服曲率 mu_phi phi_u / phi_y; fprintf(Mu%.1f kN·m phi_y%.3e phi_u%.3e mu_phi%.2f\n, ... Mu/1e6, phi_y, phi_u, mu_phi);find 取 first 的前提是下降段单调若曲线尾部有小波动先对 M_list 做 5 点滑动平均再找 0.85Mu。面积用 trapz600 个曲率步下误差小于 1%。得到的 μ_φ 是截面曲率延性不是构件位移延性两者之间还要过塑性铰长度转换悬臂构件的常用换算公式是 μ_Δ ≈ 1 3(μ_φ − 1)(l_p/L)(1 − 0.5 l_p / L)l_p 取受压区高度或截面有效高度的经验值。5.2 手算屈服曲率校验程序输出跑得通但结果错是写数值程序最常见的坑输出延性系数前先做独立手算校验。适筋梁屈服时受拉筋应变约等于 ε_sy f_sy / E_s屈服曲率近似 φ_y ≈ ε_sy / (d_s − c_y)其中 c_y 用屈服状态的平衡估算c_y ≈ (f_sy·A_s f_p0.2·A_ps) / (α_1·f_c·b)α_1 取 1.0 足够。算出 φ_y 后和等能量法结果对比偏差 ±20% 以内正常偏差超过 3 倍先查 eps_pe 的符号——预拉应变加反了是最典型的“曲线合理但数值全错”来源。校验通过后把求解器封装成输入截面参数、输出 M-φ 曲线与延性系数的独立函数文件批量跑参数表时自动标记 0.85Mu 极限点审图时一眼能看出哪根梁延性不足比在 Excel 里手动拉曲线可靠得多。本文还有配套的精品资源点击获取