1. 从一组应变片说起：为什么主应力方向值得单独写套代码

手里捏着一组应变花读数，三个方向、三个微应变值，纸面上看就是三个数字，但真正要拿去判断结构受力状态时，这三个数字本身几乎没有任何直接意义。你需要把它们换算成主应力的大小和方向，才能回答“这个测点到底是被拉还是被压、沿哪个方向最危险”这类工程上真正关心的问题。我接触过的不少现场测试报告里，应变花数据算错方向角的情况相当常见，问题往往不出在公式本身，而是出在角度基准、符号约定和象限判断这几个环节上，而这些恰恰是教科书上一笔带过、实操中却最容易翻车的地方。

这篇内容围绕一个很具体的事情展开：用 MATLAB 把应变花原始读数一路算到主应力方向，并且给出可验证的中间结果。核心关键词就是应变花数据、主应力方向、MATLAB 实现和计算验证。它适合三类人参考：一是做结构测试、需要现场快速出结果的技术人员；二是刚接触应变测量、想把公式和代码对应起来的学生或初级工程师；三是已经有现成脚本、但想核对一下自己角度约定有没有搞错的老手。整套东西不依赖任何工具箱，基础 MATLAB 就能跑，代码量不大，但每一个中间量我都会摊开讲清楚，方便你逐行对照自己的数据。

我写这套东西的出发点很朴素：应变花计算不难，难的是“算完之后你敢不敢信”。所以整篇的重点不只是把公式敲进 MATLAB，而是把验证环节做扎实——用已知解析解反算、用不同贴片角度交叉核对、用莫尔圆做几何校验，三条路互相印证。下面从整体设计思路开始拆。

2. 整体设计与思路拆解：为什么这样组织计算流程

2.1 从原始读数到主应力的链路拆解

应变花计算本质上是一条固定的数据链路：原始应变读数 → 应变分量（正应变与剪应变）→ 主应变 → 主应力 → 主方向角。每一步都有明确的数学关系，但每一步也都有容易出错的约定问题。我在组织代码时坚持一个原则：中间量全部显式保留，不做链式嵌套。也就是说，不写成一行长表达式直接出结果，而是把 εx、εy、γxy、ε1、ε2、σ1、σ2、θ 逐个算出来、逐个存变量。这样做的好处很直接——出问题时你能立刻定位是哪一环偏了，而不是面对一个错误结果干瞪眼。

具体到常见的三片式应变花（0°、45°、90° 布置），核心换算关系是这样的：设三个方向的读数分别为 ε0、ε45、ε90，那么

• εx = ε0
• εy = ε90
• γxy = 2·ε45 − ε0 − ε90

这一步是整个链路里最容易被轻视的地方。很多人以为 γxy 就是 ε45 减一减，实际上 45° 方向的应变里混着正应变和剪应变两部分贡献，系数 2 不能丢。我见过把系数漏掉导致剪应变直接减半的案例，主方向角跟着偏了将近十几度，现场还以为是贴片贴歪了，折腾半天。

拿到 εx、εy、γxy 之后，主应变按标准公式走：

• ε1,2 = (εx + εy)/2 ± sqrt( ((εx − εy)/2)² + (γxy/2)² )

主应力则通过平面应力状态的广义胡克定律转换，需要材料弹性模量 E 和泊松比 ν：

• σ1 = (E/(1−ν²))·(ε1 + ν·ε2)
• σ2 = (E/(1−ν²))·(ε2 + ν·ε1)

主方向角 θ 用下式求：

• tan(2θ) = γxy / (εx − εy)

这里必须强调：atan 出来的角度只是 2θ 的主值，落在 (−90°, 90°]，对应 θ 落在 (−45°, 45°]，它只告诉你方向线的走向，不告诉你主应力是拉还是压、指向哪一侧。这是后面象限判断那一节要专门处理的问题，也是我认为整套代码里最值得花篇幅讲清楚的部分。

2.2 为什么选 MATLAB 而不是别的工具

