从Cahn-Hilliard方程到相场模拟：旋节分解的数值实现与工程应用

简介：本资源是一套基于MATLAB实现的旋节线分解（Spinodal Decomposition）数值模拟工具，面向材料科学、物理化学及计算力学领域的研究生、科研人员与高年级本科生，用于理解并可视化多组分系统在无核化条件下的自发相分离行为。资源包共5个文件，含2个核心MATLAB脚本（主程序与拉普拉斯算子实现）、1份MIT开源许可证、1个演示视频（mp4）及1份说明文档（md），整体仅1.01MB，轻量易部署。已有146人学习下载，适合快速上手Cahn-Hilliard方程建模与求解。用户可直接运行脚本复现浓度场演化过程，结合视频直观掌握初始扰动设置、自由能驱动机制与微结构形貌演变规律，并借助代码结构清晰的模块设计（如独立laplacian.m）开展参数敏感性分析或算法改进。

1. 项目概述：从“Spinodal Decomposition”说起

最近在材料模拟和相图计算领域，一个名为“nsbalbi-Spinodal-Decomposition-v1.0-0-g7fbea5e”的项目引起了我的注意。这个项目标题看起来像是一个Git仓库的提交记录，核心指向一个经典的材料科学现象：Spinodal Decomposition，中文常译为“旋节分解”或“调幅分解”。如果你从事材料科学、物理冶金、计算材料学，或者对微观组织演化模拟感兴趣，那么这个项目及其背后的原理，绝对值得你花时间深入了解。它不是一个简单的脚本合集，而是理解合金、玻璃、高分子共混物等复杂材料在特定条件下如何自发分相、形成纳米结构的一把钥匙。

简单来说，Spinodal Decomposition描述了一种特殊的相分离过程。想象一下，一杯被剧烈摇晃的油水混合物，静置后它们会逐渐分层，这是普通的成核-生长机制。而旋节分解则不同，它发生在热力学绝对不稳定的区域内，整个体系没有任何能量壁垒，像一座本身就不稳固的山坡，任何微小的成分起伏都会自发地、连续地放大，最终导致均匀的单相分解为成分周期性调制的两相结构。这个过程没有明显的“成核”阶段，分解是全域同时发生的。这个项目，很可能就是通过数值方法（如Cahn-Hilliard方程）来模拟这一迷人的物理过程。

它能做什么？对于研究者，你可以用它来可视化不同初始条件、热力学参数下旋节分解的动力学过程，预测最终的组织形貌（比如是互联结构还是颗粒状结构）。对于工程师，理解旋节分解有助于设计具有特定性能的材料，例如高强度高韧性的调幅分解强化合金，或者具有特殊光学、电学性能的纳米复合材料。对于学生和爱好者，这是一个绝佳的学习工具，将抽象的相图理论和偏微分方程与直观的动态图像联系起来。接下来，我将带你深入这个项目的核心，拆解其背后的技术逻辑、实操要点，并分享我在复现和拓展此类模拟时的经验与教训。

2. 核心原理与模型深度解析

要真正玩转这个“Spinodal Decomposition”模拟项目，不能只停留在运行代码的层面。我们必须吃透其背后的物理模型和数学框架，这样才能在调整参数、解释结果时心中有数，而不是盲目试错。

2.1 热力学基础：自由能与相图

一切始于吉布斯自由能。对于一个二元体系（比如A-B合金），在某一温度T下，其自由能G通常是成分c（例如B组元的原子分数）的函数，曲线形状至关重要。当自由能曲线是“上凸”的（即二阶导数 ∂²G/∂c² > 0），体系是稳定的，微小起伏会被抑制。而当曲线出现“下凹”区域（∂²G/∂c² < 0）时，体系对于无限小的成分起伏是不稳定的，这个区域就是旋节区。

在相图上，旋节区位于两相平衡的化学势公切线所包围的“山脊”之内，它比传统的两相区范围要小。项目模拟的初始条件，通常会将体系的平均成分设置在这个旋节区内。理解这一点是关键：如果你的初始平均成分不在旋节区，那么模拟可能不会发生典型的旋节分解，而是需要克服能垒的成核过程，这需要用不同的模型（如相场法结合成核理论）来描述。

2.2 Cahn-Hilliard方程：动力学核心

描述旋节分解时空演化的灵魂方程是Cahn-Hilliard方程。它是一个四阶的非线性偏微分方程。别被吓到，我们可以把它拆解开来理解：

