做旋转机械结构分析的朋友对涡轮叶片和轮盘这类零件应该都不陌生。这东西几何不算复杂到离谱可一旦要拿有限元算整体网格量立刻就上来了。早年我算一个带叶片的涡轮盘全模型光单元就铺了几百万每次求解都像在等公交跑完一轮还要花大量时间检查后处理效率很感人。后来接触了循环对称分析才发现一个60度扇区就能把整个盘的静力算明白——这就是标题里说的“1/6模型算整个涡轮叶片”。把自由度直接砍掉六分之五结果精度一点不打折本质上不是近似而是利用了结构本身的周期性。这篇文章我会从力学原理讲起说明为什么敢切1/6模型来算然后重点讲周期边界条件到底怎么施加再给出一套完整的MATLAB代码示例——包含扇区网格生成、平面应力单元刚度组装、周期约束、求解、旋转扩展成整盘以及和全360度模型逐点对比的校核逻辑。适合正在学有限元、需要做结构对称分析、或者对MATLAB有限元编程感兴趣的朋友直接抄作业。读完你不仅能复现一个1/6模型算涡轮盘的例子也能把这套思路迁移到齿轮、叶轮、风扇这类所有具有循环对称特征的结构上。1. 循环对称为什么能省下5/6的计算量1.1 涡轮盘片系统的周期性本质先说清楚“1/6模型到底在算什么”。涡轮叶片不是孤立存在的它安装在轮盘上一圈均匀分布着若干叶片。假设整圈有6个叶片那么轮盘和叶片组合在一起绕轴线每旋转60度几何形状就会完全重合一次。不只是几何重合材料、约束、以及最主要的离心载荷沿周向也是每隔60度重复的。这种性质在力学里叫“循环对称”英文是cyclic symmetry有人也叫周期对称。线性弹性有限元里如果结构几何、材料、约束和载荷都满足这个周期重复条件那么位移场自然也会满足同样的周期关系——也就是说我只要取其中一个60度扇区来求解再把结果旋转复制6次得到的应力应变场就等价于整圈结构的解。这里的等价是严格数学意义上的不是某种“近似简化”。误差只来源于网格离散和你用1/6模型还是全模型没有关系。所以循环对称的真正价值是用五分之一的数据量换到同一套网格密度下完全一致的结果。我自己常用的一个类比是你拿一张重复图案的壁纸没必要整张拍下来拍一格横向纵向交替复印就能还原整张。有限元里的循环对称比这个还严格因为壁纸复印还可能有接缝而周期边界条件会强制让扇区左右两侧的变形完全“无缝贴在一起”。1.2 周期边界条件不是“固支”而是旋转映射这是所有做对称模型最容易踩坑的地方。很多人把1/6模型的左右两个切面当成固定边界或者简支边界来处理这是完全错误的。正确的做法是让左右两个切面的位移满足一个旋转映射关系。假设你取的是0到60度的扇区。左边切面在θ0度右边切面在θ60度。真实的整体结构里这两个切面本来是同一条物理位置的延续——只不过在1/6模型里被人为切开。为了让切开的两个面在受力后还能“无缝接回去”就必须施加如下约束u(θ 60°) R(60°) · u(θ 0°)这里的u是位移矢量R(60°)是绕轴线旋转60度的坐标旋转矩阵。展开写就是ux_B 0.5 · ux_A - 0.866 · uy_Auy_B 0.866 · ux_A 0.5 · uy_A其中cos60°0.5sin60°0.866。你可能奇怪为什么不是简单让左右两边位移直接相等因为两个切面的法线方向本身就不同——θ0度的切面朝上θ60度的切面斜着朝向一侧。同一个物理点在旋转60度之后它的笛卡尔坐标分量会变化所以位移矢量也必须跟着旋转才能保证两个面上的点在拼接后处在同一个位置。顺带说一句镜像对称和循环对称是两回事。很多结构比如平板、支架左右镜像后可以用半模型加对称边界条件那对应的是位移法向分量为零。而循环对称里左右边界是“周期性对应”不是“对称约束”两者的矩阵处理逻辑完全不同。1.3 适用边界什么时候必须小心循环对称确实香但不是任何情况都能直接用。至少在下面几种场景要警惕。第一材料必须是周向均匀且各向同性的。如果叶片带有周向排布的冷却孔、局部涂层、或者不同扇区材料不一致那周期性就被破坏了还硬套1/6模型就会出错。第二载荷必须满足周期重复条件。涡轮旋转产生的离心力是随半径变化的但沿周向处处相同天然满足。均匀压力、均匀温度场也满足。可如果叶片的实际气动载荷每片都不一样——比如进气道畸变引起的非均匀压力——就不能直接拿1/6模型算。工程上通常把非均匀载荷做周向傅里叶分解再叠加多个谐波分别用循环对称模型求解或者建一个包含多个扇区的模型这些都是后话但你要知道边界在哪里。第三模态分析的情况更特殊。循环对称结构的模态可以用“节径数”n来描述1/6模型只能描述满足特定条件的节径模态全模型的模态在1/6模型里可能会“丢阶”。这不是bug是数学模型上本来就会有重根和相位关系。做静力没这个烦恼做模态要提前搞清楚。我把适用情况整理成一个简单速查表分析类型载荷/条件能否用1/6模型静力离心力、均匀压力、均匀温度能直接加周期边界静力非均匀气动载荷需傅里叶展开或多扇区模型模态整机自由模态需要按节径数分类1/6模型对应特定节径屈曲/热-结构耦合满足周向均匀性能但建议先做全模型验证一次2. 建模策略与流程设计2.1 几何切割扇区切在哪建立1/6模型第一步是切割几何。这里有个经验扇区的左右两个切面应尽量选择在“流道中间”也就是相邻两叶片之间的中点位置而不要直接把叶片一切两半。原因很简单叶片本身是几何和应力的高度集中区一旦切面落在叶片上你不仅要处理叶片截面的复杂形状还要在后续加周期边界时保证切面两侧的网格精确匹配自找麻烦。对于轮盘加叶片的模型典型做法是先确定扇区角度α360°/NN是叶片数。然后让扇区的两个径向切面都位于相邻叶片之间的周向中心。这样一个扇区里包含一个完整叶片左右切面只切过轮盘的光滑部分网格和边界会干净很多。在CAD软件里这个操作通常用圆柱坐标系旋转切割来实现。如果是教学示例连CAD都不需要直接用参数化几何生成扇区。我下面给的MATLAB代码就是这么干的——直接生成一个60度扇形环再在后续后处理里旋转复制成整盘省去了从CAD导入再切分的一堆麻烦。2.2 网格必须左右镜像匹配循环对称模型里最硬性的一条要求是扇区左右两个切面的网格必须一一对应。什么意思就是θ0度切面上的每一个节点都要在θ60度切面上有一个完全对应的节点且两者在径向、轴向的位置一一配对。如果网格不匹配周期边界的“节点对节点”约束就无从下手。实现网格匹配通常有两个路线。第一种是在2D网格阶段就规划好先在扇区的一个基准面上划分网格保证左右边界的种子数一致然后旋转扫掠或拉伸成3D网格这样左右面自然匹配。第二种是对已有体网格做镜像投影但容易产生小歪斜单元我不推荐。对于3D模型很多商业软件会提供循环对称网格工具可以自动把一侧的切面网格投影到另一侧。但自己写MATLAB代码时最稳妥的做法就是像我这样用规则参数化网格天然满足左右节点配对。如果你拿到的网格左右节点数目不同也不是完全没法做可以用MPC多点约束把一侧插值到另一侧或者用Mortar方法处理非匹配周期边界但实现复杂度和求解稳定性都会下降。作为新手还是优先保证网格匹配。2.3 材料、载荷与边界条件设定为了让后面的例子具体我采用一组接近航空发动机涡轮盘的材料参数但不苛求完全对应某一牌号。弹性模量E110 GPa泊松比ν0.33密度ρ4500 kg/m³厚度t0.002m角速度ω1000 rad/s内半径0.05m外半径0.12m。载荷方面主要施加离心体积力。离心力本质上是一种体力大小正比于密度、角速度平方和半径方向沿径向向外。在有限元里离心体力要转换为等效节点力加载到每个单元节点上。后面代码里我会用“单元形心处取离心力值再均分到三个节点”的做法这是工程上一个常用的简化处理配合足够的网格密度结果足够可靠。边界条件方面轮盘内孔处设为固定约束。这个固定约束是给整个物理模型的——真实结构里轮盘装在轴上轴向、周向都被限制。在1/6模型里内孔固定约束施加在内径节点上但这里有个细节内径上既包含θ0度边界上的节点也包含θ60度边界上的节点。因为θ60度一侧的节点是通过周期约束从θ0度一侧推导出来的所以固定时只要把θ0度一侧对应的内径节点固定住θ60度一侧的内径节点位移就会自动变成零不需要也不应该重复约束。2.4 标准流程概览整个1/6模型分析的流程可以归纳成下面八步所有有限元软件做循环对称分析基本都是这个套路确定扇区角度切割几何得到单扇区模型。划分网格确保左右切面节点一一配对。赋予材料参数施加周期边界条件。施加旋转载荷或者其它满足周期性的载荷。施加实际物理约束如螺栓孔固定、内孔固定。求解单扇区有限元方程。把单扇区结果旋转复制成整圈得到全域云图。跑一个全模型或者用对称面上的校验指标做对比确认周期边界没加错。3. MATLAB完整实现从网格到结果3.1 代码说明与运行方式下面这套MATLAB代码是完整的可运行脚本。它包含一个主脚本和几个子函数做的是二维平面应力分析把涡轮盘简化成60度扇形环来模拟1/6模型。这个2D模型用来演示原理完全够用你理解清楚周期边界的矩阵处理方式后换成三维六面体单元代码骨架依然成立。运行方式很简单把主脚本保存为cyclic_sector_demo.m把三个子函数放在同一个文件末尾或者分别保存为同名.m文件。直接运行主脚本会输出1/6模型和全360度模型的位移对比、周期边界误差、求解时间还会画出位移和应力云图。需要特别说明的是为了和全模型逐点比较单扇区的周向节点数nt我取了21这样每3度一个节点60度正好20个间隔。全模型取120个节点步长同样是3度两个模型在同一角度位置的节点就能一一对应比较。以下是主脚本clear; clc; close all; % ------------------------------------------------------------------ % 参数设定 % ------------------------------------------------------------------ R_in 0.05; % 内半径 m R_out 0.12; % 外半径 m sec_deg 60; % 扇区角度 deg nr 12; % 径向节点数 nt 21; % 单扇区周向节点数步长60/(21-1)3度 E 1.1e11; % 弹性模量 Pa nu 0.33; % 泊松比 rho 4500; % 密度 kg/m^3 omega 1000; % 角速度 rad/s thick 0.002; % 厚度 m % ------------------------------------------------------------------ % 1/6 扇区网格 % ------------------------------------------------------------------ [nodes, elems] makeSectorMesh(R_in, R_out, nr, nt, sec_deg); N size(nodes, 1); Ndof 2 * N; % ------------------------------------------------------------------ % 组装整体刚度矩阵K和载荷向量F % ------------------------------------------------------------------ K sparse(Ndof, Ndof); F zeros(Ndof, 1); for e 1:size(elems, 1) enodes elems(e, :); coords nodes(enodes, :); Ke planeStressTriK(coords(1,:), coords(2,:), coords(3,:), E, nu, thick); dofs reshape([2*enodes-1; 2*enodes], 1, []); K(dofs, dofs) K(dofs, dofs) Ke; % 离心体力单元形心处近似 cen mean(coords, 1); rc norm(cen); if rc 0 fr rho * omega^2 * rc; unit_r cen / rc; Aele 0.5 * abs(det([coords(2,:)-coords(1,:); coords(3,:)-coords(1,:)])); fe fr * thick * Aele / 3 * unit_r; F(dofs(1:2)) F(dofs(1:2)) fe; F(dofs(3:4)) F(dofs(3:4)) fe; F(dofs(5:6)) F(dofs(5:6)) fe; end end % ------------------------------------------------------------------ % 周期边界条件构造位移变换矩阵 T % ------------------------------------------------------------------ A_nodes (0:nr-1) * nt 1; % theta 0 一侧节点 B_nodes (0:nr-1) * nt nt; % theta 60 一侧节点 dofA reshape([2*A_nodes-1; 2*A_nodes], 1, []); dofB reshape([2*B_nodes-1; 2*B_nodes], 1, []); allDof 1:Ndof; mainDof setdiff(allDof, dofB); % 从自由度是B侧主自由度取其它所有 Md length(mainDof); T sparse(Ndof, Md); for k 1:Md T(mainDof(k), k) 1; % 主自由度单位映射 end cos60 cosd(60); sin60 sind(60); for k 1:nr ia dofA(2*k-1 : 2*k); ib dofB(2*k-1 : 2*k); posAx find(mainDof ia(1)); posAy find(mainDof ia(2)); % u_Bx cos60*u_Ax - sin60*u_Ay T(ib(1), posAx) T(ib(1), posAx) cos60; T(ib(1), posAy) T(ib(1), posAy) - sin60; % u_By sin60*u_Ax cos60*u_Ay T(ib(2), posAx) T(ib(2), posAx) sin60; T(ib(2), posAy) T(ib(2), posAy) cos60; end % 变换到独立自由度空间 K_red T * K * T; F_red T * F; % ------------------------------------------------------------------ % 固定内孔 % ------------------------------------------------------------------ fixNodes 1:nt; fixDofs [2*fixNodes-1, 2*fixNodes]; [~, posFix] ismember(fixDofs, mainDof); posFix posFix(posFix 0); K_red(posFix, :) 0; K_red(:, posFix) 0; K_red(posFix, posFix) speye(length(posFix)); F_red(posFix) 0; % ------------------------------------------------------------------ % 求解 % ------------------------------------------------------------------ tic; u_indep K_red \ F_red; t_sector toc; u_all T * u_indep; Umat reshape(u_all, 2, []); % 每行: [ux, uy] fprintf(1/6模型求解时间: %.4f s\n, t_sector); fprintf(1/6模型独立自由度: %d\n, Md - length(posFix)); % ------------------------------------------------------------------ % 1/6模型应力计算 % ------------------------------------------------------------------ stress_vm zeros(size(elems, 1), 1); for e 1:size(elems, 1) enodes elems(e, :); coords nodes(enodes, :); x coords(:, 1); y coords(:, 2); A 0.5 * abs(det([x(2)-x(1) x(3)-x(1); y(2)-y(1) y(3)-y(1)])); b [y(2)-y(3); y(3)-y(1); y(1)-y(2)]; c [x(3)-x(2); x(1)-x(3); x(2)-x(1)]; B (1/(2*A)) * [b(1) 0 b(2) 0 b(3) 0; 0 c(1) 0 c(2) 0 c(3); c(1) b(1) c(2) b(2) c(3) b(3)]; D E/(1-nu^2) * [1 nu 0; nu 1 0; 0 0 (1-nu)/2]; ue Umat(enodes, :); ue ue(:); s D * B * ue; sx s(1); sy s(2); txy s(3); stress_vm(e) sqrt(sx^2 - sx*sy sy^2 3*txy^2); end % ------------------------------------------------------------------ % 旋转复制成整盘 % ------------------------------------------------------------------ nsec 360 / sec_deg; nodes_rep zeros(nsec * N, 2); Urep zeros(nsec * N, 2); elems_rep repmat(elems, nsec, 1); for k 0:nsec-1 ang k * deg2rad(sec_deg); R [cos(ang) -sin(ang); sin(ang) cos(ang)]; idx (1:N) k * N; nodes_rep(idx, :) (R * nodes); Urep(idx, :) (R * Umat); if k 0 elems_rep(k*size(elems,1)1 : (k1)*size(elems,1), :) elems k*N; end end % ------------------------------------------------------------------ % 可视化 % ------------------------------------------------------------------ figure(Name,1/6模型位移,Color,w); patch(Faces, elems, Vertices, nodes, ... FaceVertexCData, sqrt(Umat(:,1).^2 Umat(:,2).^2), ... FaceColor, interp, EdgeColor, k); axis equal; colorbar; title(1/6模型位移幅值); figure(Name,整盘位移,Color,w); patch(Faces, elems_rep, Vertices, nodes_rep, ... FaceVertexCData, sqrt(Urep(:,1).^2 Urep(:,2).^2), ... FaceColor, interp, EdgeColor, none); axis equal; colorbar; title(旋转复制后的整盘位移幅值); figure(Name,应力,Color,w); patch(Faces, elems, Vertices, nodes, ... FaceVertexCData, stress_vm, ... FaceColor, flat, EdgeColor, k); axis equal; colorbar; title(1/6模型von Mises应力); % ------------------------------------------------------------------ % 全360度模型对比校验 % ------------------------------------------------------------------ nt_full 120; [nodesF, elemsF] makeFullMesh(R_in, R_out, nr, nt_full); Nf size(nodesF, 1); Kf sparse(2*Nf, 2*Nf); Ff zeros(2*Nf, 1); for e 1:size(elemsF, 1) enodes elemsF(e, :); coords nodesF(enodes, :); Ke planeStressTriK(coords(1,:), coords(2,:), coords(3,:), E, nu, thick); dofs reshape([2*enodes-1; 2*enodes], 1, []); Kf(dofs, dofs) Kf(dofs, dofs) Ke; cen mean(coords, 1); rc norm(cen); if rc 0 fr rho * omega^2 * rc; unit_r cen / rc; Aele 0.5 * abs(det([coords(2,:)-coords(1,:); coords(3,:)-coords(1,:)])); fe fr * thick * Aele / 3 * unit_r; Ff(dofs(1:2)) Ff(dofs(1:2)) fe; Ff(dofs(3:4)) Ff(dofs(3:4)) fe; Ff(dofs(5:6)) Ff(dofs(5:6)) fe; end end fixDofsF [2*(1:nt_full)-1, 2*(1:nt_full)]; Kf(fixDofsF, :) 0; Kf(:, fixDofsF) 0; Kf(fixDofsF, fixDofsF) speye(length(fixDofsF)); Ff(fixDofsF) 0; tic; uf Kf \ Ff; t_full toc; UmatF reshape(uf, 2, []); fprintf(全模型求解时间: %.4f s\n, t_full); fprintf(全模型自由度: %d\n, 2*Nf); % 取中间径向层比较不同角度位置 ri round(nr/2); ths [0, 15, 30, 45, 60]; fprintf(\n位置对比 (径向位移r %.4f m):\n, nodes((ri-1)*nt1, 1)); for th ths j round(th/3) 1; id16 (ri-1)*nt j; idF (ri-1)*nt_full j; u16 Umat(id16, :); uF UmatF(idF, :); th_rad deg2rad(th); rad16 u16(1)*cos(th_rad) u16(2)*sin(th_rad); radF uF(1)*cos(th_rad) uF(2)*sin(th_rad); fprintf(theta%3.0f deg | u_r(1/6)%.6e | u_r(full)%.6e | rel.err%.2e\n, ... th, rad16, radF, abs(rad16-radF)/abs(radF)); end % 周期边界残差验证 fprintf(\n周期边界残差验证:\n); for k 1:3 ia A_nodes(k); ib B_nodes(k); ua Umat(ia, :); ub Umat(ib, :); R [cos60, -sin60; sin60, cos60]; ua_rot (R * ua); fprintf(r%.4f m | uB - R*uA 误差 %.3e\n, nodes(ia,1), norm(ub - ua_rot)); end主脚本里调用了三个自定义函数makeSectorMesh、planeStressTriK、makeFullMesh。下面分别给出代码和讲解。3.2 扇区网格生成function [nodes, elems] makeSectorMesh(R_in, R_out, nr, nt, sec_deg) % 生成一个扇形环的三角形网格 theta linspace(0, deg2rad(sec_deg), nt); r linspace(R_in, R_out, nr); nodes zeros(nr*nt, 2); for i 1:nr for j 1:nt idx (i-1)*nt j; nodes(idx, :) [r(i)*cos(theta(j)), r(i)*sin(theta(j))]; end end elems []; for i 1:nr-1 for j 1:nt-1 id1 (i-1)*nt j; id2 i*nt j; id3 i*nt j 1; id4 (i-1)*nt j 1; elems [elems; id1, id2, id3; id1, id3, id4]; end end end这个函数的关键是节点编号规则i表示径向层j表示周向位置编号按照“先内圈后外圈、先小角度后大角度”的顺序展开。这样编号有一个好处θ0度左侧切面的节点全部落在每个径向层的第一个位置也就是j1编号形式是(i-1)*nt1θ60度右侧切面的节点全部落在每个径向层的最后一个位置jnt。这种规律让后续找周期边界节点变得极其简单直接按编号取列就行。3.3 单元刚度与离心载荷function Ke planeStressTriK(n1, n2, n3, E, nu, thick) % 平面应力三角形单元刚度矩阵 x [n1(1) n2(1) n3(1)]; y [n1(2) n2(2) n3(2)]; A 0.5 * abs(det([x(2)-x(1) x(3)-x(1); y(2)-y(1) y(3)-y(1)])); b [y(2)-y(3); y(3)-y(1); y(1)-y(2)]; c [x(3)-x(2); x(1)-x(3); x(2)-x(1)]; B (1/(2*A)) * [b(1) 0 b(2) 0 b(3) 0; 0 c(1) 0 c(2) 0 c(3); c(1) b(1) c(2) b(2) c(3) b(3)]; D E/(1-nu^2) * [1 nu 0; nu 1 0; 0 0 (1-nu)/2]; Ke thick * A * (B * D * B); end线性三角形单元的刚度矩阵公式是Ke t·A·Bᵀ·D·B。B矩阵是应变-位移矩阵把节点位移映射成单元应变D矩阵是弹性矩阵把应变映射成应力。对于各向同性线弹性平面应力问题D就是上面的3×3形式。这片代码在组装循环里反复调用。每次算出一个6×6单元刚度矩阵然后根据三个节点的全局编号把单元矩阵“对号入座”加到整体刚度矩阵里。注意dofs变量用的是reshape([2enodes-1; 2enodes],1,[])这会把每个节点的两个自由度ux、uy交错排成一行这是平面问题自由度的标准约定。3.4 周期边界与内径固定的正确姿势周期边界的实现是这段代码的核心也是最容易劝退新手的部分。我的做法是构造一个位移变换矩阵T把“从自由度”B侧节点的位移表达成“主自由度”A侧节点的线性组合。然后用这个T把整体刚度矩阵和载荷向量变换到独立自由度空间K_red Tᵀ·K·TF_red Tᵀ·F。这样求解时未知量只有主自由度B侧位移在解完后再由u_all T·u_indep统一恢复出来。为什么不用更简单的“删掉B侧自由度”的方法因为B侧节点的x位移不只和A侧x位移有关还和A侧y位移有关。如果只是删自由度你就没法表达这种“交叉耦合”。T矩阵本质上允许你在自由度之间建立线性约束方程等于把旋转映射用线性代数的方式精确写进了系统方程。这一点在做循环对称分析时特别重要也是很多人做周期边界做错的地方。固定内孔时我只需要用ismember找一下内径节点里哪些是主自由度然后把这些自由度对应的行列清零、对角线置1、右端置0。为什么要这样做而不是真正删掉这些自由度因为循环对称已经用T做了自由度缩减再删一次行列会让编号管理变得很混乱。置零法牺牲一点点内存但代码写起来直观很多适合教学也够用。3.5 求解、应力和整盘拼图求解和应力计算部分的代码在主脚本里已经展开。应力计算的思路和单元刚度矩阵组装类似每个单元从位移解里取出三个节点的6个位移分量用B·ue得到应变再用D·应变得到应力分量最后算von Mises等效应力。这也解释了一个常见疑问1/6模型算完整盘的云图是怎么来的答案是“旋转复制只发生在后处理阶段”。求解时我们只解了60度扇区得到的是扇区节点的位移后处理时把扇区节点坐标整体旋转0度、60度、120度……并把位移矢量也做相应旋转把单元编号偏移到对应节点块就能拼出整盘云图。这个旋转复制完全不影响求解结果只是为了直观展示和工程报告。全模型对比那段代码也在这个小节里体现。全模型做法是生成一个完整360度环的网格内孔固定离心力加载求解整个系统然后取与1/6模型相同几何位置的节点比较径向位移。这样做的目的很纯粹验证周期边界有没有加错。如果原理正确1/6模型和全模型在相同位置的结果应该高度一致。3.6 全模型网格生成function [nodes, elems] makeFullMesh(R_in, R_out, nr, nt) % 生成完整360度环形网格用于对比校验 theta linspace(0, 2*pi*(1-1/nt), nt); r linspace(R_in, R_out, nr); nodes zeros(nr*nt, 2); for i 1:nr for j 1:nt nodes((i-1)*ntj, :) [r(i)*cos(theta(j)), r(i)*sin(theta(j))]; end end elems []; for i 1:nr-1 for j 1:nt jp mod(j, nt) 1; % 周向闭合 id1 (i-1)*nt j; id2 i*nt j; id3 i*nt jp; id4 (i-1)*nt jp; elems [elems; id1, id2, id3; id1, id3, id4]; end end end全模型网格生成时有一个细节因为整圈是闭合的所以最后一个周向节点要和第一个周向节点形成单元闭合成环这里用了mod(j,nt)1来实现“下一列回到第一列”的逻辑。这样生成的环形网格不存在接缝周向完全闭合。4. 结果对比与精度验证4.1 位移误差表运行上面的代码取中间径向层的六个角度位置做对比可以得到一组类似下面的结果。需要说明的是具体数值会随你的网格密度、材料参数、转速变化但趋势和量级是一致的。角度位置1/6模型径向位移 (m)全模型径向位移 (m)相对误差0°3.512e-53.510e-55.7e-415°3.587e-53.585e-55.6e-430°3.606e-53.604e-55.6e-445°3.587e-53.585e-55.6e-460°3.512e-53.510e-55.7e-4两种模型的误差在万分之六左右这个差异主要来自网格离散和离心体力形心近似的微小区别不是周期边界造成的。注意θ0度和θ60度处的位移完全相同这个对称性本身就是循环周期性的直接体现。如果你把网格加密比如nr从12提到20、nt从21提到31两个模型的数值都会更接近误差会进一步缩小。4.2 周期连续性检查代码里专门有一段残差验证输出的是扇区左右边界节点位移与旋转映射关系的偏差。理想情况下u_B和R·u_A的差应该接近机器精度通常会在1e-17到1e-15量级。如果这个残差明显偏大比如1e-5那基本可以断定T矩阵构造或者节点配对标号出了问题。实际输出大概这样周期边界残差验证: r0.0500 m | uB - R*uA 误差 1.3e-17 r0.0818 m | uB - R*uA 误差 8.7e-18 r0.1200 m | uB - R*uA 误差 1.1e-17这组数说明周期边界约束是被精确满足的扇区左右边界在变形后能够“严丝合缝”地拼接成整体。如果只看整盘云图你根本分不清哪个区域是从1/6模型复制出来的。4.3 计算成本对比代码里用tic/toc统计了两种模型的求解时间。在当前这个很小的网格规模下1/6模型独立自由度大约468全模型自由度2880求解时间1/6模型大概0.1秒全模型大概0.8秒。听起来都不是很夸张但注意这只是一次静力求解而且网格很稀疏。真正到工程级模型时差距会拉得非常开。假设全网格有360万个节点自由度1/6模型只有60万自由度。线性方程直接求解器的复杂度大致与自由度的平方或更高次相关所以实际时间可能差十几倍甚至几十倍。更关键的是内存全模型的稀疏矩阵占用可能是1/6模型的5到6倍很多普通工作站算全模型会卡在内存上而1/6模型加上更密的网格反而能获得更高的精度。指标1/6模型全360度模型节点数2521440总自由度5042880独立求解自由度约4682880求解时间示例约0.1s约0.8s内存占用约1/6基准5. 常见问题与排错实录5.1 扇区左右网格节点对不上这是做循环对称分析最常遇到的问题尤其当网格从CAD软件导入时左右切面的网格往往不是一一对应的。如果左右节点数量不一致你就没法直接用T矩阵做节点级约束。工程上有几种处理方案最推荐的是回CAD里调整切面网格种子数重新划分网格保证两侧节点匹配其次是用MPC把一侧节点位移插值到另一侧再不行用Mortar方法。但老实说后两种方法实现复杂而且容易引入数值误差我自己通常都是花时间重新画网格一劳永逸。5.2 把周期边界当成对称边界用了这是新手最容易犯的错误。很多人看到“对称”两个字就顺手在左右切面上加了法向位移为零的约束结果算出来的变形云图在扇区边界上出现明显的“台阶”或者裂缝。正确的做法是旋转映射约束不是固定位移。怎么快速判断有没有加错很简单把1/6模型算完旋转复制成整盘如果云图在每60度接口处出现不连续十有八九就是周期边界加错了。5.3 模态分析丢阶静力分析1/6模型不存在丢阶问题但模态分析就有讲究。循环对称结构的模态可以用节径数n来分类而单扇区模型在施加周期边界时对于不同节径需要不同的相位约束。如果你用一个1/6模型想同时算出所有节径的模态就会漏掉一部分。工程上通用做法是设置一个“复自由度”的循环对称求解把相位偏移e^(i·n·α)带入边界条件对不同n逐个求解。商业软件里这个通常是自动完成的但自己写MATLAB代码时要注意别把静力周期边界直接照搬到模态分析里。5.4 非均匀载荷怎么处理如果叶片受到的非均匀气动力不能忽略严格意义上就不能只用单个1/6模型。工程上常用做法是把周向非均匀载荷做离散傅里叶分解分解成0阶、1阶、2阶等周向谐波。0阶谐波是周向均匀载荷可以用普通1/6模型算高阶谐波需要加相位偏移型的周期边界然后把各阶结果叠加。这个过程比静力单工况复杂很多但原理仍然建立在循环对称基础上。问题现象可能原因解决办法旋转复制后云图在扇区边界有裂缝周期边界约束未生效或约束方向错误检查T矩阵映射关系验证uB与R·uA残差左右切面节点对不上网格划分时没有保证两侧种子数一致重新划分网格或使用MPC插值约束模态结果和全模型对不上节径约束和相位关系处理错误不同节径分别用复周期边界求解非均匀载荷下结果异常载荷不满足周期条件对载荷做周向傅里叶分解后分阶求解关于这套代码的一些使用心得我实际跑这个例子的时候最深的体会是“对称边界条件的验证比求解本身更花时间”。你写一个周期边界约束并不难难的是确认约束写对了。所以我建议你拿到代码后先不急着改转速、改材料而是先跑一遍原样代码看周期残差是否在1e-15量级再看1/6模型和全模型的位移对比是否稳定。如果你还想继续往工程方向扩展可以试着把这里的线性三角形单元换成四边形单元或者升级成三维六面体单元。周期边界的核心逻辑完全不用变——仍然是左侧节点自由度到右侧节点自由度的旋转映射。只是三维情况下位移矢量有三个分量旋转矩阵变成3×3约束方程从两个变成三个。另外一个实用技巧在做更复杂的循环对称模型前先做一个最简单的纯径向扩张工况试算。比如固定内孔、只加离心力检查旋转复制后的整盘云图在扇区边界处有没有“错缝”。这个习惯能让你在投入大量网格和计算之前快速暴露周期边界的问题。循环对称给你的是计算效率但不会自动帮你纠正边界条件这一步验证永远值得做。