选 MATLAB 的理由不是“它流行”，而是这个任务的特点决定的。应变花计算涉及矩阵化的批量处理——现场往往几十上百个测点，每个测点一组读数，用脚本语言逐点循环写起来啰嗦，MATLAB 的向量化天然合适；同时它内置的绘图能力让莫尔圆校验几乎零成本，画个圆、标个点，几何关系一眼就能看出来对不对。相比之下，用通用编程语言做同样的事，光是画图和数据对齐就要多写不少胶水代码。

另一个考虑是可复现性。MATLAB 脚本是纯文本，参数集中在开头，换一组数据只改几行，别人拿到你的脚本能直接跑。我习惯把材料参数、贴片角度、数据文件路径全部放在脚本顶部的配置区，下面计算区不出现任何魔法数字。这个习惯在多人协作或者隔几个月回头复查时特别值钱——你不需要回忆“这个 0.3 是哪来的”，配置区写得明明白白。

至于要不要用工具箱，我的建议是能不用就不用。应变花计算用到的全是基础线性代数和三角函数，基础 MATLAB 完全覆盖。依赖越少，脚本在别人机器上跑起来的成功率越高，这一点在跨团队交付时体现得尤其明显。

2.3 验证思路的顶层设计

前面说了，这套东西的重点在验证。我设计了三条互相独立的校验路径，任何一条对不上都说明中间某步有问题：

第一条是解析解反算。构造一组已知主应变和主方向的“标准答案”，正算出三个方向的读数，再喂给主计算流程，看能不能还原出原来的主应变和方向。这条路径验证的是公式实现本身。

第二条是角度交叉核对。同一物理状态用不同贴片角度组合（比如 0-45-90 和 0-60-120）分别计算，理论上主应力大小和方向必须一致。这条路径验证的是角度约定有没有搞错。

第三条是莫尔圆几何校验。把 εx、εy、γxy/2 画在莫尔圆上，主应变点应该在圆与横轴的交点，主方向角对应圆心到点的连线角度。这条路径用几何直觉兜底，专治“数字看着差不多但方向反了”这类隐蔽错误。

三条路径的定位不同：第一条查公式，第二条查约定，第三条查直觉。实际排查时我一般先跑第一条，因为它最快、最确定；第一条过了再跑第二条，确认角度体系没乱；最后画莫尔圆，给自己一个视觉上的确认。下面几节就按这个逻辑展开。

3. 核心细节解析与实操要点：公式背后的坑

3.1 应变花类型与角度约定的统一

应变花不止 0-45-90 一种。常见的还有 0-60-120 的等角应变花，以及各种非标准角度布置。不同厂家的贴片，丝栅方向的定义也可能有细微差别。这就带来一个现实问题：你手上的读数，到底对应哪个物理方向？

我的处理办法是在代码里把“贴片角度”做成显式输入，而不是硬编码成 0、45、90。对标准三片式，输入向量就是 [0, 45, 90]；对等角式，输入 [0, 60, 120]。然后计算剪应变时用通用公式而不是特例公式。通用形式下，任意三个角度 θa、θb、θc 的读数可以解出 εx、εy、γxy，本质上是解一个三元一次方程组，每个方程形如：

ε(θ) = εx·cos²θ + εy·sin²θ + γxy·sinθ·cosθ

三个方程三个未知数，写成矩阵形式直接左除即可。这样做的好处是一套代码通吃所有贴片类型，而且角度约定集中在输入向量里，改起来一目了然。特例公式（比如 45° 那个系数 2）虽然算得快，但一旦贴片类型换了就得改代码，反而容易出错。

注意：角度约定的正方向必须全流程统一。我习惯以 x 轴为 0°、逆时针为正，贴片角度、计算出的主方向角、绘图时的角度标注全部遵循这一套。中途换约定是方向角出错的头号原因。

3.2 剪应变符号与主方向象限判断

