简介这套MATLAB代码实现了二维离散位错动力学DDD模拟核心功能是计算滑移面上的应力分布适合从事材料力学、晶体塑性或位错动力学研究的科研人员和学生使用。压缩包共26个文件以18个.m脚本为主体辅以3个txt输入文件、2个fig与2个png图像及1份md说明文档整体仅36KB结构紧凑轻量。程序通过dd2d脚本启动读取滑移面、位错列表等输入文件后可模拟单个滑移面上的位错运动并计算沿滑移面的应力变化代码覆盖位错应力场、Peach-Koehler力、时间增量、位错源创建与滑移面绘制等关键模块同时包含fig/png结果图可供可视化参考便于二次开发与扩展至多滑移面滑移系统及晶界应力分析。目前已有149人学习适合需要快速入门位错动力学数值模拟或正在探索滑移系与晶界效应的研究者下载参考。1. DDD 这个词有两副面孔材料学里它叫离散位错动力学搜索“DDD”时十个人里有九个是在找领域驱动设计剩下一个要找的才是本文要讲的——Discrete Dislocation Dynamics离散位错动力学。纳米压痕、微柱压缩、辐照损伤这些单晶尺度问题里位错源开动、位错偶极子湮灭都是离散事件连续介质塑性模型解释不了于是就需要二维 DDD把刃型位错当作平面上的点缺陷让它们在长程应力场中相互作用和运动。“测量沿滑移面的应力分布”就是在模拟域中选定一条滑移线用位错弹性场解析叠加算出该线上各采样点的剪应力再把它喂给 Peach-Koehler 力去更新位错位置。下面的 MATLAB 脚本拆成两个 .m 文件放同目录即可直接运行新手能照着改熟手也能看到采样间距、核心截断和边界处理上的取舍。2. 二维 DDD 的应力场叠加与滑移面剪应力投影2.1 刃型位错的解析应力场二维 DDD 里每条位错都被视为垂直于建模平面的直线刃型位错弹性问题退化为平面应变问题。线弹性理论支持应力叠加任意观测点的总应力等于所有位错贡献的线性求和。这一步没有迭代求解代价是 O(N_q × N_d)采样点数乘位错条数。对一条 Burgers 矢量为 b、滑移方向沿 u 轴、滑移面法向沿 v 轴的刃型位错在局部坐标系 (u, v) 中应力分量为σuu -D0 · b · v · (3u² v²) / R⁴σvv D0 · b · v · (u² − v²) / R⁴σuv D0 · b · u · (u² − v²) / R⁴其中 R² u² v²D0 G / [2π(1−ν)]G 是剪切模量ν 是泊松比。实际位错的 Burgers 矢量方向与全局 x 轴往往有夹角 θ所以计算时先把全局坐标平移到位错处再旋转到局部系最后把局部应力张量旋转回全局系。常见实现会把 θ 作为位错结构体的一个字段而不是每次都重新推导旋转矩阵。2.2 应力张量到滑移面剪应力的投影“沿滑移面的应力分布”并不是直接把 σxx 或 σyy 画出来而是把应力张量投影到目标滑移面的滑移方向 m 和法向 n 上。设滑移方向的方位角为 θs则 m(cosθs, sinθs)n(−sinθs, cosθs)即顺时针旋转的滑移面法向则该面上的剪应力为τ mᵢ σᵢⱼ nⱼ sinθs·cosθs·(σyy − σxx) (cos²θs − sin²θs)·σxy当滑移面平行于 x 轴时 θs0τσxy滑移面垂直于 x 轴时 θsπ/2τ−σxy。用这两个极限检查投影代码是否漏项非常方便。多滑移系模拟中不同位错滑移系的 θs 不同测量每条滑移面时必须用各自的方向角不能偷懒取成位错自身的 θ。2.3 长程性、边界与平面应变假设刃型位错的远场应力按 1/R 衰减比电偶极子慢得多所以几十倍位错间距外的贡献仍然不能忽略。二维 DDD 的解析解默认是无界线弹性介质如果模拟域边界代表自由表面需要在域外放镜像位错来近似满足表面零牵引力条件否则“沿滑移面的应力分布”在边界附近会出现明显虚假值。周期性边界则用有限阶周期镜像叠加。测量之前先想清楚你要的分布究竟是无限域位错相互作用的结果还是带真实边界条件的场。两者在滑移面中段可能接近但在边界附近能差出几个数量级。3. MATLAB 代码沿滑移面测应力分布的完整脚本3.1 代码文件组织把实现拆成应力计算函数和测量主脚本两个文件disl_stress.m负责坐标变换、核心截断和应力张量旋转ddd_slip_stress.m负责构造位错阵列、生成滑移面采样线、投影剪应力和绘图。这样的拆分好处在于后续只改主脚本就能批量测不同位置的滑移面核心应力函数不用动。3.2 应力场计算函数disl_stress.mfunction [sxx, syy, sxy] disl_stress(xq, yq, disl, G, nu, rc) % 计算一组直线刃型位错在查询点处的面内应力分量 % xq, yq: 查询点坐标长度相同的列向量或矩阵 % disl: 位错结构体数组字段 x, y, b, theta % G, nu: 剪切模量与泊松比 % rc: 位错核心截断半径默认可不传 if nargin 6 rc 1e-9; end xq xq(:); yq yq(:); sxx zeros(size(xq)); syy zeros(size(xq)); sxy zeros(size(xq)); D0 G / (2 * pi * (1 - nu)); for k 1:numel(disl) dx xq - disl(k).x; dy yq - disl(k).y; % 局部坐标u 沿 Burgers 矢量方向v 沿滑移面法向 c cos(disl(k).theta); s sin(disl(k).theta); u c .* dx s .* dy; v -s .* dx c .* dy; r2 u.^2 v.^2; r2 max(r2, rc^2); % 越过位错核心的采样点按截断处理 r4 r2.^2; b disl(k).b; % 刃型位错局部应力场 suu -D0 * b .* v .* (3*u.^2 v.^2) ./ r4; svv D0 * b .* v .* ( u.^2 - v.^2) ./ r4; suv D0 * b .* u .* ( u.^2 - v.^2) ./ r4; % 二阶张量旋转回全局坐标 c2 c^2; s2 s^2; cs c*s; sxx sxx c2.*suu s2.*svv - 2*cs.*suv; syy syy s2.*suu c2.*svv 2*cs.*suv; sxy sxy cs.*(suu - svv) (c2 - s2).*suv; end end函数内部先逐条位错平移坐标再做局部旋转最后用二阶张量旋转公式把局部应力转回全局。注意r2 max(r2, rc^2)这一行它保证采样点即使落到位错核心附近也不会触发除零但代价是核心内部的应力被截断成某个有限值因此小于 rc 范围内的数值不具备物理意义。3.3 测量主脚本ddd_slip_stress.m% ddd_slip_stress.m % 测量二维离散位错动力学中沿指定滑移面的应力分布 clear; clc; % 材料参数 G 26e9; % 剪切模量单位 Pa nu 0.33; % 泊松比 rc 1.0e-9; % 核心截断半径取 Burgers 矢量约 3.5 倍 % 构造三条刃型位错坐标、带符号 Burgers 大小、滑移方向角 disl(1) struct(x, 0e-9, y, 0e-9, b, 2.86e-10, theta, 0); disl(2) struct(x, 50e-9, y, 0e-9, b, -2.86e-10, theta, 0); disl(3) struct(x, 25e-9, y, 30e-9, b, 2.86e-10, theta, pi/6); % 生成滑移面采样线测线取 y 20 nm 的水平面 samp_x linspace(-200e-9, 200e-9, 401); samp_y 20e-9 * ones(size(samp_x)); % 计算该测线上每个采样点的应力张量分量 [sxx, syy, sxy] disl_stress(samp_x, samp_y, disl, G, nu, rc); % 投影到目标滑移面这里取 θs0退化为 x 方向滑移面 theta_s 0; tau sin(theta_s) .* cos(theta_s) .* (syy - sxx) ... (cos(theta_s).^2 - sin(theta_s).^2) .* sxy; % 绘图 figure(Color, w); plot(samp_x * 1e9, tau / 1e6, b-, LineWidth, 1.2); xlabel(沿滑移面位置 (nm)); ylabel(滑移面上剪应力 (MPa)); grid on;主脚本里第 1 条和第 2 条位错构成一个刃型位错偶极子第 3 条位错以 π/6 角度倾斜。此时若把theta_s改为 π/6测的就是第三条位错所在滑移面的剪应力分布。samp_x的 401 个点用linspace生成间距约 1 nm远小于后面提到的 rc 的 1/3能保留近场陡峭峰的形状。4. 三个决定滑移面应力分布精度的参数4.1 核心截断半径 rc 怎么取rc 的取值直接影响近核心区的应力峰值。取太小网格点一旦落到 R≈0 就出现虚假的巨峰取太大真实峰会被削平。常见做法是取 3b5bb 是 Burgers 矢量大小。对大多数金属 b≈0.25 nmrc 在 0.751.25 nm 量级。应力测量时还要确认采样线离最近位错的距离应大于 rc否则该段曲线不是弹性解而是截断函数的人为产物。如果关心的不是峰值而是远场趋势rc 的影响会明显减弱。对 1/R 衰减的位错场只要测线与位错的距离大于 10rc相对误差就能降到 1% 以下。4.2 采样间距和测线长度采样间距决定应力峰能否被捕捉。位错在滑移面上产生的剪应力峰宽度大约与离位错的距离同量级如果 dx 大于峰宽的一半峰顶误差就会超过 10%。建议 dx 不超过 rc/3测线长度覆盖若干倍的平均位错间距通常取 20 倍以上。测线长度不足时曲线两端看起来像被“切断”这不是物理现象而是采样域截断。加长测线再裁剪显示区间比直接缩短测线更保险。4.3 边界条件相关的镜像参数无限域叠加结果只适用于完全无视边界的理想情况。做自由表面或异质界面附近的测量时要补镜像位错。镜像阶数越高边界零牵引力近似越好但计算量线性增加。下表给出常用参数范围参数符号建议范围对结果的影响核心截断半径rc3b5b近核心区峰值与曲线形态采样间距dx≤ rc/3峰值捕捉能力测线长度Ls≥ 20×平均位错间距远场截断误差周期镜像阶数Nimg12 阶周期边界近似误差位错总数N几十到几千计算耗时 O(Nq×N)镜像数超过 2 阶后对滑移面中段应力分布的影响通常小于 0.1%但对紧贴表面的测线仍可能有可感知的偏差。先用 1 阶算一遍再加到 2 阶对比两层结果几乎重合就说明边界影响已收敛。5. 验证与排错用解析解校核 MATLAB 应力分布5.1 单一位错的解析验证最可靠的校验是把位错数降到 1在任意不经过位错核心的测线上比较数值结果与解析式。以 σxy 为例代码应精确复现下式σxy D0 · b · x · (x² − y²) / (x² y²)²xq linspace(-200e-9, 200e-9, 501); yq 50e-9 * ones(size(xq)); disl.x 0; disl.y 0; disl.b 2.86e-10; disl.theta 0; [~, ~, sxy_num] disl_stress(xq, yq, disl, 26e9, 0.33, 0); x xq; y yq; D0 26e9 / (2 * pi * (1 - 0.33)); r2 x.^2 y.^2; sxy_ana D0 * disl.b .* x .* (x.^2 - y.^2) ./ r2.^2; max(abs(sxy_num - sxy_ana)) / max(abs(sxy_ana))这里把 rc 设成 0同时让 yq 远离位错所在的 y0 平面就不会触碰奇异点。相对误差输出应该是 1e-12 量级说明坐标旋转、张量旋转和符号全部正确。5.2 常见错误与排查顺序最容易犯的错误是局部坐标系的 u、v 定义反了。u 必须沿 Burgers 矢量方向v 是滑移面法向方向错了张量旋转会引入符号错误导致 τ 的峰出现左右不对称。第二个常见问题是 rc 与 r4 不一致有些人先计算 r2 再用 max 截断却忘了用截断后的值重新计算 r4近核心区会出现大数除以大数的跳变。第三个是角度制度MATLAB 的 sin/cos 只接受弧度把 30° 直接传入而不转成弧度投影结果不会报错但完全不可用。排查时先跑 5.1 的解析验证再把位错偶极子情况手算几个特征点。例如同号位错沿滑移面对称分布时τ 曲线应该呈对称双峰反号位错则会形成一个正负交替的偶极子响应。符号和峰值位置对不上问题几乎都出在 Burgers 矢量符号或旋转方向。5.3 热图与曲线结合的可视化曲线能告诉你沿滑移面的具体大小热图能帮你判断滑移面位置是否选得合理。做法是先在整个域计算 σxy 场xg linspace(-200e-9, 200e-9, 120); yg linspace(-200e-9, 200e-9, 120); [XX, YY] meshgrid(xg, yg); [~, ~, sxyG] disl_stress(XX(:), YY(:), disl, G, nu, rc); pcolor(XX*1e9, YY*1e9, reshape(sxyG, size(XX))); shading flat; colorbar; hold on;热图上能看到每个位错核附近的“蝴蝶形”高应力区再用plot(samp_x*1e9, tau/1e6)把滑移面曲线叠上去。曲线跨过位错偶极子所在位置时应能看到过零点和衰减尾部这比单独看曲线更容易定位是哪条位错贡献了哪个峰。6. 进阶把滑移面应力分布接到 Peach-Koehler 力上6.1 位错自身位置的应力采样分布曲线不只是为了画图它最终要转换为位错运动。Peach-Koehler 公式给出位错受力F (σ · b) × ξ。对垂直于 xy 平面的直线刃型位错 ξ(0,0,1)面内分量为Fx σxy·bx σyy·byFy −(σxx·bx σxy·by)注意这里必须排除位错自身的应力贡献。单条位错的自应力在本身位置处是奇异的物理上也不产生自驱动力。实现时在应力叠加循环里用if k ~ self跳过自己或者单独计算“外部位错对该位错位置贡献的应力”。6.2 从剪应力到位错滑移速度滑移方向的 Peach-Koehler 力密度 Fs τs·bτs 是滑移面上作用于该位错的分解剪应力b 是 Burgers 矢量大小。得到 Fs 后常见做法是用幂律滑移律更新速度v v0 · (|Fs| / F0)^n · sign(Fs)其中 v0 是参考速度F0 是参考力密度n 常取 1050。时间步长 Δt 需要保证单步内位错移动距离小于位错间距的 10%否则滑动几个采样间隔后应力分布会剧烈震荡。6.3 批处理与脚本化执行做参数扫描时不要每次都打开交互式 MATLAB 窗口。把主脚本改成函数形式例如ddd_slip_stress(theta_disl, rc_val)在系统终端里用matlab -batch ddd_slip_stress(0, 1e-9)批量执行用 Codex 这类编码代理协助改参数时脚本化的运行方式也更可靠因为无头环境下不会因图窗阻塞而死等。批处理时把采样线固定成 401 点向量所有参数组的 tau 输出就能直接堆成三维数组后面做均值、方差和极值统计都不需要重新计算应力场。本文还有配套的精品资源点击获取
离散位错动力学中的滑移面应力分布测量：MATLAB实现与参数分析
简介这套MATLAB代码实现了二维离散位错动力学DDD模拟核心功能是计算滑移面上的应力分布适合从事材料力学、晶体塑性或位错动力学研究的科研人员和学生使用。压缩包共26个文件以18个.m脚本为主体辅以3个txt输入文件、2个fig与2个png图像及1份md说明文档整体仅36KB结构紧凑轻量。程序通过dd2d脚本启动读取滑移面、位错列表等输入文件后可模拟单个滑移面上的位错运动并计算沿滑移面的应力变化代码覆盖位错应力场、Peach-Koehler力、时间增量、位错源创建与滑移面绘制等关键模块同时包含fig/png结果图可供可视化参考便于二次开发与扩展至多滑移面滑移系统及晶界应力分析。目前已有149人学习适合需要快速入门位错动力学数值模拟或正在探索滑移系与晶界效应的研究者下载参考。1. DDD 这个词有两副面孔材料学里它叫离散位错动力学搜索“DDD”时十个人里有九个是在找领域驱动设计剩下一个要找的才是本文要讲的——Discrete Dislocation Dynamics离散位错动力学。纳米压痕、微柱压缩、辐照损伤这些单晶尺度问题里位错源开动、位错偶极子湮灭都是离散事件连续介质塑性模型解释不了于是就需要二维 DDD把刃型位错当作平面上的点缺陷让它们在长程应力场中相互作用和运动。“测量沿滑移面的应力分布”就是在模拟域中选定一条滑移线用位错弹性场解析叠加算出该线上各采样点的剪应力再把它喂给 Peach-Koehler 力去更新位错位置。下面的 MATLAB 脚本拆成两个 .m 文件放同目录即可直接运行新手能照着改熟手也能看到采样间距、核心截断和边界处理上的取舍。2. 二维 DDD 的应力场叠加与滑移面剪应力投影2.1 刃型位错的解析应力场二维 DDD 里每条位错都被视为垂直于建模平面的直线刃型位错弹性问题退化为平面应变问题。线弹性理论支持应力叠加任意观测点的总应力等于所有位错贡献的线性求和。这一步没有迭代求解代价是 O(N_q × N_d)采样点数乘位错条数。对一条 Burgers 矢量为 b、滑移方向沿 u 轴、滑移面法向沿 v 轴的刃型位错在局部坐标系 (u, v) 中应力分量为σuu -D0 · b · v · (3u² v²) / R⁴σvv D0 · b · v · (u² − v²) / R⁴σuv D0 · b · u · (u² − v²) / R⁴其中 R² u² v²D0 G / [2π(1−ν)]G 是剪切模量ν 是泊松比。实际位错的 Burgers 矢量方向与全局 x 轴往往有夹角 θ所以计算时先把全局坐标平移到位错处再旋转到局部系最后把局部应力张量旋转回全局系。常见实现会把 θ 作为位错结构体的一个字段而不是每次都重新推导旋转矩阵。2.2 应力张量到滑移面剪应力的投影“沿滑移面的应力分布”并不是直接把 σxx 或 σyy 画出来而是把应力张量投影到目标滑移面的滑移方向 m 和法向 n 上。设滑移方向的方位角为 θs则 m(cosθs, sinθs)n(−sinθs, cosθs)即顺时针旋转的滑移面法向则该面上的剪应力为τ mᵢ σᵢⱼ nⱼ sinθs·cosθs·(σyy − σxx) (cos²θs − sin²θs)·σxy当滑移面平行于 x 轴时 θs0τσxy滑移面垂直于 x 轴时 θsπ/2τ−σxy。用这两个极限检查投影代码是否漏项非常方便。多滑移系模拟中不同位错滑移系的 θs 不同测量每条滑移面时必须用各自的方向角不能偷懒取成位错自身的 θ。2.3 长程性、边界与平面应变假设刃型位错的远场应力按 1/R 衰减比电偶极子慢得多所以几十倍位错间距外的贡献仍然不能忽略。二维 DDD 的解析解默认是无界线弹性介质如果模拟域边界代表自由表面需要在域外放镜像位错来近似满足表面零牵引力条件否则“沿滑移面的应力分布”在边界附近会出现明显虚假值。周期性边界则用有限阶周期镜像叠加。测量之前先想清楚你要的分布究竟是无限域位错相互作用的结果还是带真实边界条件的场。两者在滑移面中段可能接近但在边界附近能差出几个数量级。3. MATLAB 代码沿滑移面测应力分布的完整脚本3.1 代码文件组织把实现拆成应力计算函数和测量主脚本两个文件disl_stress.m负责坐标变换、核心截断和应力张量旋转ddd_slip_stress.m负责构造位错阵列、生成滑移面采样线、投影剪应力和绘图。这样的拆分好处在于后续只改主脚本就能批量测不同位置的滑移面核心应力函数不用动。3.2 应力场计算函数disl_stress.mfunction [sxx, syy, sxy] disl_stress(xq, yq, disl, G, nu, rc) % 计算一组直线刃型位错在查询点处的面内应力分量 % xq, yq: 查询点坐标长度相同的列向量或矩阵 % disl: 位错结构体数组字段 x, y, b, theta % G, nu: 剪切模量与泊松比 % rc: 位错核心截断半径默认可不传 if nargin 6 rc 1e-9; end xq xq(:); yq yq(:); sxx zeros(size(xq)); syy zeros(size(xq)); sxy zeros(size(xq)); D0 G / (2 * pi * (1 - nu)); for k 1:numel(disl) dx xq - disl(k).x; dy yq - disl(k).y; % 局部坐标u 沿 Burgers 矢量方向v 沿滑移面法向 c cos(disl(k).theta); s sin(disl(k).theta); u c .* dx s .* dy; v -s .* dx c .* dy; r2 u.^2 v.^2; r2 max(r2, rc^2); % 越过位错核心的采样点按截断处理 r4 r2.^2; b disl(k).b; % 刃型位错局部应力场 suu -D0 * b .* v .* (3*u.^2 v.^2) ./ r4; svv D0 * b .* v .* ( u.^2 - v.^2) ./ r4; suv D0 * b .* u .* ( u.^2 - v.^2) ./ r4; % 二阶张量旋转回全局坐标 c2 c^2; s2 s^2; cs c*s; sxx sxx c2.*suu s2.*svv - 2*cs.*suv; syy syy s2.*suu c2.*svv 2*cs.*suv; sxy sxy cs.*(suu - svv) (c2 - s2).*suv; end end函数内部先逐条位错平移坐标再做局部旋转最后用二阶张量旋转公式把局部应力转回全局。注意r2 max(r2, rc^2)这一行它保证采样点即使落到位错核心附近也不会触发除零但代价是核心内部的应力被截断成某个有限值因此小于 rc 范围内的数值不具备物理意义。3.3 测量主脚本ddd_slip_stress.m% ddd_slip_stress.m % 测量二维离散位错动力学中沿指定滑移面的应力分布 clear; clc; % 材料参数 G 26e9; % 剪切模量单位 Pa nu 0.33; % 泊松比 rc 1.0e-9; % 核心截断半径取 Burgers 矢量约 3.5 倍 % 构造三条刃型位错坐标、带符号 Burgers 大小、滑移方向角 disl(1) struct(x, 0e-9, y, 0e-9, b, 2.86e-10, theta, 0); disl(2) struct(x, 50e-9, y, 0e-9, b, -2.86e-10, theta, 0); disl(3) struct(x, 25e-9, y, 30e-9, b, 2.86e-10, theta, pi/6); % 生成滑移面采样线测线取 y 20 nm 的水平面 samp_x linspace(-200e-9, 200e-9, 401); samp_y 20e-9 * ones(size(samp_x)); % 计算该测线上每个采样点的应力张量分量 [sxx, syy, sxy] disl_stress(samp_x, samp_y, disl, G, nu, rc); % 投影到目标滑移面这里取 θs0退化为 x 方向滑移面 theta_s 0; tau sin(theta_s) .* cos(theta_s) .* (syy - sxx) ... (cos(theta_s).^2 - sin(theta_s).^2) .* sxy; % 绘图 figure(Color, w); plot(samp_x * 1e9, tau / 1e6, b-, LineWidth, 1.2); xlabel(沿滑移面位置 (nm)); ylabel(滑移面上剪应力 (MPa)); grid on;主脚本里第 1 条和第 2 条位错构成一个刃型位错偶极子第 3 条位错以 π/6 角度倾斜。此时若把theta_s改为 π/6测的就是第三条位错所在滑移面的剪应力分布。samp_x的 401 个点用linspace生成间距约 1 nm远小于后面提到的 rc 的 1/3能保留近场陡峭峰的形状。4. 三个决定滑移面应力分布精度的参数4.1 核心截断半径 rc 怎么取rc 的取值直接影响近核心区的应力峰值。取太小网格点一旦落到 R≈0 就出现虚假的巨峰取太大真实峰会被削平。常见做法是取 3b5bb 是 Burgers 矢量大小。对大多数金属 b≈0.25 nmrc 在 0.751.25 nm 量级。应力测量时还要确认采样线离最近位错的距离应大于 rc否则该段曲线不是弹性解而是截断函数的人为产物。如果关心的不是峰值而是远场趋势rc 的影响会明显减弱。对 1/R 衰减的位错场只要测线与位错的距离大于 10rc相对误差就能降到 1% 以下。4.2 采样间距和测线长度采样间距决定应力峰能否被捕捉。位错在滑移面上产生的剪应力峰宽度大约与离位错的距离同量级如果 dx 大于峰宽的一半峰顶误差就会超过 10%。建议 dx 不超过 rc/3测线长度覆盖若干倍的平均位错间距通常取 20 倍以上。测线长度不足时曲线两端看起来像被“切断”这不是物理现象而是采样域截断。加长测线再裁剪显示区间比直接缩短测线更保险。4.3 边界条件相关的镜像参数无限域叠加结果只适用于完全无视边界的理想情况。做自由表面或异质界面附近的测量时要补镜像位错。镜像阶数越高边界零牵引力近似越好但计算量线性增加。下表给出常用参数范围参数符号建议范围对结果的影响核心截断半径rc3b5b近核心区峰值与曲线形态采样间距dx≤ rc/3峰值捕捉能力测线长度Ls≥ 20×平均位错间距远场截断误差周期镜像阶数Nimg12 阶周期边界近似误差位错总数N几十到几千计算耗时 O(Nq×N)镜像数超过 2 阶后对滑移面中段应力分布的影响通常小于 0.1%但对紧贴表面的测线仍可能有可感知的偏差。先用 1 阶算一遍再加到 2 阶对比两层结果几乎重合就说明边界影响已收敛。5. 验证与排错用解析解校核 MATLAB 应力分布5.1 单一位错的解析验证最可靠的校验是把位错数降到 1在任意不经过位错核心的测线上比较数值结果与解析式。以 σxy 为例代码应精确复现下式σxy D0 · b · x · (x² − y²) / (x² y²)²xq linspace(-200e-9, 200e-9, 501); yq 50e-9 * ones(size(xq)); disl.x 0; disl.y 0; disl.b 2.86e-10; disl.theta 0; [~, ~, sxy_num] disl_stress(xq, yq, disl, 26e9, 0.33, 0); x xq; y yq; D0 26e9 / (2 * pi * (1 - 0.33)); r2 x.^2 y.^2; sxy_ana D0 * disl.b .* x .* (x.^2 - y.^2) ./ r2.^2; max(abs(sxy_num - sxy_ana)) / max(abs(sxy_ana))这里把 rc 设成 0同时让 yq 远离位错所在的 y0 平面就不会触碰奇异点。相对误差输出应该是 1e-12 量级说明坐标旋转、张量旋转和符号全部正确。5.2 常见错误与排查顺序最容易犯的错误是局部坐标系的 u、v 定义反了。u 必须沿 Burgers 矢量方向v 是滑移面法向方向错了张量旋转会引入符号错误导致 τ 的峰出现左右不对称。第二个常见问题是 rc 与 r4 不一致有些人先计算 r2 再用 max 截断却忘了用截断后的值重新计算 r4近核心区会出现大数除以大数的跳变。第三个是角度制度MATLAB 的 sin/cos 只接受弧度把 30° 直接传入而不转成弧度投影结果不会报错但完全不可用。排查时先跑 5.1 的解析验证再把位错偶极子情况手算几个特征点。例如同号位错沿滑移面对称分布时τ 曲线应该呈对称双峰反号位错则会形成一个正负交替的偶极子响应。符号和峰值位置对不上问题几乎都出在 Burgers 矢量符号或旋转方向。5.3 热图与曲线结合的可视化曲线能告诉你沿滑移面的具体大小热图能帮你判断滑移面位置是否选得合理。做法是先在整个域计算 σxy 场xg linspace(-200e-9, 200e-9, 120); yg linspace(-200e-9, 200e-9, 120); [XX, YY] meshgrid(xg, yg); [~, ~, sxyG] disl_stress(XX(:), YY(:), disl, G, nu, rc); pcolor(XX*1e9, YY*1e9, reshape(sxyG, size(XX))); shading flat; colorbar; hold on;热图上能看到每个位错核附近的“蝴蝶形”高应力区再用plot(samp_x*1e9, tau/1e6)把滑移面曲线叠上去。曲线跨过位错偶极子所在位置时应能看到过零点和衰减尾部这比单独看曲线更容易定位是哪条位错贡献了哪个峰。6. 进阶把滑移面应力分布接到 Peach-Koehler 力上6.1 位错自身位置的应力采样分布曲线不只是为了画图它最终要转换为位错运动。Peach-Koehler 公式给出位错受力F (σ · b) × ξ。对垂直于 xy 平面的直线刃型位错 ξ(0,0,1)面内分量为Fx σxy·bx σyy·byFy −(σxx·bx σxy·by)注意这里必须排除位错自身的应力贡献。单条位错的自应力在本身位置处是奇异的物理上也不产生自驱动力。实现时在应力叠加循环里用if k ~ self跳过自己或者单独计算“外部位错对该位错位置贡献的应力”。6.2 从剪应力到位错滑移速度滑移方向的 Peach-Koehler 力密度 Fs τs·bτs 是滑移面上作用于该位错的分解剪应力b 是 Burgers 矢量大小。得到 Fs 后常见做法是用幂律滑移律更新速度v v0 · (|Fs| / F0)^n · sign(Fs)其中 v0 是参考速度F0 是参考力密度n 常取 1050。时间步长 Δt 需要保证单步内位错移动距离小于位错间距的 10%否则滑动几个采样间隔后应力分布会剧烈震荡。6.3 批处理与脚本化执行做参数扫描时不要每次都打开交互式 MATLAB 窗口。把主脚本改成函数形式例如ddd_slip_stress(theta_disl, rc_val)在系统终端里用matlab -batch ddd_slip_stress(0, 1e-9)批量执行用 Codex 这类编码代理协助改参数时脚本化的运行方式也更可靠因为无头环境下不会因图窗阻塞而死等。批处理时把采样线固定成 401 点向量所有参数组的 tau 输出就能直接堆成三维数组后面做均值、方差和极值统计都不需要重新计算应力场。本文还有配套的精品资源点击获取