相场模拟——枝晶生长（COMSOL模拟雪花形成）
一、介绍
1.1 物理含义
雪花是人们最常见的枝晶。枝晶生长是一种生长的不稳定现象，常起因于过冷的液体，或晶体的生长速度受限于溶质原子向固体表面的扩散速度。造成以上条件的原因，可以是液相中负的温度梯度，也可以是负的浓度梯度。在这种模式下，晶体倾向于在其棱角处优先生长，从而形成突触状结构。
这篇博文会介绍用相场模拟的方法，模拟雪花生长的过程。
1.2 相场模拟介绍
在相场法中，使用连续变量描述微观结构特征。这些变量有两种形式：代表物理性质的守恒变量，如原子浓度或材料密度；描述材料微观结构(包括晶粒和不同相)的非守恒序参数(order parameters)。这些连续变量的演化用自由能的函数表达，可以定义为一个PDE偏微分方程系统。该相场模型还可以与其他物理相耦合，例如力学或热传导。这些模型方程可以用多种方法求解，包括有限差分法、谱方法和有限元法。
相场PDE是各种变量的演化方程，是自由能泛函F的函数。
对于守恒变量（conserved variables），其满足Cahn-Hilliard方程：
对于非保守量，其满足 Allen-Cahn方程：
自由能函数包括局部自由能，梯度自由能，以及系统的其它额外自由能（比如电势能，弹性势能等等）。可以用如下公式计算：
二、枝晶生长模型
模型包含两个变量: 相场θ(r, t)，温度场T(r, t)
θ = 0 表示液态；θ = 1 表示固态；
PDE：
其它方程：
三、COMSOL建模
1. 建立一个2D，两个变量的General Form PDE模型(自变量为theta, T)，瞬态求解。
2. Global Definition-> Parameters
K 1.8 1.8
TAU 0.0003 3E-4
EPS 0.01 0.01
DELTA 0.02 0.02
ANGLE0 1.57 1.57
ANISO 6 6
ALPHA 0.9 0.9
GAMMA 10.0 10
TEQ 1.0 1
4. Model -> Definition -> Variables
epsilon EPS*(1+DELTA*cos(ANISO*(angle-ANGLE0)))
epsilonp -EPS*ANISO*DELTA*sin(ANISO*(angle-ANGLE0))
epsilon2 epsilon*epsilon
m ALPHA/pi*atan(GAMMA*(TEQ-T)) rad
angle if(abs(thetax)>eps,atan2(thetay,thetax)+pi*(atan2(thetay,thetax)<0),pi/2*((thetay>0)+3*(thetay<0))) rad
5. Geometry
建立一个关于原点对称的9×9矩形
6. PDE设置
初始值
7. Mesh
建立Mapped型网格，Maximum element size为0.1(可适当加密，0.06~0.08左右效果比较好)
8. Study求解时间为range(0,0.1,0.3)
导出的部分matlab代码：
function out = model
%
% crystalGrowth.m
%
% Model exported on Jul 16 2023, 16:41 by COMSOL 6.0.0.318.
import com.comsol.model.*
import com.comsol.model.util.*
model = ModelUtil.create('Model');
model.component.create('comp1', true);
model.component('comp1').geom.create('geom1', 2);
model.component('comp1').mesh.create('mesh1');
model.component('comp1').physics.create('g', 'GeneralFormPDE', {'theta' 'T'});
model.study.create('std1');
model.study('std1').create('time', 'Transient');
model.study('std1').feature('time').setSolveFor('/physics/g', true);
model.param.set('K', '1.8');
model.param.set('TAU', '0.0003');
model.param.set('EPS', '0.01');
model.param.set('DELTA', '0.02');
model.param.set('ANGLE0', '1.57');
model.param.set('ANISO', '6');
model.param.set('ALPHA', '0.9');
model.param.set('GAMMA', '10.0');
model.param.set('TEQ', '1.0');
model.func.create('step1', 'Step');
model.func('step1').set('location', 0.09);
model.func('step1').set('from', 1);
model.func('step1').set('to', 0);
model.func('step1').set('smooth', 0.05);
model.variable.create('var1');
model.variable.remove('var1');
model.component('comp1').variable.create('var1');
model.component('comp1').geom('geom1').run;
model.component('comp1').variable('var1').set('epsilon', 'EPS*(1+DELTA*cos(ANISO*(angle-ANGLE0)))');
model.component('comp1').variable('var1').set('epsilonp', '-EPS*ANISO*DELTA*sin(ANISO*(angle-ANGLE0))');
model.component('comp1').variable('var1').set('epsilon2', 'epsilon*epsilon');
model.component('comp1').variable('var1').set('m', 'ALPHA/pi*atan(GAMMA*(TEQ-T))');
model.component('comp1').variable('var1').set('angle', 'if(abs(thetax)>eps,atan2(thetay,thetax)+pi*(atan2(thetay,thetax)<0),pi/2*((thetay>0)+3*(thetay<0)))');
model.component('comp1').geom('geom1').create('r1', 'Rectangle');
model.component('comp1').geom('geom1').feature('r1').set('size', [9 9]);
model.component('comp1').geom('geom1').runPre('fin');
model.component('comp1').geom('geom1').run;
model.component('comp1').physics('g').feature('gfeq1').setIndex('Ga', {'-thetax*epsilon2+epsilon*epsilonp*thetay' '-thetay'}, 0);
model.component('comp1').physics('g').feature('gfeq1').setIndex('Ga', {'-thetax*epsilon2+epsilon*epsilonp*thetay' '-thetay*epsilon2-epsilon*epsilonp*thetax'}, 0);
model.component('comp1').physics('g').feature('gfeq1').setIndex('f', 'd(epsilon2,x)*thetax+d(epsilon2,y)*thetay+theta*(1-theta)*(theta-0.5+m)', 0);
model.component('comp1').physics('g').feature('gfeq1').setIndex('f', 0, 1);
model.component('comp1').physics('g').feature('gfeq1').setIndex('da', 'TAU', 0);
model.component('comp1').physics('g').feature('gfeq1').setIndex('da', '-K', 1);
model.component('comp1').physics('g').feature('init1').set('theta', 'step1(x^2+y^2)');
model.component('comp1').mesh('mesh1').create('size1', 'Size');
model.component('comp1').mesh('mesh1').feature.remove('size1');
model.component('comp1').mesh('mesh1').create('map1', 'Map');
model.component('comp1').mesh('mesh1').feature('size').set('hmax', 0.1);
model.component('comp1').mesh('mesh1').feature('size').set('hmin', '6.75e-4');
model.component('comp1').mesh('mesh1').feature('size').set('hgrad', 1.2);
model.component('comp1').mesh('mesh1').feature('size').set('hcurve', 0.25);
model.component('comp1').mesh('mesh1').run;
model.sol.create('sol1');
model.sol('sol1').study('std1');
model.study('std1').feature('time').set('notlistsolnum', 1);
model.study('std1').feature('time').set('notsolnum', 'auto');
model.study('std1').feature('time').set('listsolnum', 1);
model.study('std1').feature('time').set('solnum', 'auto');
model.sol('sol1').create('st1', 'StudyStep');
model.sol('sol1').feature('st1').set('study', 'std1');
model.sol('sol1').feature('st1').set('studystep', 'time');
model.sol('sol1').create('v1', 'Variables');
model.sol('sol1').feature('v1').set('control', 'time');
model.sol('sol1').create('t1', 'Time');
model.sol('sol1').feature('t1').set('tlist', 'range(0,0.1,1)');
model.sol('sol1').feature('t1').set('plot', 'off');
model.sol('sol1').feature('t1').set('plotgroup', 'Default');
model.sol('sol1').feature('t1').set('plotfreq', 'tout');
model.sol('sol1').feature('t1').set('probesel', 'all');
model.sol('sol1').feature('t1').set('probes', {});
model.sol('sol1').feature('t1').set('probefreq', 'tsteps');
model.sol('sol1').feature('t1').set('atolglobalvaluemethod', 'factor');
model.sol('sol1').feature('t1').set('endtimeinterpolation', true);
model.sol('sol1').feature('t1').set('control', 'time');
model.sol('sol1').feature('t1').create('fc1', 'FullyCoupled');
model.sol('sol1').feature('t1').feature('fc1').set('linsolver', 'dDef');
model.sol('sol1').feature('t1').feature.remove('fcDef');
model.sol('sol1').attach('std1');
model.result.create('pg1', 'PlotGroup2D');
model.result('pg1').set('data', 'dset1');
model.result('pg1').create('surf1', 'Surface');
model.result('pg1').feature('surf1').set('expr', 'theta');
model.sol('sol1').runAll;
四、结果
参考：
https://mooseframework.inl.gov/modules/phase_field/Phase_Field_Equations.html
https://mooseframework.inl.gov/modules/phase_field/Derivation_explanations.html
https://blog.sina.com.cn/s/blog_4a0a8b5d01011spl.html