这是整套计算里我最想强调的部分。前面提到，atan 给出的 θ 落在 (−45°, 45°]，它只描述方向线的走向。但主应力是有拉压之分的，方向角应该能区分“沿这个方向是拉”还是“沿这个方向是压”。举个具体例子：假设算出 θ = 30°，那么主应力方向线是 30° 这条线，但 σ1 可能沿 30° 方向，也可能沿 30°+90° = 120° 方向，取决于哪个是较大主应力。

正确的处理方式是：先算出两个主应力 σ1、σ2（σ1 取较大者），然后判断 σ1 对应的方向。用公式 θ1 = 0.5·atan2(γxy, εx − εy) 可以得到一个候选角，但 atan2 的象限处理和 2θ 的周期性需要小心。更稳妥的做法是回代验证：把候选角 θ 代入应变转换公式，算出该方向的法向应变，看它是否等于 ε1。如果等于 ε2，说明方向差了 90°，加或减 90° 修正即可。这种“算出来再代回去验”的思路看起来笨，但它对符号约定不敏感，几乎不会错。

我在实际项目里踩过的坑是：某批数据剪应变符号整体反了（贴片接线极性搞反），atan 算出来的角度看着“挺合理”，但回代一验就露馅——该方向的法向应变对不上 ε1。所以回代验证这一步我建议无条件保留，成本极低，收益极高。

3.3 材料参数与平面应力假设的边界

主应力换算依赖 E 和 ν，这两个参数填错，应力值整体偏移，但方向角不受影响（因为方向只由应变分量比例决定）。这一点值得说清楚：方向角对材料参数不敏感，应力大小对材料参数高度敏感。所以如果你的关注点是方向，材料参数可以先用标称值；如果关注应力绝对值，参数必须用实测或可靠来源。

平面应力假设也有边界。应变花贴在结构表面，表面法向应力为零，平面应力假设成立，这是标准做法。但如果测点附近有应力集中、或者结构很厚导致表面附近仍有明显三向应力，平面应力假设就会引入误差。这种情况下应变花本身也给不出完整的三向应力状态，需要补充其他测量手段。我在代码注释里会明确写一句“本流程基于平面应力假设”，提醒后来者别在不适用场景硬套。

3.4 单位与量纲的统一

应变读数通常是微应变（με），而 E 常用 GPa，两者量纲差着数量级。我的做法是在数据入口统一换算成无量纲应变（με 除以 1e6），E 统一换算成 Pa，这样中间计算全程 SI 单位，最后输出应力时再按需要转成 MPa。单位混乱是数值结果离谱的常见原因，尤其是从不同来源拼数据的时候。代码里我会在配置区写一行注释标明“输入应变单位：με；内部已转无量纲”，避免后来者重复换算。

4. 实操过程与核心环节实现：从脚本到结果

4.1 配置区与数据入口的写法

脚本开头我固定放一个配置区，把所有外部输入集中管理。结构大致是这样：材料参数（E、ν）、贴片角度向量、数据文件路径、输出单位选项。数据我用一个简单的文本或 CSV 文件读入，每行一个测点，列依次是三个方向的读数。读入后用矩阵运算一次性处理所有测点，避免逐点循环。

% ===== 配置区 ===== E = 210e9; % 弹性模量，Pa nu = 0.30; % 泊松比 ang = [0, 45, 90]; % 贴片角度，度，逆时针为正 dataFile = 'strain_data.csv'; % 每行一个测点，三列读数，单位 με outUnit = 'MPa'; % 输出应力单位 % ===== 数据读入 ===== raw = readmatrix(dataFile); % N x 3，单位 με epsMeas = raw * 1e-6; % 转无量纲应变

这段代码没什么花哨的，但“配置区与计算区分离”这个习惯我强烈建议保留。现场改参数只动上面几行，下面计算逻辑一个字不用碰，出错概率大幅下降。

4.2 应变分量求解的向量化实现

贴片角度通用的情况下，每个测点要解一个 3×3 线性方程组。向量化的写法是把系数矩阵预先构造好，然后对整批数据做一次矩阵左除。系数矩阵的每一行对应一个贴片角度，元素是 [cos²θ, sin²θ, sinθ·cosθ]，未知量是 [εx, εy, γxy]。