∂c/∂t = M ∇² (δF/δc)

其中：

• c(r, t)是位置r和时间t处的成分。
• M是原子迁移率，与扩散系数相关，通常假设为常数。
• ∇²是拉普拉斯算子（散度的梯度）。
• δF/δc是总自由能F对成分c的变分导数，可以理解为“化学势”。

总自由能F通常包含两部分：体自由能和梯度能。F = ∫_V [f(c) + κ|∇c|²] dV

1. 体自由能 f(c)：就是前面提到的均匀体系的自由能密度。在模拟中，为了简化，常采用双阱势，例如f(c) = - (A/2) c² + (B/4) c⁴。这里的A和B是正的温度相关参数。这个函数在c=0和c=1处有两个极小值（两相），在c=0附近有一个极大值，完美描述了旋节区的不稳定性。
2. 梯度能项 κ|∇c|²：这是Cahn和Hilliard的关键贡献。它惩罚成分在空间上的剧烈变化，反映了相界面具有额外能量的事实。κ是梯度能系数，为正数。正是这一项，使得分解形成的两相之间不是锐利的界面，而是有一定宽度的扩散界面，并且限制了分解的波长不会无限小。

将F的表达式代入Cahn-Hilliard方程，经过运算，我们得到常用的形式：∂c/∂t = M ∇² [f'(c) - 2κ ∇² c]

这里f'(c)是体自由能对c的一阶导数。方程右边括号内的项就是广义的化学势差。这个方程清晰地表明：物质流动（导致成分变化）正比于化学势梯度的散度，而化学势又由体自由能的驱动力和梯度能的抑制作用共同决定。

2.3 数值离散化：从连续到网格

计算机无法处理连续的方程，我们必须将其离散化。项目通常采用有限差分法在规则的二维或三维网格上进行。

• 空间离散：将模拟区域划分为Nx×Ny的网格，每个格点(i, j)赋予一个成分值c_{i,j}。拉普拉斯算子∇²可以用中心差分格式近似，例如在二维情况下：(∇² c)_{i,j} ≈ (c_{i+1,j} + c_{i-1,j} + c_{i,j+1} + c_{i,j-1} - 4c_{i,j}) / (Δx)²其中Δx是网格间距。注意，对于四阶项∇²(∇² c)，需要应用两次拉普拉斯算子。

• 时间离散：采用显式欧拉法最简单，但稳定性要求时间步长Δt非常小，计算效率低。c_{new} = c_{old} + Δt * RHS_{old}其中RHS是Cahn-Hilliard方程的右端项。更稳健的方法是采用半隐式或全隐式格式，例如对线性项（-2κM∇⁴ c）做隐式处理，可以允许更大的Δt。这个项目中可能采用了傅里叶谱方法，这是处理此类周期性边界条件和线性高阶导数的利器。

注意：选择显式还是隐式格式，是精度与计算成本之间的权衡。对于学习和小规模二维模拟，显式法简单直观；但对于大规模三维模拟或追求长期稳定性，半隐式或谱方法是更专业的选择。你需要检查项目代码中的时间迭代部分来确认。

2.4 初始条件与边界条件

• 初始条件：为了触发旋节分解，初始成分场通常设置为均匀平均成分c0加上一个微小的随机扰动，例如c(x,y,0) = c0 + η * rand(-0.01, 0.01)。c0必须位于旋节区内（例如对于对称双阱势，c0=0.5）。扰动η要足够小，以确保分解是由热力学不稳定性驱动，而非大的初始起伏。
• 边界条件：最常用的是周期性边界条件。这意味着模拟区域的左边界和右边界相连，上边界和下边界相连，像一个环面。这消除了复杂表面效应的影响，专注于体相分解行为，并且特别适合使用傅里叶谱方法求解。

3. 代码结构与关键模块拆解

假设“nsbalbi-Spinodal-Decomposition-v1.0”项目结构相对清晰，我们可以推断其核心模块构成。一个典型的旋节分解模拟代码会包含以下几个部分：

3.1 参数定义与初始化模块

这是脚本的开头部分，所有“旋钮”都在这里设置。你需要重点关注以下参数：

# 模拟参数 Nx, Ny = 256, 256 # 网格数。越大分辨率越高，计算越慢。256x256是平衡点。 dx, dy = 1.0, 1.0 # 网格间距（无量纲）。通常设为1。 dt = 0.1 # 时间步长。**这是关键易错点！** 必须满足稳定性条件。 nsteps = 10000 # 总模拟步数。决定模拟多久。 # 物理参数 c0 = 0.5 # 平均成分。0.5对应对称双阱势的旋节区中心。 A = 1.0 # 双阱势参数A，控制不稳定性深度。 B = 1.0 # 双阱势参数B，控制两相平衡成分。 kappa = 0.5 # 梯度能系数κ。影响界面宽度和特征波长。 M = 1.0 # 迁移率M。影响分解动力学速度。 # 初始化成分场 c = np.ones((Nx, Ny)) * c0 c += 0.01 * (np.random.rand(Nx, Ny) - 0.5) # 添加微小随机扰动

实操心得：dt的选择至关重要。对于显式格式，稳定性要求dt < (dx^4) / (常数 * M * kappa)，这个常数与离散格式有关。一个经验法则是先设一个非常小的dt（如0.01）确保运行，然后逐步增大，观察结果是否出现数值发散（成分值爆炸或变成NaN）。如果采用谱方法，稳定性限制会宽松很多。

3.2 核心计算循环模块

这是模拟的心脏，一个巨大的时间循环。每一步都包含：

1. 计算体自由能导数f_prime = -A*c + B*c**3（对于双阱势f(c) = -A/2 c^2 + B/4 c^4）。
2. 计算化学势mu = f_prime - 2*kappa*laplacian(c)。这里需要实现一个拉普拉斯算子函数。
3. 计算化学势的拉普拉斯laplacian_mu。
4. 更新成分场c_new = c_old + dt * M * laplacian_mu。

如果项目使用了傅里叶谱方法，代码会看起来截然不同。它会将成分场c变换到傅里叶空间（使用FFT），在那里拉普拉斯算子简单地变为乘以-k^2（k是波矢），从而高效处理高阶导数。更新方程在傅里叶空间中可能是一个标量乘法。

# 伪代码示例（谱方法思路） for step in range(nsteps): # 将c变换到傅里叶空间 c_hat = np.fft.fft2(c) # 计算波数网格 kx, ky # 在傅里叶空间计算更新：c_hat_new = c_hat * (某个与k相关的因子) + ... # 逆变换回实空间得到新的c c = np.fft.ifft2(c_hat_new).real # 可能还需要处理数值误差导致c超出[0,1]范围的情况

3.3 可视化与输出模块

模拟的价值在于观察。这个模块负责定期（比如每100步）将成分场c以图像或视频的形式输出。

• 静态快照：使用matplotlib.pyplot.imshow(c, cmap='seismic', vmin=0, vmax=1)可以生成伪彩色图，清晰显示两相区域。seismic色图在显示正负波动时非常直观。
• 动态演化：将每一步的快照保存到列表，最后用matplotlib.animation或imageio库生成视频。这能让你直观看到分解的早期、中期和后期阶段。
• 定量分析：除了看图，更深入的分析包括：
  • 计算结构因子 S(k, t)：对成分场进行傅里叶变换并取模平方|FFT(c)|^2。它可以揭示主导波长及其随时间的变化。在旋节分解早期，S(k)会在某个特征波数k_max处出现峰值，且峰值随时间增长。
  • 计算相分数：统计c > c0的格点比例，可以跟踪两相的比例变化。
  • 计算界面面积/总能量：监控体系总自由能F随时间下降的过程，验证模拟的物理正确性。

3.4 工具函数模块

一些独立的、可复用的函数会被放在这里，例如：

• laplacian_2d(arr)：计算二维数组的离散拉普拉斯。
• free_energy_density(c, A, B)：计算体自由能密度。
• calc_total_free_energy(c, kappa, dx)：计算整个体系的总自由能（包含梯度能）。

一个结构良好的项目会将这些函数模块化，方便调试和复用。

4. 完整复现与模拟实操指南

现在，让我们抛开对现有项目的依赖，从头开始构建一个基础的二维旋节分解模拟器。这个过程能让你彻底掌握每一个环节。

4.1 环境准备与依赖安装

我们使用Python，因为它有强大的科学计算和可视化库。

# 创建并激活虚拟环境（推荐） python -m venv spinodal_env source spinodal_env/bin/activate # Linux/Mac # spinodal_env\Scripts\activate # Windows # 安装核心库 pip install numpy matplotlib scipy imageio # 如果追求更快的FFT，可以安装pyfftw，但numpy.fft对于入门已足够。

4.2 基础版显式有限差分模拟实现

以下是基于显式有限差分法的完整代码框架，包含详细注释：

import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation import imageio def laplacian_2d(c, dx=1.0): """计算二维数组的5点中心差分拉普拉斯。使用周期性边界条件。""" c_top = np.roll(c, shift=1, axis=0) c_bottom = np.roll(c, shift=-1, axis=0) c_left = np.roll(c, shift=1, axis=1) c_right = np.roll(c, shift=-1, axis=1) laplacian = (c_top + c_bottom + c_left + c_right - 4 * c) / (dx**2) return laplacian def update_cahn_hilliard_explicit(c, dt, M, A, B, kappa, dx): """显式更新一步Cahn-Hilliard方程。""" # 1. 计算体自由能导数 (对于 f(c) = -A/2 c^2 + B/4 c^4) f_prime = -A * c + B * c**3 # 2. 计算化学势 mu = f' - 2*kappa*laplacian(c) mu = f_prime - 2 * kappa * laplacian_2d(c, dx) # 3. 计算化学势的拉普拉斯 laplacian_mu = laplacian_2d(mu, dx) # 4. 更新成分场 c_new = c + dt * M * laplacian_mu # 5. (可选) 简单的数值稳定性处理，防止成分过度偏离 # c_new = np.clip(c_new, 0.0, 1.0) return c_new # 主模拟参数 Nx, Ny = 128, 128 # 初始可以用小网格测试 dx = 1.0 dt = 0.01 # 必须很小！对于显式格式，这是关键。 M = 1.0 A = 1.0 B = 1.0 kappa = 0.5 nsteps = 5000 save_interval = 100 # 每100步保存一帧 # 初始化 c0 = 0.5 c = c0 + 0.01 * (np.random.rand(Nx, Ny) - 0.5) # 存储快照用于生成动画 snapshots = [c.copy()] # 主模拟循环 print("开始模拟...") for step in range(1, nsteps+1): c = update_cahn_hilliard_explicit(c, dt, M, A, B, kappa, dx) if step % save_interval == 0: snapshots.append(c.copy()) print(f"已完成第 {step}/{nsteps} 步") print("模拟完成！")

4.3 可视化与结果分析

模拟完成后，我们生成动画并分析结构因子。

# 1. 生成动画 fig, ax = plt.subplots() im = ax.imshow(snapshots[0], cmap='seismic', vmin=0, vmax=1, animated=True) ax.set_title(f"Spinodal Decomposition (t={0})") plt.colorbar(im, ax=ax) def update_frame(frame): im.set_array(snapshots[frame]) ax.set_title(f"Spinodal Decomposition (t={frame*save_interval})") return [im] ani = FuncAnimation(fig, update_frame, frames=len(snapshots), interval=50, blit=True) ani.save('spinodal_decomposition.gif', writer='pillow', fps=10) plt.close() # 2. 计算并绘制最终的结构因子（功率谱） final_c = snapshots[-1] # 减去均值，关注涨落 c_fluctuation = final_c - np.mean(final_c) # 二维傅里叶变换 c_hat = np.fft.fft2(c_fluctuation) power_spectrum = np.abs(np.fft.fftshift(c_hat))**2 # 将零频移到中心 # 绘制结构因子 kx = np.fft.fftshift(np.fft.fftfreq(Nx, dx)) * 2 * np.pi ky = np.fft.fftshift(np.fft.fftfreq(Ny, dx)) * 2 * np.pi KX, KY = np.meshgrid(kx, ky) fig2, ax2 = plt.subplots(1, 2, figsize=(12, 5)) ax2[0].imshow(final_c, cmap='seismic', vmin=0, vmax=1) ax2[0].set_title("Final Composition Field") ax2[0].set_xlabel("x") ax2[0].set_ylabel("y") # 用对数坐标显示功率谱，更清晰 im_ps = ax2[1].imshow(np.log10(power_spectrum + 1e-10), extent=[kx.min(), kx.max(), ky.min(), ky.max()], cmap='viridis') ax2[1].set_title("Structure Factor S(k) (log scale)") ax2[1].set_xlabel("Wave number kx") ax2[1].set_ylabel("Wave number ky") plt.colorbar(im_ps, ax=ax2[1]) plt.tight_layout() plt.savefig('final_state_and_structure_factor.png', dpi=150) plt.show()

运行这段代码，你将看到一个动态的分解过程：初始均匀的灰色区域（成分0.5）迅速出现蓝红相间的斑点（代表两相），这些斑点逐渐粗化（Coarsening），界面变得清晰。结构因子图会显示一个圆环状的亮环，其半径的倒数对应着分解的特征波长。

5. 参数影响分析与高级话题探讨

仅仅能运行模拟是不够的，我们需要理解每个物理参数如何影响最终结果，并探索更复杂的场景。

5.1 关键参数的作用与影响

通过系统性地改变参数并重新模拟，我们可以总结出以下规律：

5.2 从旋节分解到粗化（Ostwald Ripening）

模拟运行足够长时间后，你会观察到在相分离完成后，小颗粒溶解、大颗粒长大的现象，这就是粗化。粗化过程由界面曲率驱动，旨在降低总的界面能。描述粗化动力学的经典定律是Lifshitz-Slyozov-Wagner (LSW)理论，它预测平均颗粒半径R随时间t的增长满足R^3 ∝ t。

你可以在模拟后期，通过图像分析技术（如阈值分割、连通域标记）来统计“颗粒”的平均尺寸，并绘制其随时间变化的对数图，验证是否与LSW理论的预测相符。这是将模拟与经典理论对照的绝佳练习。

5.3 引入弹性应变能

在实际材料中，新相与母相之间往往存在晶格失配，从而产生弹性应变能。这会显著改变分解形貌。例如，在立方各向异性弹性条件下，析出相可能沿着特定的软弹性方向排列，形成棋盘状或条带状结构，而不是各向同性的斑点。

在模型中引入弹性应变能，需要在总自由能F中加入一项弹性应变能密度(1/2) C_{ijkl} ε_{ij}^{el} ε_{kl}^{el}，其中ε^{el}是弹性应变，与成分场引起的本征应变（ε^0 ∝ (c - c0)）相关。这会使控制方程变得更加复杂（耦合了力学平衡方程），通常需要采用相场法的框架，并使用微弹性理论（Khachaturyan）来高效计算长程弹性相互作用。

经验分享：从纯扩散控制的Cahn-Hilliard模型升级到包含弹性效应的相场模型，是一个质的飞跃。计算量会大幅增加，因为需要求解额外的力学平衡方程（通常是线性弹性问题）。这时，傅里叶谱方法的优势更加凸显，因为弹性格林函数在傅里叶空间中有简洁的解析形式。如果你要做这方面的研究，准备好学习更多的连续介质力学和数值分析知识。

5.4 三维模拟的挑战与技巧

二维模拟直观且计算快，但真实世界是三维的。三维旋节分解会产生复杂的互联海绵状结构，其粗化动力学也与二维有所不同。

进行三维模拟时，挑战主要在于：

1. 计算量：网格点从N²变为N³，内存和计算时间激增。N=256的三维模拟网格点数是二维的256倍！
2. 可视化：三维标量场的可视化比二维复杂。可以使用等值面渲染（如mayavi,pyvista库）来展示两相界面。

优化建议：

• 使用谱方法，其在三维中的优势比有限差分法更明显。
• 利用并行计算。傅里叶变换（通过pyfftw或scipy.fft）和许多线性代数操作可以很好地并行化。
• 从较小的三维网格（如64³）开始，并考虑使用自适应网格细化（AMR）技术，只在界面附近使用细网格。
• 输出数据时，可以考虑输出切片（二维截面）进行快速检查，完整三维等值面渲染可以每隔很多步做一次。

6. 常见问题、调试技巧与性能优化

在实际操作中，你一定会遇到各种问题。下面是我踩过的一些坑和总结的解决方案。

6.1 数值不稳定（发散）

这是最常见的问题，表现为成分值迅速变得极大或变成NaN。

• 症状：图像出现彩色斑点（值超出正常范围），或程序报错“NaN encountered”。
• 原因与排查：
  1. 时间步长dt太大：这是首要怀疑对象。对于显式有限差分，稳定性条件极其苛刻。解决方案：将dt减小为原来的1/10、1/100再试。一个经验公式是dt < (dx^4) / (32 * M * kappa)，可以作为起点。
  2. 物理参数组合导致“刚度”：当M、A、κ的乘积很大时，方程右端项极大，即使很小的dt也会导致更新步长过大。解决方案：尝试对参数进行无量纲化。通常将长度单位设为界面宽度，时间单位设为扩散时间，从而将M和κ归一化到1附近的数量级。
  3. 初始扰动太大：如果初始随机扰动的幅度η过大，可能导致一开始就进入高度非线性区，引发不稳定。解决方案：确保η很小（如0.001 * c0）。
• 根本解决之道：放弃显式格式，实现半隐式格式或使用傅里叶谱方法。半隐式格式对线性高阶导数项进行隐式处理，允许dt提高几个数量级。谱方法则因其指数收敛性和对周期性问题的天然适配性，成为此类模拟的工业标准。

6.2 结果不物理或未发生分解

• 症状：模拟结束后，成分场几乎没变化，或者变化模式很奇怪（如出现条纹而非斑点）。
• 排查清单：
  1. 检查平均成分c0：确认它是否真的位于旋节区内。计算f''(c0) = -A + 3B*c0^2，如果大于0，则体系稳定，不会发生旋节分解。
  2. 检查梯度能系数κ：κ过大可能会过度抑制所有波长模式的增长，导致分解被“冻结”。尝试减小κ。
  3. 检查模拟时间是否足够：分解，尤其是后期粗化，是一个非常缓慢的过程。你可能需要增加nsteps（总步数）。可以通过监控总自由能F是否在持续下降来判断过程是否在继续。
  4. 检查边界条件：确保你的拉普拉斯算子实现正确使用了周期性边界条件（例如用np.roll）。错误的边界条件会引入虚假的反射或约束。
  5. 可视化中间过程：不要只看最终结果。保存并查看每一步或每若干步的快照，观察是否在早期有细微的结构出现后又消失。

6.3 性能瓶颈与优化

当网格变大或模拟步数增多时，速度会成为问题。

• 瓶颈分析：使用cProfile或line_profiler工具找出代码中最耗时的函数。在显式有限差分法中，通常是计算拉普拉斯算子的循环；在谱方法中，是FFT计算。
• 优化策略：
  1. 向量化操作：确保使用NumPy的数组整体运算，绝对避免在Python层写for循环遍历网格。我们的laplacian_2d函数使用np.roll就是向量化的典范。
  2. 升级算法：如前所述，傅里叶谱方法是解决此类周期性边界问题的最快方法。将显式差分升级为谱方法，可能带来数十倍甚至百倍的性能提升，并且稳定性更好。
  3. 使用更快的FFT库：用pyfftw替代numpy.fft，它可以调用更优化的FFTW库。
  4. GPU加速：对于超大规模三维模拟，可以考虑使用cupy或jax库在GPU上运行。FFT和数组运算在GPU上并行效率极高。
  5. 减少I/O和可视化开销：不要每一步都保存数据或绘图。可以每隔成百上千步保存一次，或者只在最后生成动画时从内存中读取缓存的数据。

6.4 复现他人工作的注意事项

当你拿到像“nsbalbi-Spinodal-Decomposition-v1.0”这样的项目代码时，想复现其结果：

1. 仔细阅读README和注释：了解作者使用的参数、算法和单位制。
2. 检查参数的无量纲化：这是最大的坑。作者可能使用了一套无量纲化方案，使得代码中的A、B、κ、M、dt、dx与物理值对应关系不明。尝试在论文或文档中寻找无量纲化公式。
3. 运行基准测试：如果作者提供了某个参数下的标准结果（如图像），先用完全相同的参数运行，对比是否一致。
4. 理解随机种子：初始随机扰动会影响分解的具体图案。如果作者固定了随机种子（如np.random.seed(42)），你也应该固定，才能得到完全一致的微观结构演化路径。如果不固定，宏观统计规律（如结构因子）应该一致，但具体图案会不同。
5. 版本依赖：检查所需的Python库及其版本。不同版本的库（特别是numpy的FFT默认值）可能导致细微的差异。

通过这个从原理到实践，从基础到进阶的完整梳理，相信你已经对“Spinodal Decomposition”模拟有了立体而深入的理解。这个项目就像一座桥梁，连接了抽象的热力学理论与生动的微观组织图像。动手去调整参数，观察现象，验证理论，你收获的将不仅仅是代码运行的结果，更是对材料相变这一复杂过程直观而深刻的洞察力。

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