简介面向金属塑性变形与位错动力学研究者的二维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 DynamicsDDD和有限元、分子动力学都不一样它把材料看作一根根可滑移的位错线集合通过应力场驱动位错运动。这里的 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 生成位错对时会带符号字段读列表时第一件事就是核对这一列。下面这张表把应力计算链上的几个函数串起来方便定位问题函数文件计算目标关键输入输出dislocationStressField.m单根位错在观测点的应力张量x, y, b, mu, nusigma(xx, yy, xy)forcePeachKoehler.m位错受力sigma, bx, by, xiFx, FyprojectVector.m应力向滑移方向投影sigma, 方向向量标量应力sortDislocations.m位错排序加速累加位错结构数组排序后的数组2.3 多体位错应力叠加与自应力排除DDD 中一个位错感受到的是其余所有位错应力场的叠加dislocationStressField.m 每调用一次只贡献一根位错的应力外层循环需要把位错数组完整遍历一遍。常见实现是排两根位错一根作为源位错、一根作为受力位错内层从 j i1 开始可以省掉一半重复计算。下面是应力累加和受力计算的核心逻辑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.txtx0 y0 angle length滑移面起点、方向角、长度角度制转换dislList.txtid x y b sgn初始位错几何与 Burgers 符号符号错误则受力反向dsourceList.txtid x y tau_c span位错源位置、临界应力和发射间距初始间距要大于 rc三个文件的单位必须一致如果 slipPlane.txt 按微米写、dislList.txt 按纳米写应力分布会差三个数量级。项目文档没明确单位时先跑一次 drawSimulation 看几何比例滑移面线段和位错点是否落在同一尺度内比翻代码找单位换算更直接。4. dd2d.m 主循环粘性拖拽、时间增量与位错源释放4.1 从 Peach-Koehler 力到位错速度的粘性拖拽模型位错在晶体中运动时受到声子阻尼和电子阻尼的共同作用宏观表现为一个随速度线性增长的阻力。因此主循环不直接解牛顿第二定律而是让滑移速度与 Peach-Koehler 力成正比v F / BB 是拖拽系数。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_relc 一般取 0.1 到 0.5。最终步长取这两个结果的最小值再对 dt 加一个下限防止除零。三个时间增量函数的职责可以从文件名区分函数文件判断依据解决的问题timeIncrement.m全局基准步长模拟启动阶段的默认值dislocationTimeIncrement.m单根位错最大速度限制每个位错的滑移距离dislocationPairTimeIncrement.m位错对间距与相对速度防止偶极子对穿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 导致应力分布张冠李戴。本文还有配套的精品资源点击获取