Cahn-Hillard Equation

initialization (base)

FDLaplace[g_]:=Table[(
If[ix==1 || ix==gn || iy==1 || iy==gn,
g[[If[ix>1,ix-1,gn],iy]]+g[[If[ix<gn,ix+1,1],iy]]+
g[[ix,If[iy>1,iy-1,gn]]]+g[[ix,If[iy<gn,iy+1,1]]]
,
g[[ix-1,iy]]+g[[ix+1,iy]]+
g[[ix,iy-1]]+g[[ix,iy+1]]] 
-4*g[[ix,iy]])*((gn)/(gs))^2,{ix,1,gn},{iy,1,gn}];

Module[{g,gs,gn},
CFDLaplace=Compile[{{g,_Real,2},{gs,_Real},{gn,_Integer}},
Table[(
If[ix==1 || ix==gn || iy==1 || iy==gn,
g[[If[ix>1,ix-1,gn],iy]]+g[[If[ix<gn,ix+1,1],iy]]+
g[[ix,If[iy>1,iy-1,gn]]]+g[[ix,If[iy<gn,iy+1,1]]]
,
g[[ix-1,iy]]+g[[ix+1,iy]]+
g[[ix,iy-1]]+g[[ix,iy+1]]] 
-4*g[[ix,iy]]),{ix,1,gn},{iy,1,gn}]*((gn)/(gs))^2
]]

ff[x_]:=x^3-x;
Module[{m1},
FDCahnHillard[g_,eps_,D_]:=D FDLaplace[Table[ff[g[[ix,iy]]],{ix,1,gn},{iy,1,gn}]-eps*FDLaplace[g]];
]

initialization(optimized)

Module[{g,gs,gn},
CFDLaplace=Compile[{{g,_Real,2},{du,_Real}},
(RotateLeft[g,1]+RotateLeft[g,-1]+RotateLeft[g,{0,1}]+RotateLeft[g,{0,-1}]-4*g)/du^2
]]

CFDGradSquare=Compile[{{g,_Real,2},{du,_Real}},
((RotateLeft[g,1]-RotateLeft[g,-1])^2+(RotateLeft[g,{0,1}]-RotateLeft[g,{0,-1}])^2)/(4du^2)
]

Module[{g,gs,gn},
gl2=1./3;
CFDLaplace2=Compile[{{g,_Real,2},{du,_Real}},
(
0.5*gl2*(RotateLeft[g,{1,1}]+RotateLeft[g,{1,-1}]+RotateLeft[g,{-1,1}]+RotateLeft[g,{-1,-1}])+
(1-gl2)(RotateLeft[g,1]+RotateLeft[g,-1]+RotateLeft[g,{0,1}]+RotateLeft[g,{0,-1}])-(4(1-gl2)+2gl2)*g)/du^2
]]

du=1;
Module[{fCH},
fCH=Compile[{{g,_Real,2}},(g*g-1)*g];
gCH=Compile[{{g,_Real,2}},(g*g-1)^2/4];
CFDCahnHillard[g_,eps_]:=CFDLaplace2[fCH[g]-eps*CFDLaplace2[g,du],du];  
CFDCahnHillardG[g_,eps_]:=Total[Flatten[gCH[g]+eps/2*CFDGradSquare[g,du]  ]]*du^2;
]

dynamic

random initialization of the grid

gs=1;gn=50;du=gs/gn; cr=1.0;
g=Table[(RandomReal[]-0.5)*2,{ix,1,gn},{iy,1,gn}];
lg=CFDLaplace[g,du];
l2g=CFDLaplace[CFDLaplace[g,du],du];
ch=CFDCahnHillard[g,0.5^2];
xg=g;

