基于MATLAB的离散位错动力学模拟：应力场与Peach-Koehler力解析

简介：面向金属塑性变形与位错动力学研究者的二维DDD（离散位错动力学）MATLAB工具包，用于计算滑移面上的应力分布并模拟位错运动。当前版本聚焦单个滑移面上的位错运动，作者正计划扩展为一系列相同滑移面的模拟，适合材料专业学生或科研人员快速搭建位错应力场分析环境。压缩包共26个文件，以18个.m脚本为主体，覆盖主程序dd2d、位错创建/读取、应力场计算、Peach-Koehler力求解和绘图等模块；另有3个txt输入文件、2个fig图形、2个png示意及1个README说明文档，整体仅36KB，轻量易读。已有149人学习下载。通过运行dd2d.m即可复现结果，配合slipPlane.txt、dislList.txt、dsourceList.txt可自定义滑移面位置、位错源与初始位错列表；fig与png图则可直观对照滑移面应力分布和Weertman构型，适合在此基础上继续开发多滑移系与晶界应力分析。

1. DDD 滑移面模拟的骨架：这个 MATLAB 项目怎么连起来

离散位错动力学（Discrete Dislocation Dynamics，DDD）和有限元、分子动力学都不一样，它把材料看作一根根可滑移的位错线集合，通过应力场驱动位错运动。这里的 DDD 指位错动力学，不是开发圈常说的领域驱动设计，同名的缩写经常被搜索引擎混为一谈，读代码时先分清。这个 MATLAB 项目的切入点很收敛：只在单个滑移面上模拟位错运动，同时把滑移面上的应力分布完整算出来，最后输出 simulated.png 和 weertman.png 两张对照图。

整套代码的入口只有一个：dd2d.m。它读入 slipPlane.txt、dislList.txt、dsourceList.txt 三个文本文件，初始化时就建立滑移面、位错数组、位错源三个对象，主循环里交替完成应力场叠加、Peach-Koehler 力计算、时间增量判定和位错源释放。对想从连续介质模拟转向微观位错机制的工程师来说，这种能直接改参数看曲线变化的 MATLAB 工程，比完整的三维 DDD 框架好拆得多，也更容易把应力场公式和代码对应起来。

2. 位错应力场与 Peach-Koehler 力：DDD 内核中的两个关键计算

2.1 二维刃型位错的解析应力场如何落到 MATLAB 函数

在二维 DDD 中，位错线被简化成 z 方向无限长的直线，位错运动只发生在 xy 平面。一根 Burgers 矢量沿 x 方向的刃型位错，在距位错核心 ((x, y)) 的点产生的平面应力分量有闭式解，其中对位错滑移最关键的是切应力 (\sigma_{xy})。系数 (D = \mu b / [2\pi(1-\nu)]) 把剪切模量 (\mu)、泊松比 (\nu) 和 Burgers 矢量模长 (b) 压成一个标量，避免每个应力分量重复写同一串物理量。项目里的 dislocationStressField.m 负责把这个解析场数值化，代码核心是以下三行：

function [sigma] = dislocationStressField(x, y, b, mu, nu) % 输入: x,y 观测点相对位错核心的坐标偏移，可传向量 % b: Burgers 矢量模长, mu: 剪切模量, nu: 泊松比 % 输出: 结构体 sigma, 含 xx/yy/xy 三个平面应力分量 D = mu * b / (2 * pi * (1 - nu)); r2 = x.^2 + y.^2; r4 = r2.^2; sigma.xx = -D * y .* (3 * x.^2 + y.^2) ./ r4; sigma.yy = D * y .* (x.^2 - y.^2) ./ r4; sigma.xy = D * x .* (x.^2 - y.^2) ./ r4; end

代码把 x、y 当作数组处理，用.*和./做元素级运算，所以调用一次可以同时算出滑移面上几百个离散点的应力快照。位错芯附近的 r2 趋于零，应力公式发散，实际使用要加一个截断半径，常见做法是取 rc = b，对小于 rc 的观测点强制赋有限值。这既避免 MATLAB 里出现 inf 或 NaN，也给位错对的短程排斥提供一个最简单的物理截断。

提示：截断半径不要取到滑移面长度的 0.1 倍以上，否则应力峰会被人为抹平，沿滑移面的应力分布曲线会丢失位错堆积的关键特征。

2.2 Peach-Koehler 力的展开：为什么滑移方向只看 σ_xy

