手搓非线性有限元的那些事儿
钢和混凝土结构的非线性有限元求解器，MATLAB编程，含 15 个示例文件： 1-Example_Concrete_Compressive_DC.m: 采用位移控制的受压混凝土截面 2-Example_Concrete_Flexure_AC：使用弧长控制的受弯混凝土梁 3-Example_Concrete_Flexure_Brick_DC.m：使用位移控制的受弯混凝土梁（元素类型为砖） 4-Example_Concrete_Flexure_Tetra_DC.m：使用位移控制的受弯混凝土梁（元素类型为四面体） 5-Example_Concrete_Flexure_Wedge_DC.m：使用位移控制的受弯混凝土梁（元素类型为楔形） 6-Example_Concrete_Tensile_DC.m：使用位移控制受拉的混凝土截面（元素类型为楔形 7-Example_Elastic_Flexure_LC.m：使用荷载控制的受弯弹性梁 8-Example_Reinforced_Concrete_Bond_Flexure_DC.m：使用位移控制的钢筋混凝土梁受弯（假设钢筋和混凝土之间存在粘结滑移） 9-Example_Reinforced_Concrete_Rigid_Flexure_DC.m：使用位移控制的钢筋混凝土梁受弯（假定钢筋与混凝土之间存在刚性结合） 10-Example_Steel_Flexure_AC.m：使用弧长控制的受弯钢梁 11-示例_Steel_Flexure_DC.m：采用位移控制的受弯钢梁 12-Example_Steel_Flexure_LC.m：使用荷载控制的受弯钢梁 13-Example_Steel_Tensile_DC.m：使用位移控制的受拉钢梁 14-Example_Steel_Flexure_Brick_DC：使用位移控制受弯的细长钢梁（元素类型为砖） 15-Example_Steel_Flexure_Brick_TL_DC：使用位移控制受弯的细长钢梁（元素类型为非线性几何[总拉格朗日]砖）
最近在折腾钢结构和混凝土的非线性有限元求解器，MATLAB写起来真是酸爽。别被"非线性"吓到，其实拆开来看都是体力活。今天拿几个例子唠唠，代码和翻车经验都在这了。
位移控制才是真香定律
看第一个例子ExampleConcreteCompressive_DC.m，混凝土受压截面分析。位移控制的核心在迭代步长控制，这里用牛顿-拉夫逊法的变形版本：
while abs(R)/R0 > tol Kt = assembleTangentStiffness(u); % 组装切线刚度 delta_u = -Kt\R; % 求解位移增量 u = u + alpha*delta_u; % 带松弛因子的更新 [R, F_int] = computeResidual(u); % 计算残余力 end这里有个坑——混凝土压溃后的软化段，刚度矩阵可能奇异。解决办法是引入阻尼因子alpha，我通常从0.5开始试探。试过纯牛顿法直接崩，分分钟教你做人。
弧长法：走钢丝的艺术
ExampleSteelFlexure_AC.m用弧长法处理钢梁屈曲。关键在约束方程：
lambda = sqrt(delta_s^2 - delta_u'*delta_u); % 球面弧长约束 if sign(det(Kt)) ~= sign_prev % 判断平衡路径转向 lambda = -lambda; end钢梁屈曲后承载力下降，荷载控制直接失效。这里用行列式符号判断平衡路径转向，比单纯判断内积更稳定。但要注意当多个屈曲模态耦合时，这招可能翻车。
钢和混凝土结构的非线性有限元求解器，MATLAB编程，含 15 个示例文件： 1-Example_Concrete_Compressive_DC.m: 采用位移控制的受压混凝土截面 2-Example_Concrete_Flexure_AC：使用弧长控制的受弯混凝土梁 3-Example_Concrete_Flexure_Brick_DC.m：使用位移控制的受弯混凝土梁（元素类型为砖） 4-Example_Concrete_Flexure_Tetra_DC.m：使用位移控制的受弯混凝土梁（元素类型为四面体） 5-Example_Concrete_Flexure_Wedge_DC.m：使用位移控制的受弯混凝土梁（元素类型为楔形） 6-Example_Concrete_Tensile_DC.m：使用位移控制受拉的混凝土截面（元素类型为楔形 7-Example_Elastic_Flexure_LC.m：使用荷载控制的受弯弹性梁 8-Example_Reinforced_Concrete_Bond_Flexure_DC.m：使用位移控制的钢筋混凝土梁受弯（假设钢筋和混凝土之间存在粘结滑移） 9-Example_Reinforced_Concrete_Rigid_Flexure_DC.m：使用位移控制的钢筋混凝土梁受弯（假定钢筋与混凝土之间存在刚性结合） 10-Example_Steel_Flexure_AC.m：使用弧长控制的受弯钢梁 11-示例_Steel_Flexure_DC.m：采用位移控制的受弯钢梁 12-Example_Steel_Flexure_LC.m：使用荷载控制的受弯钢梁 13-Example_Steel_Tensile_DC.m：使用位移控制的受拉钢梁 14-Example_Steel_Flexure_Brick_DC：使用位移控制受弯的细长钢梁（元素类型为砖） 15-Example_Steel_Flexure_Brick_TL_DC：使用位移控制受弯的细长钢梁（元素类型为非线性几何[总拉格朗日]砖）
单元类型选型指南
对比三个混凝土梁弯曲案例：
- 砖单元（Brick）适合规则形状，但计算量大
- 四面体（Tetra）自动剖分方便，但精度差
- 楔形（Wedge）折中方案，适合梁状结构
看看楔形单元刚度矩阵组装：
for zeta = 1:numGauss [B, detJ] = computeStrainDisplacement(zeta,coords); Ke = Ke + B'*D*B * detJ * weight(zeta); % 单元刚度累加 end楔形单元积分点布置要特别注意，z方向建议用2点高斯积分，否则剪切自锁让你怀疑人生。试过全积分反而精度下降，这反直觉的现象就是传说中的剪力锁死。
材料模型里的门道
钢筋粘结滑移模型ReinforcedConcreteBondFlexureDC.m中，处理粘结滑移的关键：
slip = u_steel - u_concrete; % 钢筋与混凝土位移差 bond_stress = f_slip(slip); % 滑移本构关系这里用了三线性滑移模型，拐点处理要加滞回阻尼。对比刚性粘结模型，计算时间直接翻三倍，但能捕捉到裂缝开展时的刚度退化。建议调试时先用刚性模型，跑通再加滑移。
总拉格朗日格式的钢梁SteelFlexureBrickTLDC用格林应变：
E = 0.5*(F'*F - eye(3)); % 格林应变 S = constitutiveModel(E); % 第二类P-K应力 P = F*S; % 名义应力大变形问题必须用这个，但注意迭代步长要小。试过0.01的步长都能发散，后来改用自动步长调整才稳下来。
这些示例文件里藏着二十多个我踩过的坑，从单元锁死到迭代发散，从材料软化到接触收敛。非线性有限元就像谈恋爱——你以为按流程走就行，实际处处是惊喜（吓）。建议跑案例时准备好咖啡和耐心，毕竟，看见收敛曲线的那一刻，比游戏通关还带劲。