% 构造系数矩阵 A（3x3），与测点无关，只与贴片角度有关 th = deg2rad(ang(:)); A = [cos(th).^2, sin(th).^2, sin(th).*cos(th)]; % 对每个测点求解：epsMeas 是 N x 3，转置后左除 % A * [ex; ey; gxy] = [e_a; e_b; e_c] comp = (A \ epsMeas.').'; % N x 3，列依次为 ex, ey, gxy ex = comp(:,1); ey = comp(:,2); gxy = comp(:,3);

这里有个细节值得说：A \ epsMeas.'里 A 是 3×3，右边是 3×N，左除一次就把所有测点解完了，比循环快得多，代码也更短。解出来的 γxy 是工程剪应变（不是张量剪应变），后面莫尔圆校验时要注意除以 2 才是张量分量，这个区别在绘图那节会再提。

4.3 主应变、主应力与方向角的计算

拿到 εx、εy、γxy 之后，主应变和主应力都是逐元素的向量运算，写起来很直接。方向角我按前面说的“atan2 候选 + 回代验证”两步走。

% 主应变 epsAvg = (ex + ey) / 2; epsRad = sqrt(((ex - ey)/2).^2 + (gxy/2).^2); e1 = epsAvg + epsRad; e2 = epsAvg - epsRad; % 主应力（平面应力） coef = E / (1 - nu^2); s1 = coef * (e1 + nu * e2); s2 = coef * (e2 + nu * e1); % 方向角候选（弧度），atan2 处理象限 thetaCand = 0.5 * atan2(gxy, ex - ey); % 回代验证：算候选方向的法向应变，与 e1 比对 epsAtCand = ex.*cos(thetaCand).^2 + ey.*sin(thetaCand).^2 ... + gxy.*sin(thetaCand).*cos(thetaCand); % 若与 e2 更接近，说明方向差 90 度，修正 needFlip = abs(epsAtCand - e2) < abs(epsAtCand - e1); theta1 = thetaCand; theta1(needFlip) = theta1(needFlip) + pi/2; theta1 = rad2deg(theta1);

回代验证那几行是整个脚本里我最看重的部分。它不依赖任何符号约定的假设，纯粹用“这个方向算出来的法向应变应该等于较大主应变”这一物理事实做判断，鲁棒性很好。needFlip是逻辑向量，批量处理时一次修正所有需要翻转的测点，不用循环。

4.4 一个可复现的验证算例

光有代码不够，得有标准答案。我构造一个已知状态：设 ε1 = 800 με、ε2 = −200 με、主方向 θ1 = 30°（即较大主应变沿 30° 方向）。用应变转换公式正算出 0°、45°、90° 三个方向的读数，作为“实测数据”喂给上面的流程，看能否还原。

正算过程：ε(θ) = ε1·cos²(θ−θ1) + ε2·sin²(θ−θ1)。代入 θ = 0°、45°、90°，得到三个读数（具体数值跑一遍代码就有，这里不列死数，避免你直接抄结果而不理解过程）。把这三个数写进数据文件，跑主流程，预期输出 ε1 ≈ 800 με、ε2 ≈ −200 με、θ1 ≈ 30°。如果方向角出来是 −60° 或 120°，说明象限处理有问题；如果 ε1、ε2 对调，说明主次判断反了。这个算例我建议每个改过代码的人都跑一遍，它是回归测试的基线。

4.5 莫尔圆校验的绘制方法

莫尔圆校验的价值在于几何直觉。以 ε 为横轴、γ/2 为纵轴，圆心在 (εx+εy)/2，半径就是前面算的 epsRad。主应变点是圆与横轴的两个交点，主方向角对应从圆心到 (εx, γxy/2) 这条半径与横轴的夹角的两倍关系。

% 取第一个测点画莫尔圆 k = 1; cx = (ex(k) + ey(k)) / 2; R = sqrt(((ex(k)-ey(k))/2)^2 + (gxy(k)/2)^2); phi = linspace(0, 2*pi, 200); plot(cx + R*cos(phi), R*sin(phi), 'LineWidth', 1.2); hold on; plot([ex(k), ey(k)], [gxy(k)/2, -gxy(k)/2], 'o'); plot([cx-R, cx+R], [0, 0], 'x'); axis equal; grid on;

画出来之后，两个主应变交点是否落在横轴上、测点是否在圆上、半径是否等于 epsRad，一眼可辨。我遇到过剪应变符号反了的情况，莫尔圆上的点跑到横轴另一侧，视觉上非常明显，比盯数字快得多。这张图我一般会存进报告，作为计算可信度的佐证。

5. 常见问题与排查技巧实录

5.1 方向角差 90 度：最常见的“假错误”

新手最容易慌的情况是：算出来的方向角和预期差了 90°。先别急着改公式，这大概率不是错误，而是主次方向的定义问题。主应力方向是一条线，θ 和 θ+90° 描述的是同一条线的两个方向，哪个是 σ1 取决于你取哪个为主。判断方法就是回代验证：看该方向法向应变对应的是 ε1 还是 ε2。如果对应 ε2，加 90° 就对了。我在代码里用needFlip自动处理，但人工核对时也要理解这个逻辑，否则看到 90° 的差异会误以为代码坏了。

5.2 剪应变整体偏小或偏大：系数与单位问题

如果主应力大小系统性偏离预期，但方向角正常，优先查两处：一是 45° 方向剪应变公式里的系数 2 有没有漏（用通用矩阵解法时这个坑自动规避）；二是单位换算有没有重复或遗漏，比如数据已经是无量纲却被又乘了一次 1e-6。方向角正常说明应变分量的比例关系大致对，问题多半出在绝对量纲上，顺着单位链条捋一遍基本能定位。

5.3 个别测点结果离谱：数据质量与贴片问题

整批数据里个别点结果异常，先别怀疑算法，查数据本身。常见原因包括：读数超出量程、贴片脱胶导致某方向读数漂移、接线极性问题导致某一路符号反了。排查办法是把该测点的三个原始读数单独拎出来，手工算一遍 εx、εy、γxy，再画个莫尔圆，看三个点是否构成合理的几何关系。如果某个读数明显不合群，基本可以判定是测量环节的问题，而不是计算环节。

5.4 常见问题速查表

最后一行值得多说一句：当 εx ≈ εy 且 γxy ≈ 0 时，主方向在数学上是不确定的（莫尔圆退化成点），此时方向角的微小扰动会被放大，算出来的角度跳变是正常现象，不代表代码有问题。这种测点报告里应注明“方向不确定”，而不是硬报一个角度。

5.5 几条踩坑换来的实操心得

第一条，先跑验证算例再碰真实数据。我现在的习惯是每次改动脚本，先跑 4.4 节那个已知答案的算例，通过了才处理现场数据。这个习惯帮我挡掉过好几次手滑改错公式的事故。

第二条，中间量存下来别扔。εx、εy、γxy、ε1、ε2 全部保留在结果表里，报告里不一定全放，但排查时它们是第一手证据。只存最终应力和方向，出了问题等于把线索全丢了。

第三条，角度基准写在报告里。方向角是相对哪个基准、正方向怎么定义，必须在报告里白纸黑字写清楚。我见过两份报告因为基准不同，同一个测点的方向角差了 90°，对接时扯了半天，其实两边都没算错。

第四条，批量数据先看分布再下结论。几十个测点的方向角，先看整体分布是否合理，再看个别异常。如果整体分布就乱，多半是系统性约定问题；如果整体合理个别跳，多半是数据质量问题。这个判断顺序能省很多排查时间。

这套流程我用了挺久，从单点手算到批量脚本，最大的体会是：应变花计算的数学部分其实很轻，真正花时间的是约定统一和结果验证。把这两块做扎实，剩下的就是改改参数跑数据的事。后面如果测点规模再大，可以考虑把验证算例做成自动化回归测试，每次跑数据前自动执行一遍，进一步降低人为失误的空间。