尧图网站建设 
简介面向材料计算与相场模拟学习者这份案例包聚焦相场模型在微观组织演变及相变模拟中的数值实现适用于材料科学、计算物理等方向的科研与工程实践。压缩包共含10个文件以MATLAB源码为主7个.m覆盖自由能计算、有限差分求解、晶粒形核与生长等核心模块并附2个.avi模拟演示视频和1个.inp输入文件整体仅1.01MB轻量便携。目前已有221人学习下载适合正在入门相场模型或面临建模困难的科研人员与工程师。通过解析源码与视频可直观理解相场变量的定义、自由能函数构建、动力学方程离散化及边界条件设置快速复现25晶粒演化案例进而迁移到凝固、腐蚀、裂纹扩展等复杂相变场景中。1. 一个 zip 里的相场模型从文件名到物理问题解压case_study_2.zip之后最先看到的是free_energ_fd_ca_v1.m、fd_ca_v1.m、grain_25.inp以及两个 avi 视频很多人会误以为这只是一个普通的作业附件。但把fd_ca拆开看fd是有限差分Finite Differenceca是 Cahn-Hilliard 方程的缩写这其实是一套完整的相场模型Phase Field Model模拟程序。它解决的问题很具体在一个二维计算区域内多个晶粒如何通过界面迁移和曲率驱动粗化最终形成更大尺寸的组织。相比直接追踪尖锐界面的方法相场模型用一个连续变量c描述相状态界面处c从 0 连续过渡到 1因此不需要显式标记界面位置天然可以处理晶粒合并、消失和拓扑变化。这套代码非常适合刚接触计算材料学、又想在 MATLAB 里快速跑通一个相场案例的科研人员也适合需要把模拟结果输出成 VTK 进行可视化后处理的工程师。视频fd_2g_1c.avi和fd_25g_1c.avi分别展示了 2 晶粒和 25 晶粒的演化从启动文件到最终画面能直观感受到界面能驱动的粗化过程。2. Cahn-Hilliard 方程与有限差分离散laplacian.m 和稳定性边界在拆代码之前必须先说清楚这个模型描述的是什么。Cahn-Hilliard 方程是一类守恒场动力学方程被广泛用于相分离、粗化和晶粒生长模拟。与 Allen-Cahn 方程的关键差异在于Cahn-Hilliard 要求整个计算域内浓度场c的总量守恒这一点可以从方程的形式直接看出∂c/∂t M ∇² μ其中M是迁移率μ是化学势定义为自由能泛函的变分导数μ δF/δc f(c) - κ ∇² c自由能泛函包含两项局部自由能密度f(c)和界面梯度能(κ/2)|∇c|²。经典的局部自由能密度通常是双阱势f(c) A c² (1 - c)²这样c0和c1是两个能量极小值对应两个相c0.5是能量最高点所以界面会被压缩成一个弥散带。κ控制界面宽度A控制能量势垒高度。这个模型之所以适合晶粒粗化是因为曲率越大、凸界面处的化学势越高物质会从高曲率区域向低曲率区域扩散导致小晶粒缩小、大晶粒长大。代码中对应的离散部分非常直接。laplacian.m实现了标准五点差分格式这是整个有限差分方案的基石。下面是一个符合该资源命名习惯的实现包含周期性边界条件function lap laplacian(field, dx, dy) % field: 二维标量场 % dx, dy: 网格间距 [ny, nx] size(field); lap zeros(ny, nx); % 内部网格点使用五点差分 lap(2:end-1, 2:end-1) ... (field(1:end-2, 2:end-1) - 2*field(2:end-1, 2:end-1) field(3:end, 2:end-1)) / dx^2 ... (field(2:end-1, 1:end-2) - 2*field(2:end-1, 2:end-1) field(2:end-1, 3:end)) / dy^2; % 周期性边界把边缘看成与对侧相邻 lap(1, 2:end-1) lap(1, 2:end-1) ... (field(end, 2:end-1) - 2*field(1, 2:end-1) field(2, 2:end-1)) / dy^2; lap(end, 2:end-1) lap(end, 2:end-1) ... (field(end-1, 2:end-1) - 2*field(end, 2:end-1) field(1, 2:end-1)) / dy^2; lap(2:end-1, 1) lap(2:end-1, 1) ... (field(2:end-1, end) - 2*field(2:end-1, 1) field(2:end-1, 2)) / dx^2; lap(2:end-1, end) lap(2:end-1, end) ... (field(2:end-1, end-1) - 2*field(2:end-1, end) field(2:end-1, 1)) / dx^2; end这个函数把二阶导数离散成了三点模板中间点权重为-2左右两个相邻点权重为1再除以网格间距的平方。边界处的处理是为了满足周期假设适合模拟一个大块材料内部晶粒演化的情形。如果换成实验中的有限样品表面则需要改为零通量边界对应 Neumann 条件下的差分处理。在显式时间推进下完整的更新格式是mu df_energy(c) - kappa * laplacian(c, dx, dy); dc M * laplacian(mu, dx, dy); c c dt * dc;这样做虽然简单但稳定性限制极其苛刻。Cahn-Hilliard 方程的显式格式要求时间步长满足Δt dx^4 / (4 M κ)也就是说如果你把网格间距缩小一半时间步长需要缩小到原来的十六分之一。这个约束让v1版本在较大网格下寸步难行。v2版本通常会对κ∇⁴c这一项做半隐式处理或者在傅里叶空间更新才能用更大的Δt完成长时间演化。下表列出了几个关键符号的含义后续调试参数时会反复用到符号含义典型取值范围c相场变量0 和 1 代表两相0.0 ~ 1.0M迁移率控制扩散速率0.1 ~ 10.0A双阱势高度控制势垒0.5 ~ 4.0κ梯度能系数控制界面宽度0.1 ~ 2.0dx, dy网格间距0.5 ~ 2.0Δt时间步长受稳定性约束1e-5 ~ 1e-3理解这些符号后再看fd_ca_v1.m和fd_ca_v2.m的差异就不会迷惑。v1往往是为了教学把每一步都写成显式便于检查v2则在保证守恒量的前提下牺牲一些实现简洁度换取更高的效率。后文拆解代码时我会以v1为主线指出v2的改进位置。3. MATLAB 代码模块拆解从 init_grain_micro.m 到 write_vtk_grid_values.m这个压缩包里的文件划分非常清晰init_grain_micro.m负责初始微结构生成free_energ_fd_ca_v1.m负责计算自由能导数laplacian.m负责空间二阶导数fd_ca_v1.m负责时间推进主循环write_vtk_grid_values.m负责输出可视化数据。先看整体文件功能对照表文件作用输入输出init_grain_micro.m初始化晶粒位置和浓度场网格尺寸、晶粒数初始c场free_energ_fd_ca_v1.m计算局部自由能导数df/dc浓度场c、参数A导数场laplacian.m计算二维拉普拉斯算子标量场、网格间距二阶导数场fd_ca_v1.m主程序读入参数并推进输入文件名或参数结构体时间序列、视频数据fd_ca_v2.m改进版可采用半隐式推进同上更稳定的演化结果write_vtk_grid_values.m写出 VTK 结构化数据场变量、网格信息、文件名.vti或.vtk文件init_grain_micro.m的常见做法是先随机生成晶粒种子点再把每个像素归属到最近的种子点最后通过一个平滑函数把边界变成弥散界面。下面这段代码模拟了这个思路function c init_grain_micro(nx, ny, n_grains, interface_width) % 生成 n_grains 个随机种子点 rng(42); % 固定随机种子保证结果可复现 seeds rand(n_grains, 2); seeds(:, 1) seeds(:, 1) * nx; seeds(:, 2) seeds(:, 2) * ny; % 计算每个网格点到最近种子点的欧氏距离 [X, Y] meshgrid(1:nx, 1:ny); dist zeros(ny, nx); for i 1:n_grains d sqrt((X - seeds(i, 1)).^2 (Y - seeds(i, 2)).^2); if i 1 dist d; else dist min(dist, d); end end % 用双曲正切函数生成弥散界面界面厚度由 interface_width 控制 c 0.5 * (1 - tanh(dist / interface_width)); end这里的逻辑是距离种子点越近c越接近 0距离越远c越接近 1。所以晶粒内部和外部分别对应两个不同的相晶界处c在约三倍interface_width范围内逐渐变化。固定随机种子非常重要否则每次初始化都不同后面对比参数影响时会引入额外噪声。interface_width一般取 2~4 个网格间距太大会让初始界面过厚无法分辨单个晶粒太小则容易在后续演化中产生数值振荡。接下来是自由能导数部分。free_energ_fd_ca_v1.m的核心是双阱势的导数一次性返回f(c)而不是返回自由能本身。原因是化学势计算中只需要导数function df free_energ_fd_ca_v1(c, A) % c: 0 到 1 之间的场变量 % A: 双阱势系数 df 2 * A * c .* (1 - c) .* (1 - 2*c); end这个式子的来源是对f(c) A c²(1-c)²求导。可以看到c0和c1处导数为 0对应两个稳定相c0.5处导数也为 0但这是不稳定平衡点。在数值模拟中如果初始场中含有一点噪声那么浓度场会在双阱势的驱动下自动分离成两相。主程序fd_ca_v1.m把这些模块组织在一起。核心循环大约是这样% 读入参数初始化场 c init_grain_micro(params.nx, params.ny, params.n_grains, params.interface_width); M params.M; A params.A; kappa params.kappa; dx params.dx; dy params.dy; dt params.dt; n_steps params.n_steps; save_every params.save_every; for step 1:n_steps % 化学势 f(c) - kappa*laplacian(c) mu free_energ_fd_ca_v1(c, A) - kappa * laplacian(c, dx, dy); % 浓度通量 M * laplacian(mu) dc M * laplacian(mu, dx, dy); % 显式欧拉时间推进 c c dt * dc; % 数值修正确保 c 始终在允许范围内 c max(0, min(1, c)); % 周期性保存可视化数据 if mod(step, save_every) 0 filename sprintf(output_%05d.vtk, step); write_vtk_grid_values(filename, c, dx, dy); end end注意这里的c max(0, min(1, c))只是兜底手段如果频繁触发说明时间步长过大或界面宽度与网格尺度不匹配应该回头调整参数而不是依赖压制。显式欧拉法每步只需要计算两次拉普拉斯一次用于化学势一次用于通量所以成本不高但稳定性局限明显。v2版本通常会把c从实空间变换到傅里叶空间用(i k)²代替拉普拉斯算子然后对线性项做隐式推进。这样Δt可以放大十倍到上百倍但代价是代码复杂度上升。write_vtk_grid_values.m是输出部分它把二维场写成 VTK 格式方便在 ParaView 中查看。标准 VTK 结构化点数据格式可以用下面这种方式生成function write_vtk_grid_values(filename, field, dx, dy) [ny, nx] size(field); fid fopen(filename, w); fprintf(fid, # vtk DataFile Version 3.0\n); fprintf(fid, phase field output\n); fprintf(fid, ASCII\n); fprintf(fid, DATASET STRUCTURED_POINTS\n); fprintf(fid, DIMENSIONS %d %d 1\n, nx, ny); fprintf(fid, SPACING %f %f 1.0\n, dx, dy); fprintf(fid, ORIGIN 0 0 0\n); fprintf(fid, POINT_DATA %d\n, nx * ny); fprintf(fid, SCALARS concentration double 1\n); fprintf(fid, LOOKUP_TABLE default\n); % 按行优先输出与 MATLAB 的列优先不同 fprintf(fid, %f\n, transpose(field(:))); fclose(fid); end这段代码的关键在于维度顺序和transpose处理。VTK 的STRUCTURED_POINTS按 x 方向最快变化MATLAB 是列优先存储所以直接输出field(:)会得到 y 方向优先的数据导致图像转置。用transpose转成行优先后再输出才能得到正确朝向。如果你在 ParaView 里看到的界面方向反了先检查这里而不是去怀疑物理模型。4. grain_25.inp 实战从 2 晶粒到 25 晶粒的参数与运行在主程序跑起来之前需要从grain_25.inp读入参数。这个文件从命名上可以看出是为 25 晶粒准备的输入文件而对应的结果视频是fd_25g_1c.avi。inp文件不是 MATLAB 原生格式所以需要自己写一个轻量解析函数。下面是一个典型的参数解析实现function params read_inp(filename) fid fopen(filename, r); % 先读全部行再按关键字解析 raw textscan(fid, %s %s, CommentStyle, #); fclose(fid); keys raw{1}; values raw{2}; params struct(); for i 1:length(keys) switch lower(keys{i}) case nx, params.nx str2double(values{i}); case ny, params.ny str2double(values{i}); case dx, params.dx str2double(values{i}); case dy, params.dy str2double(values{i}); case M, params.M str2double(values{i}); case A, params.A str2double(values{i}); case kappa, params.kappa str2double(values{i}); case dt, params.dt str2double(values{i}); case n_steps, params.n_steps str2double(values{i}); case n_grains, params.n_grains str2double(values{i}); case save_every, params.save_every str2double(values{i}); case interface_width, params.interface_width str2double(values{i}); otherwise, warning(未知参数: %s, keys{i}); end end end对应的grain_25.inp内容大致如下nx 256 ny 256 dx 1.0 dy 1.0 M 1.0 A 2.0 kappa 1.0 dt 1e-4 n_steps 20000 n_grains 25 save_every 100 interface_width 2.5上面save_every表示每 100 步写一个 VTKn_steps为总步数。这样在 20000 步后会生成 200 个帧足以看清楚粗化过程。如果你想和fd_2g_1c.avi做对照把n_grains改成 2其他参数不动再跑一遍就能看到同样的初始条件、不同的晶粒密度对演化行为的影响。运行方式在 MATLAB 中很简单切换到解压目录然后在命令行执行params read_inp(grain_25.inp); fd_ca_v1(params);如果fd_ca_v1.m的入口函数没有接受结构体参数也可以把它当作脚本直接修改文件末尾的硬编码参数。但我的经验是尽量用输入文件驱动仿真这样每次实验的参数都能留档避免改完代码改不回原样。v1版本可以打印一些关键量来观察演化状态比如在每个保存步计算总浓度sum(c(:))如果sum相对初始值漂移超过 1%就说明数值误差过大需要调小Δt或改用半隐式格式。我在跑 25 晶粒的时候最初用dt1e-3结果几步后c直接出现棋盘式振荡最终变成了类似椒盐噪声的图案。这是因为显式格式不满足Δt dx^4/(4Mκ)而dx1, M1, κ1下临界值大约是0.25但实际求解中交叉项会进一步压缩稳定区间所以1e-3在前面几步可以不崩时间长了必然发散。换成1e-4之后视频里的界面迁移变得平滑可以看到典型的长程粗化小晶粒逐渐缩小大晶粒缓慢长大最终到达一个准静态状态。参数调优方面重点关注kappa与interface_width的关系。理论上 Cahn-Hilliard 的弥散界面宽度正比于sqrt(kappa / A)。如果你把kappa调大、A调小界面会变宽反之界面会变锐利。但界面不能比网格间距小太多否则离散误差会让界面像锯齿一样移动。我的建议是让界面宽度保持在 3~5 个网格点上也就是sqrt(kappa / A) / dx ≈ 2 ~ 4例如dx1.0, A2.0时kappa取0.5到2.0都是安全范围。下表给出了一些测试组合的实际效果n_grainskappaAdt界面宽度格点现象21.02.01e-4约 3单个界面缓慢平直化251.02.01e-4约 3多晶粗化晶粒数量逐渐减少250.52.05e-5约 2界面更薄但波动性增加252.02.02e-4约 4界面更宽稳定但细节分辨率差如果在演化的中后期发现晶粒形状变成奇怪的多边形而不是平滑的弧线通常不是物理问题而是界面宽度只占了一个网格点各向异性来自差分模板的方向偏好。这时把kappa调大一点或者改用九点差分格式可以显著减轻网格各向异性。5. 进阶技巧能量校验、界面宽度标定与谱方法加速比起盯着视频看界面是否平滑数值验证更重要。相场模型有两个硬性约束质量守恒和能量单调递减。前者可以通过sum(c(:))在任何时刻检查后者需要计算系统总自由能。总自由能包括局部自由能和梯度能两部分离散形式如下function [total_energy, bulk_energy, gradient_energy] phase_field_energy(c, dt, dx, dy, A, kappa) bulk A * c.^2 .* (1 - c).^2; [gx, gy] gradient(c, dx, dy); grad_energy 0.5 * kappa * (gx.^2 gy.^2); total_energy sum(bulk(:)) sum(grad_energy(:)); end如果total_energy在某个时间段上升尤其是初期急剧上升几乎可以肯定是数值不稳定导致浓度场越界。正确的做法是先用一个小dt跑几十步确认能量曲线单调下降再逐步放大dt。不要相信“只跑一步”的结果因为显式格式可以在第一步保持守恒但后续振荡会悄悄累积。界面宽度标定也是一个实用技巧。在初始化时可以通过测量c0.25到c0.75之间的距离得到界面厚度。MATLAB 里可以用contourc提取等值线c_contour contourc(c, [0.25 0.75]);然后计算这两条等值线之间的平均距离。均匀界面的条件下理论界面厚度约为2.44 * sqrt(kappa / A)如果测量结果远大于理论值说明网格分辨率不足远小于理论值说明界面处采样点太少。这个测量结果可以作为后续网格加密或kappa调整的依据。如果想让模拟跑得更快一个直接的办法是放弃纯有限差分改用傅里叶谱方法。Cahn-Hilliard 方程中的拉普拉斯算子在谱空间是对角的两步更新可以合并成一次快速傅里叶变换。这里给出一个半隐式的谱方法实现片段% kx, ky 为频率向量k2 kx.^2 ky.^2 c_hat fft2(c); mu_hat fft2(free_energ_fd_ca_v1(c, A)) kappa * k2 .* c_hat; c_hat_new (c_hat - dt * M * k2 .* mu_hat) ./ (1 dt * M * kappa * k2.^2); c real(ifft2(c_hat_new));这里将-κ∇²c的线性项隐式化把f(c)的非线性项留在显式部分从而把稳定性约束从Δt ~ dx^4放宽到Δt ~ dx²。实际使用时可以比显式格式大几十倍代价是代码不够直观对边界条件也没有周期假设之外的选择。如果在fd_ca_v2.m中看到类似的结构说明这套案例已经切到更高效率的实现上了。用上述能量校验函数配合谱方法可以在几分钟内完成原本需要数小时的大规模晶粒粗化模拟。本文还有配套的精品资源点击获取 Cahn-Hilliard相场模型与有限差分：MATLAB晶粒粗化模拟实战
简介面向材料计算与相场模拟学习者这份案例包聚焦相场模型在微观组织演变及相变模拟中的数值实现适用于材料科学、计算物理等方向的科研与工程实践。压缩包共含10个文件以MATLAB源码为主7个.m覆盖自由能计算、有限差分求解、晶粒形核与生长等核心模块并附2个.avi模拟演示视频和1个.inp输入文件整体仅1.01MB轻量便携。目前已有221人学习下载适合正在入门相场模型或面临建模困难的科研人员与工程师。通过解析源码与视频可直观理解相场变量的定义、自由能函数构建、动力学方程离散化及边界条件设置快速复现25晶粒演化案例进而迁移到凝固、腐蚀、裂纹扩展等复杂相变场景中。1. 一个 zip 里的相场模型从文件名到物理问题解压case_study_2.zip之后最先看到的是free_energ_fd_ca_v1.m、fd_ca_v1.m、grain_25.inp以及两个 avi 视频很多人会误以为这只是一个普通的作业附件。但把fd_ca拆开看fd是有限差分Finite Differenceca是 Cahn-Hilliard 方程的缩写这其实是一套完整的相场模型Phase Field Model模拟程序。它解决的问题很具体在一个二维计算区域内多个晶粒如何通过界面迁移和曲率驱动粗化最终形成更大尺寸的组织。相比直接追踪尖锐界面的方法相场模型用一个连续变量c描述相状态界面处c从 0 连续过渡到 1因此不需要显式标记界面位置天然可以处理晶粒合并、消失和拓扑变化。这套代码非常适合刚接触计算材料学、又想在 MATLAB 里快速跑通一个相场案例的科研人员也适合需要把模拟结果输出成 VTK 进行可视化后处理的工程师。视频fd_2g_1c.avi和fd_25g_1c.avi分别展示了 2 晶粒和 25 晶粒的演化从启动文件到最终画面能直观感受到界面能驱动的粗化过程。2. Cahn-Hilliard 方程与有限差分离散laplacian.m 和稳定性边界在拆代码之前必须先说清楚这个模型描述的是什么。Cahn-Hilliard 方程是一类守恒场动力学方程被广泛用于相分离、粗化和晶粒生长模拟。与 Allen-Cahn 方程的关键差异在于Cahn-Hilliard 要求整个计算域内浓度场c的总量守恒这一点可以从方程的形式直接看出∂c/∂t M ∇² μ其中M是迁移率μ是化学势定义为自由能泛函的变分导数μ δF/δc f(c) - κ ∇² c自由能泛函包含两项局部自由能密度f(c)和界面梯度能(κ/2)|∇c|²。经典的局部自由能密度通常是双阱势f(c) A c² (1 - c)²这样c0和c1是两个能量极小值对应两个相c0.5是能量最高点所以界面会被压缩成一个弥散带。κ控制界面宽度A控制能量势垒高度。这个模型之所以适合晶粒粗化是因为曲率越大、凸界面处的化学势越高物质会从高曲率区域向低曲率区域扩散导致小晶粒缩小、大晶粒长大。代码中对应的离散部分非常直接。laplacian.m实现了标准五点差分格式这是整个有限差分方案的基石。下面是一个符合该资源命名习惯的实现包含周期性边界条件function lap laplacian(field, dx, dy) % field: 二维标量场 % dx, dy: 网格间距 [ny, nx] size(field); lap zeros(ny, nx); % 内部网格点使用五点差分 lap(2:end-1, 2:end-1) ... (field(1:end-2, 2:end-1) - 2*field(2:end-1, 2:end-1) field(3:end, 2:end-1)) / dx^2 ... (field(2:end-1, 1:end-2) - 2*field(2:end-1, 2:end-1) field(2:end-1, 3:end)) / dy^2; % 周期性边界把边缘看成与对侧相邻 lap(1, 2:end-1) lap(1, 2:end-1) ... (field(end, 2:end-1) - 2*field(1, 2:end-1) field(2, 2:end-1)) / dy^2; lap(end, 2:end-1) lap(end, 2:end-1) ... (field(end-1, 2:end-1) - 2*field(end, 2:end-1) field(1, 2:end-1)) / dy^2; lap(2:end-1, 1) lap(2:end-1, 1) ... (field(2:end-1, end) - 2*field(2:end-1, 1) field(2:end-1, 2)) / dx^2; lap(2:end-1, end) lap(2:end-1, end) ... (field(2:end-1, end-1) - 2*field(2:end-1, end) field(2:end-1, 1)) / dx^2; end这个函数把二阶导数离散成了三点模板中间点权重为-2左右两个相邻点权重为1再除以网格间距的平方。边界处的处理是为了满足周期假设适合模拟一个大块材料内部晶粒演化的情形。如果换成实验中的有限样品表面则需要改为零通量边界对应 Neumann 条件下的差分处理。在显式时间推进下完整的更新格式是mu df_energy(c) - kappa * laplacian(c, dx, dy); dc M * laplacian(mu, dx, dy); c c dt * dc;这样做虽然简单但稳定性限制极其苛刻。Cahn-Hilliard 方程的显式格式要求时间步长满足Δt dx^4 / (4 M κ)也就是说如果你把网格间距缩小一半时间步长需要缩小到原来的十六分之一。这个约束让v1版本在较大网格下寸步难行。v2版本通常会对κ∇⁴c这一项做半隐式处理或者在傅里叶空间更新才能用更大的Δt完成长时间演化。下表列出了几个关键符号的含义后续调试参数时会反复用到符号含义典型取值范围c相场变量0 和 1 代表两相0.0 ~ 1.0M迁移率控制扩散速率0.1 ~ 10.0A双阱势高度控制势垒0.5 ~ 4.0κ梯度能系数控制界面宽度0.1 ~ 2.0dx, dy网格间距0.5 ~ 2.0Δt时间步长受稳定性约束1e-5 ~ 1e-3理解这些符号后再看fd_ca_v1.m和fd_ca_v2.m的差异就不会迷惑。v1往往是为了教学把每一步都写成显式便于检查v2则在保证守恒量的前提下牺牲一些实现简洁度换取更高的效率。后文拆解代码时我会以v1为主线指出v2的改进位置。3. MATLAB 代码模块拆解从 init_grain_micro.m 到 write_vtk_grid_values.m这个压缩包里的文件划分非常清晰init_grain_micro.m负责初始微结构生成free_energ_fd_ca_v1.m负责计算自由能导数laplacian.m负责空间二阶导数fd_ca_v1.m负责时间推进主循环write_vtk_grid_values.m负责输出可视化数据。先看整体文件功能对照表文件作用输入输出init_grain_micro.m初始化晶粒位置和浓度场网格尺寸、晶粒数初始c场free_energ_fd_ca_v1.m计算局部自由能导数df/dc浓度场c、参数A导数场laplacian.m计算二维拉普拉斯算子标量场、网格间距二阶导数场fd_ca_v1.m主程序读入参数并推进输入文件名或参数结构体时间序列、视频数据fd_ca_v2.m改进版可采用半隐式推进同上更稳定的演化结果write_vtk_grid_values.m写出 VTK 结构化数据场变量、网格信息、文件名.vti或.vtk文件init_grain_micro.m的常见做法是先随机生成晶粒种子点再把每个像素归属到最近的种子点最后通过一个平滑函数把边界变成弥散界面。下面这段代码模拟了这个思路function c init_grain_micro(nx, ny, n_grains, interface_width) % 生成 n_grains 个随机种子点 rng(42); % 固定随机种子保证结果可复现 seeds rand(n_grains, 2); seeds(:, 1) seeds(:, 1) * nx; seeds(:, 2) seeds(:, 2) * ny; % 计算每个网格点到最近种子点的欧氏距离 [X, Y] meshgrid(1:nx, 1:ny); dist zeros(ny, nx); for i 1:n_grains d sqrt((X - seeds(i, 1)).^2 (Y - seeds(i, 2)).^2); if i 1 dist d; else dist min(dist, d); end end % 用双曲正切函数生成弥散界面界面厚度由 interface_width 控制 c 0.5 * (1 - tanh(dist / interface_width)); end这里的逻辑是距离种子点越近c越接近 0距离越远c越接近 1。所以晶粒内部和外部分别对应两个不同的相晶界处c在约三倍interface_width范围内逐渐变化。固定随机种子非常重要否则每次初始化都不同后面对比参数影响时会引入额外噪声。interface_width一般取 2~4 个网格间距太大会让初始界面过厚无法分辨单个晶粒太小则容易在后续演化中产生数值振荡。接下来是自由能导数部分。free_energ_fd_ca_v1.m的核心是双阱势的导数一次性返回f(c)而不是返回自由能本身。原因是化学势计算中只需要导数function df free_energ_fd_ca_v1(c, A) % c: 0 到 1 之间的场变量 % A: 双阱势系数 df 2 * A * c .* (1 - c) .* (1 - 2*c); end这个式子的来源是对f(c) A c²(1-c)²求导。可以看到c0和c1处导数为 0对应两个稳定相c0.5处导数也为 0但这是不稳定平衡点。在数值模拟中如果初始场中含有一点噪声那么浓度场会在双阱势的驱动下自动分离成两相。主程序fd_ca_v1.m把这些模块组织在一起。核心循环大约是这样% 读入参数初始化场 c init_grain_micro(params.nx, params.ny, params.n_grains, params.interface_width); M params.M; A params.A; kappa params.kappa; dx params.dx; dy params.dy; dt params.dt; n_steps params.n_steps; save_every params.save_every; for step 1:n_steps % 化学势 f(c) - kappa*laplacian(c) mu free_energ_fd_ca_v1(c, A) - kappa * laplacian(c, dx, dy); % 浓度通量 M * laplacian(mu) dc M * laplacian(mu, dx, dy); % 显式欧拉时间推进 c c dt * dc; % 数值修正确保 c 始终在允许范围内 c max(0, min(1, c)); % 周期性保存可视化数据 if mod(step, save_every) 0 filename sprintf(output_%05d.vtk, step); write_vtk_grid_values(filename, c, dx, dy); end end注意这里的c max(0, min(1, c))只是兜底手段如果频繁触发说明时间步长过大或界面宽度与网格尺度不匹配应该回头调整参数而不是依赖压制。显式欧拉法每步只需要计算两次拉普拉斯一次用于化学势一次用于通量所以成本不高但稳定性局限明显。v2版本通常会把c从实空间变换到傅里叶空间用(i k)²代替拉普拉斯算子然后对线性项做隐式推进。这样Δt可以放大十倍到上百倍但代价是代码复杂度上升。write_vtk_grid_values.m是输出部分它把二维场写成 VTK 格式方便在 ParaView 中查看。标准 VTK 结构化点数据格式可以用下面这种方式生成function write_vtk_grid_values(filename, field, dx, dy) [ny, nx] size(field); fid fopen(filename, w); fprintf(fid, # vtk DataFile Version 3.0\n); fprintf(fid, phase field output\n); fprintf(fid, ASCII\n); fprintf(fid, DATASET STRUCTURED_POINTS\n); fprintf(fid, DIMENSIONS %d %d 1\n, nx, ny); fprintf(fid, SPACING %f %f 1.0\n, dx, dy); fprintf(fid, ORIGIN 0 0 0\n); fprintf(fid, POINT_DATA %d\n, nx * ny); fprintf(fid, SCALARS concentration double 1\n); fprintf(fid, LOOKUP_TABLE default\n); % 按行优先输出与 MATLAB 的列优先不同 fprintf(fid, %f\n, transpose(field(:))); fclose(fid); end这段代码的关键在于维度顺序和transpose处理。VTK 的STRUCTURED_POINTS按 x 方向最快变化MATLAB 是列优先存储所以直接输出field(:)会得到 y 方向优先的数据导致图像转置。用transpose转成行优先后再输出才能得到正确朝向。如果你在 ParaView 里看到的界面方向反了先检查这里而不是去怀疑物理模型。4. grain_25.inp 实战从 2 晶粒到 25 晶粒的参数与运行在主程序跑起来之前需要从grain_25.inp读入参数。这个文件从命名上可以看出是为 25 晶粒准备的输入文件而对应的结果视频是fd_25g_1c.avi。inp文件不是 MATLAB 原生格式所以需要自己写一个轻量解析函数。下面是一个典型的参数解析实现function params read_inp(filename) fid fopen(filename, r); % 先读全部行再按关键字解析 raw textscan(fid, %s %s, CommentStyle, #); fclose(fid); keys raw{1}; values raw{2}; params struct(); for i 1:length(keys) switch lower(keys{i}) case nx, params.nx str2double(values{i}); case ny, params.ny str2double(values{i}); case dx, params.dx str2double(values{i}); case dy, params.dy str2double(values{i}); case M, params.M str2double(values{i}); case A, params.A str2double(values{i}); case kappa, params.kappa str2double(values{i}); case dt, params.dt str2double(values{i}); case n_steps, params.n_steps str2double(values{i}); case n_grains, params.n_grains str2double(values{i}); case save_every, params.save_every str2double(values{i}); case interface_width, params.interface_width str2double(values{i}); otherwise, warning(未知参数: %s, keys{i}); end end end对应的grain_25.inp内容大致如下nx 256 ny 256 dx 1.0 dy 1.0 M 1.0 A 2.0 kappa 1.0 dt 1e-4 n_steps 20000 n_grains 25 save_every 100 interface_width 2.5上面save_every表示每 100 步写一个 VTKn_steps为总步数。这样在 20000 步后会生成 200 个帧足以看清楚粗化过程。如果你想和fd_2g_1c.avi做对照把n_grains改成 2其他参数不动再跑一遍就能看到同样的初始条件、不同的晶粒密度对演化行为的影响。运行方式在 MATLAB 中很简单切换到解压目录然后在命令行执行params read_inp(grain_25.inp); fd_ca_v1(params);如果fd_ca_v1.m的入口函数没有接受结构体参数也可以把它当作脚本直接修改文件末尾的硬编码参数。但我的经验是尽量用输入文件驱动仿真这样每次实验的参数都能留档避免改完代码改不回原样。v1版本可以打印一些关键量来观察演化状态比如在每个保存步计算总浓度sum(c(:))如果sum相对初始值漂移超过 1%就说明数值误差过大需要调小Δt或改用半隐式格式。我在跑 25 晶粒的时候最初用dt1e-3结果几步后c直接出现棋盘式振荡最终变成了类似椒盐噪声的图案。这是因为显式格式不满足Δt dx^4/(4Mκ)而dx1, M1, κ1下临界值大约是0.25但实际求解中交叉项会进一步压缩稳定区间所以1e-3在前面几步可以不崩时间长了必然发散。换成1e-4之后视频里的界面迁移变得平滑可以看到典型的长程粗化小晶粒逐渐缩小大晶粒缓慢长大最终到达一个准静态状态。参数调优方面重点关注kappa与interface_width的关系。理论上 Cahn-Hilliard 的弥散界面宽度正比于sqrt(kappa / A)。如果你把kappa调大、A调小界面会变宽反之界面会变锐利。但界面不能比网格间距小太多否则离散误差会让界面像锯齿一样移动。我的建议是让界面宽度保持在 3~5 个网格点上也就是sqrt(kappa / A) / dx ≈ 2 ~ 4例如dx1.0, A2.0时kappa取0.5到2.0都是安全范围。下表给出了一些测试组合的实际效果n_grainskappaAdt界面宽度格点现象21.02.01e-4约 3单个界面缓慢平直化251.02.01e-4约 3多晶粗化晶粒数量逐渐减少250.52.05e-5约 2界面更薄但波动性增加252.02.02e-4约 4界面更宽稳定但细节分辨率差如果在演化的中后期发现晶粒形状变成奇怪的多边形而不是平滑的弧线通常不是物理问题而是界面宽度只占了一个网格点各向异性来自差分模板的方向偏好。这时把kappa调大一点或者改用九点差分格式可以显著减轻网格各向异性。5. 进阶技巧能量校验、界面宽度标定与谱方法加速比起盯着视频看界面是否平滑数值验证更重要。相场模型有两个硬性约束质量守恒和能量单调递减。前者可以通过sum(c(:))在任何时刻检查后者需要计算系统总自由能。总自由能包括局部自由能和梯度能两部分离散形式如下function [total_energy, bulk_energy, gradient_energy] phase_field_energy(c, dt, dx, dy, A, kappa) bulk A * c.^2 .* (1 - c).^2; [gx, gy] gradient(c, dx, dy); grad_energy 0.5 * kappa * (gx.^2 gy.^2); total_energy sum(bulk(:)) sum(grad_energy(:)); end如果total_energy在某个时间段上升尤其是初期急剧上升几乎可以肯定是数值不稳定导致浓度场越界。正确的做法是先用一个小dt跑几十步确认能量曲线单调下降再逐步放大dt。不要相信“只跑一步”的结果因为显式格式可以在第一步保持守恒但后续振荡会悄悄累积。界面宽度标定也是一个实用技巧。在初始化时可以通过测量c0.25到c0.75之间的距离得到界面厚度。MATLAB 里可以用contourc提取等值线c_contour contourc(c, [0.25 0.75]);然后计算这两条等值线之间的平均距离。均匀界面的条件下理论界面厚度约为2.44 * sqrt(kappa / A)如果测量结果远大于理论值说明网格分辨率不足远小于理论值说明界面处采样点太少。这个测量结果可以作为后续网格加密或kappa调整的依据。如果想让模拟跑得更快一个直接的办法是放弃纯有限差分改用傅里叶谱方法。Cahn-Hilliard 方程中的拉普拉斯算子在谱空间是对角的两步更新可以合并成一次快速傅里叶变换。这里给出一个半隐式的谱方法实现片段% kx, ky 为频率向量k2 kx.^2 ky.^2 c_hat fft2(c); mu_hat fft2(free_energ_fd_ca_v1(c, A)) kappa * k2 .* c_hat; c_hat_new (c_hat - dt * M * k2 .* mu_hat) ./ (1 dt * M * kappa * k2.^2); c real(ifft2(c_hat_new));这里将-κ∇²c的线性项隐式化把f(c)的非线性项留在显式部分从而把稳定性约束从Δt ~ dx^4放宽到Δt ~ dx²。实际使用时可以比显式格式大几十倍代价是代码不够直观对边界条件也没有周期假设之外的选择。如果在fd_ca_v2.m中看到类似的结构说明这套案例已经切到更高效率的实现上了。用上述能量校验函数配合谱方法可以在几分钟内完成原本需要数小时的大规模晶粒粗化模拟。本文还有配套的精品资源点击获取