Manipulate[Dynamic[fig=Show[
Switch[what,"c(x,y)",
ListDensityPlot[xg,
ColorFunction->(ColorData["ThermometerColors"][(#1/cr+1)/2]&),ColorFunctionScaling->False,
DataRange->{{0,1},{0,1}},PlotRange->1.1{-cr,cr},FrameLabel->{"x","y"}],
"(Subscript[Δg, 0])(x,y)",
ListDensityPlot[gCH[xg],
ColorFunction->(ColorData["BeachColors"][(1-#1)]&),
DataRange->{{0,1},{0,1}},PlotRange->All,FrameLabel->{"x","y"}],
"(|∇g|)^2",
ListDensityPlot[CFDGradSquare[xg,du],
ColorFunction->(ColorData["BeachColors"][(1-#1)]&),
DataRange->{{0,1},{0,1}},PlotRange->All,FrameLabel->{"x","y"}]
],
Graphics[{Text["ΔG="<>ToString[SetPrecision[CFDCahnHillardG[xg,e0^2],2]],Offset[{-5,-5},{1,1}],{1,1},Background->White]}],
ImageSize->6*72
]
],
Style["Initialization",Bold],
{{n0,50,"Grid points"},25,250,5,Appearance->"Labeled"},
{{rng,0.5,"Random noise strength"},0.01,1,0.1},
{{phx,0,"Balance between phases"},-0.5,0.5,0.1},
Button["reset",
xg=g=Table[phx+(RandomReal[]*2-1)rng *(1-Abs[phx]),{ix,1,n0},{iy,1,n0}]; du=gs/n0;
,ImageSize->80],
Style["Runtime options",Bold],
{{e0,0.03,"Interface range"},Dynamic[du/2],0.5,Appearance->"Labeled"},
{{ldt,-0.5,"Time Step"},-3,0,Row[{Manipulator[##],"  ",Dynamic[Round[10^ldt,0.1],Appearance->"Frameless",BaseStyle->"Label"]}]&},
{{mxst,1000,"Number of Steps"},{10,100,1000,10000,100000}},
{{what,"c(x,y)","Show"},{"c(x,y)","(Subscript[Δg, 0])(x,y)","(|∇g|)^2"}},
Button["Go!",
xg=g;
Do[
sdt=10^ldt *du^4/(32e0^2);
ch=CFDCahnHillard[g,e0^2];
g+=ch*sdt;
If[Flatten[g].Flatten[g]+(Flatten[xg].Flatten[xg])-2(Flatten[g].Flatten[xg])>1*^-3,xg=g];
,{j,1,10*mxst}]
,ImageSize->80,Method->"Queued"],
Initialization:>{
CFDGradSquare=Compile[{{g,_Real,2},{du,_Real}},
((RotateLeft[g,1]-RotateLeft[g,-1])^2+(RotateLeft[g,{0,1}]-RotateLeft[g,{0,-1}])^2)/(4du^2)
];
gs=1.0;
gl2=1./3;cr=1.0;
xg=g=Table[(RandomReal[]-0.5)*2rng,{ix,1,n0},{iy,1,n0}]; du=gs/n0;
CFDLaplace2=Compile[{{g,_Real,2},{du,_Real}},
(
0.5*gl2*(RotateLeft[g,{1,1}]+RotateLeft[g,{1,-1}]+RotateLeft[g,{-1,1}]+RotateLeft[g,{-1,-1}])+
(1-gl2)(RotateLeft[g,1]+RotateLeft[g,-1]+RotateLeft[g,{0,1}]+RotateLeft[g,{0,-1}])-(4(1-gl2)+2gl2)*g)/du^2
];
fCH=Compile[{{g,_Real,2}},(g*g-1)*g];
gCH=Compile[{{g,_Real,2}},(g*g-1)^2/4];
CFDCahnHillard[g_,eps_]:=CFDLaplace2[fCH[g]-eps*CFDLaplace2[g,du],du];  
CFDCahnHillardG[g_,eps_]:=Total[Flatten[gCH[g]+eps/2*CFDGradSquare[g,du]  ]]*du^2;
}
]

Export["/tmp/ch-7.png",Show[fig,ImageSize->4*72],"PNG",ImageResolution->200]

GIBBS SURFACE EXCESS

initialization

fe[{x0_,s_,ca_,cb_},x_]:=cb+(ca-cb)/(1+Exp[(x-x0)/s])
ife[{x0_,s_,ca_,cb_},{xa_,xb_}]:=ca*(-xa+xb)+(ca-cb)*s*(Log[E^(x0/s)+E^(xa/s)]-Log[E^(x0/s)+E^(xb/s)])

dynamic

{xa,xb}={-5,5};
Manipulate[Show[
Plot[{fe[{x1,s1,ca1,cb1},x],
fe[{x2,s2,ca2,cb2},x]},
{x,xa,s},Filling->{1->ca1,2->ca2},PlotStyle->None,FillingStyle->Directive[Black,Opacity[0.1]]],
Plot[{fe[{x1,s1,ca1,cb1},x],
fe[{x2,s2,ca2,cb2},x]},
{x,s,xb},Filling->{1->cb1,2->cb2},PlotStyle->{RGBColor[0.2,0.5,1],RGBColor[0.8,0,0]}],
Graphics[{Thick,Black,Line[{{s,0},{s,1}}],
g1=SetAccuracy[ife[{x1,s1,ca1,cb1},{xa,xb}]-(s-xa)ca1-(xb-s)cb1,3];
g2=SetAccuracy[ife[{x2,s2,ca2,cb2},{xa,xb}]-(s-xa)ca2-(xb-s)cb2,3];
g21=SetAccuracy[g2-g1(ca2-cb2)/(ca1-cb1),3];
Text["\!\(\*SubsuperscriptBox[\(Γ\), \(1\), \(\\ \)]\)="<>ToString[If[Abs[g1]<=0.01,0,g1]],Offset[{-5,-15},{xb,1}],{1,-1}],
Text["\!\(\*SubsuperscriptBox[\(Γ\), \(2\), \(\\ \)]\)="<>ToString[If[Abs[g2]<=0.01,0,g2]],Offset[{-5,-35},{xb,1}],{1,-1}],
Text["\!\(\*SubsuperscriptBox[\(Γ\), \(2\), \(\(\\ \)\(1\)\)]\)="<>ToString[If[Abs[g21]<=0.01,0,g21]],Offset[{-5,-55},{xb,1}],{1,-1}]
}],
Plot[
{
fe[{x1,s1,ca1,cb1},x],
fe[{x2,s2,ca2,cb2},x]
},{x,xa,xb},PlotStyle->{RGBColor[0.2,0.5,1],RGBColor[0.8,0,0]}
],PlotRange->{{xa,xb},{0,1}},Frame->True,Axes->False,FrameLabel->{"x","c"},ImageSize->7*72
], 
Style["Concentration profiles",Bold],
{{ca1,1,"Solvent concentration in the first phase"},0.001,1,Appearance->"Labeled"},
{{cb1,0.1,"Solvent concentration in the second phase"},0.001,1,Appearance->"Labeled"},
{{x1,0,"Solvent position"},-3,3,Appearance->"Labeled"},
{{s1,0.5,"Solvent layer width"},0.1,1,Appearance->"Labeled"},
{{ca2,0.3,"Solute concentration in the first phase"},0.001,1,Appearance->"Labeled"},
{{cb2,0.02,"Solute concentration in the second phase"},0.001,1,Appearance->"Labeled"},
{{x2,1,"Solute position"},-3,3,Appearance->"Labeled"},
{{s2,0.25,"Solute layer width"},0.1,1,Appearance->"Labeled"},
Style["Gibbs Dividing Surface",Bold],
{{s,0.2,"Dividing Surface position"},-1,1,Appearance->"Labeled"},SaveDefinitions->True
]

Regular solution model

initialization

G[{mua_,mub_,Om_,RT_},x_]:=x(mua+RT Log[x]) +(1-x) (mub +RT Log[1-x]) +Om x (1-x)
dG[{mua_,mub_,Om_,RT_},x_]:=mua-mub+Om (1-x)-Om x-RT Log[1-x]+RT Log[x]

dynamic

{x0,x1}={x0,x1}/.FindRoot[{G[{mua,mub ,Om,RT},x1]-x1 dG[{mua,mub,Om,RT},x1]==G[{mua,mub,Om,RT},x0]-x0 dG[{mua,mub,Om,RT},x0],
 dG[{mua,mub,Om,RT},x1]== dG[{mua,mub,Om,RT},x0]},{{x0,0.01},{x1,0.99}}];

Manipulate[
Module[{rl,dx0,dx1,x1,x0,xa,xb,x0t,x1t,pd},
Off[FindRoot::lstol];Off[FindRoot::jsing];
{x0,x1}={x0,x1}/.FindRoot[{G[{mua,mub ,Om,RT},x1]-x1 dG[{mua,mub,Om,RT},x1]==G[{mua,mub,Om,RT},x0]-x0 dG[{mua,mub,Om,RT},x0],
 dG[{mua,mub,Om,RT},x1]== dG[{mua,mub,Om,RT},x0]},{{x0,0.01},{x1,0.99}}];
fig=GraphicsRow[Flatten[{
Show[
Plot[{
x mua+(1-x) mub,
G[{mua,mub,Om,RT},x]},{x,0,1},PlotRange->{{-0.01,1.01},All},Frame->True,Axes->False,FrameLabel->{"x","G"},FrameStyle->Directive[Black,Dashing[{}]],PlotStyle->{Directive[Thin,Dashed,Red],Directive[{Thick,RGBColor[0.5,0,0]}]} ],
If[Abs[x0-x1]>1*^-3,
Show[
Plot[G[{mua,mub,Om,RT},x0]+dG[{mua,mub,Om,RT},x0](x-x0),{x,0,1},PlotRange->{{-0.01,1.01},All},Frame->True,Axes->False,FrameLabel->{"x","G"},PlotStyle->Directive[{Thick,Dotted,RGBColor[0.9,0.3,0]}]],
Graphics[{Red,PointSize[0.025],Point[{x0,G[{mua,mub,Om,RT},x0]}],Point[{x1,G[{mua,mub,Om,RT},x1]}]}]
],
Graphics[{}]],ImageSize->6*72
],
pd=If[dopd,
Flatten[Table[
{xa,xb}={x0t,x1t}/.FindRoot[{G[{mua,mub ,Om,t},x1t]-x1t dG[{mua,mub,Om,t},x1t]==G[{mua,mub,Om,t},x0t]-x0t dG[{mua,mub,Om,t},x0t],
 dG[{mua,mub,Om,t},x1t]== dG[{mua,mub,Om,t},x0t]},{{x0t,0.01},{x1t,0.99}}];
If[Abs[xa-xb]>1*^-3,Re[{{xa,t},{xb,t}}],{}]
,{t,0.01,10,0.05}],1]
,{}];
xpd=pd;
If[pd≠{},Show[
ListPlot[Sort[pd],
Joined->True,Frame->True,Axes->False,FrameLabel->{"x","RT"},PlotStyle->{Directive[{Thick,Dashing[{}],RGBColor[0.1,0.3,0.8]}]} ],
Graphics[{Directive[Thin,Dashed,Black],Line[{{0,RT},{1,RT}}]}]
,PlotRange->All],{}]
}]]
], 
{{mua,1.2,"(Subscript[μ, A])"},0,10,Appearance->"Labeled"},
{{mub,1,"(Subscript[μ, B])"},0,10,Appearance->"Labeled"},
{{Om,2.5,"Ω"},-5,5,Appearance->"Labeled"},
{{RT,1,"RT"},0.01,10,Appearance->"Labeled"},
{{dopd,False,"Plot phase diagram"},{True,False}},
SaveDefinitions->True
]