位错能否沿滑移面滑动，取决于它受到的 Peach-Koehler 力。线力密度通用形式是 (\mathbf{F} = (\sigma\cdot\mathbf{b})\times\boldsymbol{\xi})，其中 (\boldsymbol{\xi}) 是位错线单位切向。二维模拟里 (\boldsymbol{\xi}) 取 ((0,0,\pm1))，把叉乘展开后，滑移面内的分力化简成两个分量。对 Burgers 矢量沿 x 方向的纯刃型位错 (\mathbf{b}=(b,0,0))，公式进一步简化为：

Fx = sigma_xy * b Fy = -sigma_xx * b

也就是说，位错沿滑移方向的驱动力完全由切应力 (\sigma_{xy}) 决定，正应力 (\sigma_{xx}) 则倾向于把位错推向相邻滑移面。forcePeachKoehler.m 就是把这段展开写成独立函数，我一般让它接收应力张量、Burgers 矢量和线方向三个参数，方便读入不同符号的位错。Burgers 矢量符号由 dislList.txt 里的数据决定，符号写反则受力方向整体反向，位错偶极子会从相互吸引变成相互排斥。

function [Fx, Fy] = forcePeachKoehler(sigma, bx, by, xi) % sigma: 位错核心位置处的 2x2 应力张量 % bx, by: Burgers 矢量分量, xi: 位错线切向, +1 或 -1 A_x = sigma(1,1) * bx + sigma(1,2) * by; A_y = sigma(2,1) * bx + sigma(2,2) * by; Fx = A_y * xi; Fy = -A_x * xi; end

注意xi的符号：同一根 Burgers 矢量如果 z 方向反向，Peach-Koehler 力也取反，这正是刃型位错偶极子两端符号相反、相互吸引的数学来源。createDislocation.m 和 createDislocationSource.m 生成位错对时会带符号字段，读列表时第一件事就是核对这一列。下面这张表把应力计算链上的几个函数串起来，方便定位问题：

2.3 多体位错应力叠加与自应力排除

DDD 中一个位错感受到的是其余所有位错应力场的叠加，dislocationStressField.m 每调用一次只贡献一根位错的应力，外层循环需要把位错数组完整遍历一遍。常见实现是排两根位错，一根作为源位错、一根作为受力位错，内层从 j = i+1 开始可以省掉一半重复计算。下面是应力累加和受力计算的核心逻辑：

for i = 1:nDisl sx = 0; sy = 0; sxy = 0; for j = 1:nDisl if i == j continue; % 跳过自应力，否则位错会被自己推走 end dx = disl(i).x - disl(j).x; dy = disl(i).y - disl(j).y; sig = dislocationStressField(dx, dy, disl(j).b, mu, nu); sx = sx + sig.xx; sy = sy + sig.yy; sxy = sxy + sig.xy; end sigmaAtI = [sx, sxy; sxy, sy]; [disl(i).Fx, disl(i).Fy] = forcePeachKoehler(sigmaAtI, ... disl(i).bx, disl(i).by, disl(i).xi); end

自应力这行continue是 DDD 计算里最常见的错误来源：如果不跳过，位错会受到自身奇异场的巨大伪力，模拟几步就会发散。sortDislocations.m 按 x 坐标先排序，再在相邻区间内配对，可以把严格的双重循环剪枝成近邻搜索，位错数量到几百根时计算量差距非常明显。这个项目目前只处理单个滑移面上的几十根位错，O(N²) 还能接受，但源码层面保留排序函数，显然是在为往多滑移面扩展做准备。

3. 输入文件与初始化：slipPlane.txt、dislList.txt、dsourceList.txt 的装配

3.1 slipPlane.txt 的几何定义与 readSlipPlane 解析

滑移面的几何在这个模型里就是一条线段：起点坐标、方向角、长度构成最基本的三要素。slipPlane.txt 里每一行对应一个滑移面，常见列格式是 x0 y0 angle length，角度单位用度还是弧度，要看 readSlipPlane.m 里是否乘了 pi/180，这一步很容易踩坑。我拿到代码的第一件事是在命令窗口打印 readSlipPlane('slipPlane.txt') 的返回结果，确认角度没有差 180 倍。

滑移面初始化之后，createSlipPlane.m 会按一定间距把线段离散成大量观测点，这些点就是后面计算滑移面应力分布的位置取样。离散间距直接影响应力曲线的光滑度：间距太大峰值被平均掉，间距太小则容易把数值噪声放大。合理的起点是把滑移面长度分成 500 到 1000 段，既能看清位错堆积的应力峰，又不至于让投影计算太慢。readSlipPlane 返回的结构体里，除了四个几何量，还应该带离散点数组供后面重复使用。

3.2 dislList.txt 的位错数组：符号、坐标和滑移约束

dislList.txt 描述的是初始时刻已经存在的位错。每一行至少包含位错编号、x 坐标、y 坐标、Burgers 矢量模长 b 和符号 sgn。符号表示位错沿滑移面的滑移方向，正号位错的 Burgers 矢量沿 +x，负号沿 -x。readDislocationList.m 读进来之后通常转成 MATLAB 结构体数组，字段就是 disl(i).x、disl(i).y、disl(i).b、disl(i).sgn。

位错必须被约束在滑移面上运动，所以每次更新位置后，常见做法是把新坐标投影回滑移线，这就是 projectVector.m 的用途之一。如果读入的位错坐标不落在滑移线附近，程序不会报错，但应力分布会出现奇怪的偏移。初始化时可以先做一次位置校验：计算每个位错到滑移线所在直线的距离，大于一个容忍值就打印警告。

固定位错也会出现在 dislList.txt 里，用符号或额外 flag 字段标出。固定位错不参与时间步更新，只参与应力场计算，这个区分在程序里要特别小心。许多看起来像"位错穿过了晶界"的反常模拟图，本质上是把固定位错也推进了动力学循环。

3.3 dsourceList.txt 与位错源的发射机制

位错源是 Frank-Read 源的二维简化。dsourceList.txt 里记录源的位置、临界切应力和发射间距，createDislocationSource.m 在每个时间步检查位错源处的局部切应力是否超过临界值，一旦超过就在源点两侧生成一对 Burgers 符号相反的位错偶极子，并给它们一个微小初始间距，模拟 Frank-Read 源不断弓出位错环的过程。

这里有一个物理上容易忽略的点：新发射位错的初始间距必须大于位错芯截断半径，否则偶极子内部相互作用力趋于无穷大，时间步进一推就飞出去。常见做法是让初始间距等于 2 到 5 倍的截断半径。读入 dsourceList 后可以先把源信息整理成下面的表核对一遍：

三个文件的单位必须一致：如果 slipPlane.txt 按微米写、dislList.txt 按纳米写，应力分布会差三个数量级。项目文档没明确单位时，先跑一次 drawSimulation 看几何比例，滑移面线段和位错点是否落在同一尺度内，比翻代码找单位换算更直接。

4. dd2d.m 主循环：粘性拖拽、时间增量与位错源释放

4.1 从 Peach-Koehler 力到位错速度的粘性拖拽模型

位错在晶体中运动时受到声子阻尼和电子阻尼的共同作用，宏观表现为一个随速度线性增长的阻力。因此主循环不直接解牛顿第二定律，而是让滑移速度与 Peach-Koehler 力成正比：v = F / B，B 是拖拽系数。dd2d.m 在计算完 forcePeachKoehler 后会立刻用这个关系算出速度向量，再交给时间增量部分决定走多远。

initializeSimulation.m 里的 sim.B 就是干这个用的。常见金属在室温下的 B 量级是 1e-5 到 1e-4 Pa·s，设置太大会让位错几乎不动，太小则容易数值振荡。这个参数和材料、温度都强相关，项目 README 里如果没有给默认值，先按这个量级试跑，再用第 5 章提到的 Weertman 解对照校准。

4.2 dislocationTimeIncrement 与 dislocationPairTimeIncrement 互补定步长

时间步长是这里最敏感的参数。固定 dt 不可取：位错速度较大时，一步可能跨过好几个位错间距，位错直接跑到滑移面外。dislocationTimeIncrement.m 的常见设计是统计当前所有位错速度的最大值，给定一个最大允许位移 dx_max，令 dt = dx_max / v_max：

% dislocationTimeIncrement 的一种常见实现 % disl: 位错数组, dxMax: 单个时间步内允许的最大位移 % 返回全局时间步长 dt vMax = 0; for i = 1:numel(disl) v = sqrt(disl(i).Fx^2 + disl(i).Fy^2) / sim.B; vMax = max(vMax, v); end if vMax > 0 dt = dxMax / vMax; else dt = inf; % 系统静止时不再推进，由外层逻辑终止 end

只靠这个公式不够，因为两个反向位错快速靠近时，单根位错的位移限制挡不住二者互相冲过对方。dislocationPairTimeIncrement.m 的作用就是额外检查相邻位错对的距离 d 和相对速度 v_rel，取 dt_pair = c * d / v_rel，c 一般取 0.1 到 0.5。最终步长取这两个结果的最小值，再对 dt 加一个下限防止除零。三个时间增量函数的职责可以从文件名区分：

4.3 dd2d.m 主循环的装配与位错源冷却逻辑

把前面所有模块接起来，主循环是一个四步往复的过程：算总应力场、算每个位错受力、定时间步长并更新位置、检查位错源是否发射。位错源释放时调用 createDislocationSource.m，传入该源位置的局部切应力和临界值，发射完成后再重新排序，排序在下一个循环里自然生效。这里有一个工程细节容易被忽略：位错源不能在同一个时间步里重复发射，常见做法是在源结构体上加冷却时间字段。

% dd2d.m 主循环的简化骨架 for istep = 1:sim.nsteps [Fx, Fy] = computeAllForces(disl, sim); % 内部函数，绕大循环调 forcePeachKoehler dtStep = dislocationTimeIncrement(disl, sim); % 按最大速度定 dt dtPair = dislocationPairTimeIncrement(disl, sim); dt = min(dtStep, dtPair); for i = 1:numel(disl) disl(i).x = disl(i).x + Fx(i) / sim.B * dt; disl(i).y = disl(i).y + Fy(i) / sim.B * dt; end disl = sortDislocations(disl); % 每次更新后重排 for s = 1:numel(sources) if sources(s).coolTime == 0 tauLocal = computeLocalShear(disl, sources(s), sim); if abs(tauLocal) > sources(s).tauCritical disl = createDislocationSource(disl, sources(s), sim); sources(s).coolTime = sources(s).coolInterval; end else sources(s).coolTime = sources(s).coolTime - 1; end end if mod(istep, sim.snapStep) == 0 drawSimulation(disl, sim, istep); end end

冷却时间按整数步计数是个实用技巧：新位错刚发射时，source 的临界应力检查被冷却封住，防止同一源在 dt 很小时连续吐出一排位错。dt 本身又受 dislocationPairTimeIncrement 约束，因此新产生的偶极子不会在一步内叠加出非物理的巨大力。代码里的 computeAllForces 和 computeLocalShear 是拆分后的辅助函数，原项目可能把这两段逻辑直接写在 dd2d.m 或对应文件名里，只要职责划分一致，后续改成 GPU 数组或 C-MEX 加速时改动面都很小。

5. 沿滑移面的应力提取、可视化和 Weertman 解对照

5.1 用 projectVector 把全局应力张量投影到滑移方向

滑移面不一定是全局坐标的 x 轴，所以沿滑移面的切应力要做投影。一套常见做法是：先取滑移面的切向 (\mathbf{t}=(\cos\theta,\sin\theta))，再和张量做双线性投影 (\tau(s)=\mathbf{t}\cdot\sigma\cdot\mathbf{n})，其中 n 是滑移面法向。展开后，(\tau(s)=\sigma_{xy}\cos2\theta+(\sigma_{yy}-\sigma_{xx})/2\cdot\sin2\theta)，只有当滑移面与全局 x 轴平行时，(\tau) 才等于 (\sigma_{xy})。

projectVector.m 就是做这个投影的最小函数，输入 2x2 应力张量和法向单位向量，输出标量。slipPlaneStressDistribution.m 更进一步：遍历滑移线上所有离散点，调用 dislocationStressField 累加所有位错贡献，再逐点投影，最终返回一条一维分布 sDist = [s; τ(s)]，其中 s 是距滑移面起点的弧长。伪代码框架如下：

% slipPlaneStressDistribution 的逐点投影框架 % slipPts: 2xN 离散点, disl: 位错数组, theta: 滑移面方向角 tauArr = zeros(1, N); for k = 1:N sigma2d = zeros(2, 2); for j = 1:numel(disl) dx = slipPts(1,k) - disl(j).x; dy = slipPts(2,k) - disl(j).y; sig = dislocationStressField(dx, dy, disl(j).b, sim.mu, sim.nu); sigma2d = sigma2d + [sig.xx sig.xy; sig.xy sig.yy]; end nVec = [-sin(theta); cos(theta)]; % 滑移面法向单位向量 tVec = [cos(theta); sin(theta)]; % 滑移面切向单位向量 tauArr(k) = tVec' * sigma2d * nVec; end sDist = [slipPts(1, :); tauArr];

投影里的方向角 theta 要和 dislList 里的 Burgers 矢量方向一致：theta 是滑移方向与 x 轴的夹角，如果 txt 里给的是法向角，读入时就要换算。投影结果中位错芯位置会出现应力尖峰，这些尖峰附近数值没有连续介质意义，因为位错芯本身就是奇异点，观察时应跳过。要区分位错堆积导致的应力集中和数值尖刺，可以看尖峰两侧是否平滑衰减，后者通常是截断半径或离散间距设置不合理。

5.2 输出文件 simulated 和 weertman：用解析解做交叉验证

项目输出里的 simulated.fig/png 和 weertman.fig/png，分别画的是数值模拟结果和 Weertman 解析解。Weertman 给出的是位错塞积群的闭式应力解，在远离塞积头部的位置渐近成立。把两条曲线放在同一坐标系里看趋势，比只比峰值更可靠。具体做法是让横坐标用塞积群长度 L 归一化，纵坐标用外加剪应力归一化，这样不同材料参数下也能叠加对比。

drawSimulation.m 负责画位错位置和滑移面几何，plotSlipPlaneStress.m 画滑移切应力曲线。MATLAB 里的 color map 不建议用 jet，会把应力尖峰的视觉权重放大，改用 parula 更客观。位错位置画成散点，滑移面画成直线，应力曲线用独立坐标轴叠加在同一张图上。如果模拟结果的应力峰位置和 Weertman 解的峰位错开超过 20%，优先怀疑 Burgers 矢量符号约定，而不是解析解本身。

5.3 排错：曲线失控时第一个检查点是什么

如果 plotSlipPlaneStress 的输出剧烈振荡，先从三个量查起：dislocationStressField 的 r2 是否被零除、slipPlane 的离散间距是否小于位错芯截断半径、时间步是否过大导致位错穿过了离散取样点。这三种病态的曲线形态完全不同，第一种是尖刺型振荡，第二种是整条曲线毛刺化，第三种是位错位置上方出现台阶跳变，看一眼就能区分。

另一种隐蔽问题是符号约定冲突。当 dislList.txt 的 sgn 列方向和 createDislocation.m 里的默认符号不一致时，Peach-Koehler 力反向，曲线会在滑移面两端出现对称的镜像峰。验证方法很直接：把 dislList.txt 的符号列全部取反再跑一次，若生成的曲线也镜像翻转，说明符号是唯一误差源，物理参数没有问题。

6. 往多滑移面扩展前，先动手做的三个检查

6.1 扫描位错源临界切应力，验证模型灵敏度

位错源临界切应力 tau_c 是模型里最直观的旋钮。把 dsourceList.txt 里的 tau_c 按 0.8、1.0、1.2 倍各跑一遍，对比滑移面应力曲线的最大应力和峰位。如果最大值随 tau_c 单调变化，说明位错源的"发射-松弛"机制整体自洽；如果曲线出现跳变而不是渐变，多半是新偶极子生成后与邻近位错瞬间强烈交互导致的步进不稳定，此时应调小 dislocationPairTimeIncrement 里的比例系数 c，而不是去改物理参数。

6.2 检查固定位错在排序和更新中的角色

把固定位错、可动位错都放进排序数组，但只在更新阶段跳过固定位错的坐标增量。一个容易漏掉的细节是：sortDislocations.m 排序后，固定位错的下标变了，如果源代码里用固定位错编号做索引，排序必须在更新前完成，且固定位错的坐标不允许出现在 drawSimulation 的增量箭头里。常用的解决办法是给固定位错加一个 flag 字段，在位移更新和受力累加两个入口同时做判断。

6.3 用 Weertman 解给后续扩展留一个回归基准

在修改任何参数之前，先把当前单滑移面的模拟曲线和 weertman.png 导出成一组 baseline 数据。之后每改动一次代码，就跑一遍同样输入，比较曲线的积分绝对误差而非峰值误差。这样在扩展成多滑移面时，如果应力分布出现趋势性偏移，能立刻区分是新增平面之间的交叉项算错，还是原有单面计算被改坏。扩展时在 slipPlane.txt 里按行添加不同角度的滑移面，记得给每个平面分配独立的位错符号约定和投影方向，避免多平面共用同一个全局 theta 导致应力分布张冠李戴。

本文还有配套的精品资源，点击获取