In[]:=
ClearAll["Global`*"]
In[]:=
Needs["NDSolve`FEM`"]
In[]:=
(*inputData*)
In[]:=
(*Xcomp-listofmolefractionsofeachelements(Mg,Fe,Ca,andMn).XMn+XMg+XFe+XCa=1,sothereareonlythreeindependentvariable.IchooseMg,FeandCaasindendentvariables*)​​Xcomp[t_,x_]:={xMg[t,x],xFe[t,x],xCa[t,x],1-xMg[t,x]-xFe[t,x]-xCa[t,x]};​​(*Dself0-listofself-diffusioncoefficientsfor4elementsforstandardcondition.Itisusedforcalculationofself-diffusioncoefficientsinintrestingcondition-Dself*)​​Dself0={9.64*
-7
10
,9.11*
-8
10
,8.37*
-10
10
,1.51*
-8
10
};​​(*Q1bar,delV,R,T,P-parameterswhichisusedforcalculationofself-diffusioncoefficientsininterestingcondition-Dself*)​​Q1bar={304300,277600,255700,249300};​​delV={0.53,0.56,0.60,0.60};​​R=8.31446261815324;​​T=1625.15;​​P=63900;​​(*Dself-self-difussioncoefficientswhichareusedforcalculationofmatrixofinterdiffusioncoeffcients*)​​Dself=TableDself0[[i]]*Exp-
Q1bar[[i]]+((P-1)*delV[[i]])
R*T
,{i,1,Length[Dself0]};​​(*Dm-matrixofinterdiffusioncoefficients.Lastcolumnandrowaremadeupofdependentvariables(Mn-dependentvariables).So,weusematrix3x3(withoutlastcolumnandrow).See-Dij.Isimplifytheseformorecompactvariant.*)​​Dm=FullSimplifyTableDself[[i]]*KroneckerDelta[i,j]-
Dself[[i]]*Xcomp[t,x][[i]]
Length[Dself]
∑
l=1
Xcomp[t,x][[l]]*Dself[[l]]
*(Dself[[j]]-Dself[[Length[Dself]]]),{i,1,Length[Dself]},{j,1,Length[Dself]};​​Dij=Dm[[1;;Dimensions[Dm][[2]]-1,1;;Dimensions[Dm][[2]]-1]];​​(*xbeginandxend-coordinatesofintervalwhichisusedinNDSolve*)​​xbegin=0.00188396;​​xend=0.0038575;​​(*xrim-coordinatewhichisusedforinitialconditionsofmolefracturesofelements*)​​xrim=xend-0.0005;​​(*FirstthreevaluesaremolefracturesforMg,Fe,andCaatt=0ininterval[xbegin,xrim).SecondthreevaluesaremolefracutreforMg,FeandCaatt=oininterval[xrim,xend].ItisusedforintialconditionseeinitFEMandinitTPG*)​​list={0.72996,0.20652,0.05635,0.72449,0.21915,0.04910};
In[]:=
(*equationsandconditionsforFEM*)
In[]:=
(*eqnFEM-threePDEsinInactiveforms,becauseorderofPDEsishigherthan1.DijmatrixofdiffusioncoefficientsmultiplybyIdentityMatrix,becauseGradgivesvector.Infact,IthasformasIwroteinmessage.IsawguideforHeattransferwhenwrotetheseequations.*)​​eqnFEM=TableD[Xcomp[t,x][[i]],t]+Inactive[Div]
3
∑
j=1
(-Dij[[i,j]]*IdentityMatrix[1]).Inactive[Grad][Xcomp[t,x][[j]],{x}],{x}0,{i,1,3};​​(*bound1FEMandbound2FEM-boundaryconditionswhichsaythatmolefractionsatxbeginandxendareconstant.Infact,wecanwritetheseconditionwithoutDirichletCondition,butasIunderstandNDSolveneedDirichletConditionforselectingFEMasMethod*)​​bound1FEM=Table[DirichletCondition[Xcomp[t,x][[i]]list[[i]],xxbegin],{i,1,3}];​​bound2FEM=Table[DirichletCondition[Xcomp[t,x][[i]]list[[i+3]],xxend],{i,1,3}];​​(*initFEM-intialconditionatt=0asdescribedincommentforlist(seeabove)*)​​initFEM=Table[Xcomp[0,x][[i]]==Piecewise[{{list[[i+3]],x≥xrim},{list[[i]],x<xrim}}],{i,1,3}];​​(*NumericalsolutionofeqnFEMwithbound1FEMandbound2FEMasboundaryconditionsandinitFEMasinitialcondition*)​​z1=NDSolve[{eqnFEM,bound1FEM,bound2FEM,initFEM},{xMg,xFe,xCa},{t,0,30000000000},{x,xbegin,xend}]
NDSolve
:PDE parsing error of Div{{-2.96527×
-19
10
-2.45938×
-36
10
xCaPower[2]}}.Grad$101862+{{0.-2.60715×
-37
10
xCaPower[2]}}.Grad$101863+{{0.+1.31815×
-36
10
xCaPower[2]}}.Grad$101861. Inconsistent equation dimensions.
Out[]=
NDSolve
∇
{x}
·-
1.11983×
-34
10
-1.08118×
-34
10
xCa[t,x]-1.14614×
-35
10
xFe[t,x]
8.59047×
-18
10
-8.29394×
-18
10
xCa[t,x]-8.79228×
-19
10
xFe[t,x]+4.44529×
-18
10
xMg[t,x]
.
∇
{x}
xMg[t,x]+0.-
1.08118×
-34
10
xMg[t,x]
8.59047×
-18
10
-8.29394×
-18
10
xCa[t,x]-8.79228×
-19
10
xFe[t,x]+4.44529×
-18
10
xMg[t,x]
.
∇
{x}
xCa[t,x]+0.-
1.14614×
-35
10
xMg[t,x]
8.59047×
-18
10
-8.29394×
-18
10
xCa[t,x]-8.79228×
-19
10
xFe[t,x]+4.44529×
-18
10
xMg[t,x]
.
∇
{x}
xFe[t,x]+
(1,0)
xMg
[t,x]0,
∇
{x}
·0.-
6.39565×
-35
10
xFe[t,x]
8.59047×
-18
10
-8.29394×
-18
10
xCa[t,x]-8.79228×
-19
10
xFe[t,x]+4.44529×
-18
10
xMg[t,x]
.
∇
{x}
xCa[t,x]+-7.71124×
-18
10
-
6.77993×
-36
10
xFe[t,x]
8.59047×
-18
10
-8.29394×
-18
10
xCa[t,x]-8.79228×
-19
10
xFe[t,x]+4.44529×
-18
10
xMg[t,x]
.
∇
{x}
xFe[t,x]+0.+
3.42787×
-35
10
xFe[t,x]
8.59047×
-18
10
-8.29394×
-18
10
xCa[t,x]-8.79228×
-19
10
xFe[t,x]+4.44529×
-18
10
xMg[t,x]
.
∇
{x}
xMg[t,x]+
(1,0)
xFe
[t,x]0,
∇
{x}
·-2.96527×
-19
10
-
2.45938×
-36
10
xCa[t,x]
8.59047×
-18
10
-8.29394×
-18
10
xCa[t,x]-8.79228×
-19
10
xFe[t,x]+4.44529×
-18
10
xMg[t,x]
.
∇
{x}
xCa[t,x]+0.-
2.60715×
-37
10
xCa[t,x]
8.59047×
-18
10
-8.29394×
-18
10
xCa[t,x]-8.79228×
-19
10
xFe[t,x]+4.44529×
-18
10
xMg[t,x]
.
∇
{x}
xFe[t,x]+0.+
1.31815×
-36
10
xCa[t,x]
8.59047×
-18
10
-8.29394×
-18
10
xCa[t,x]-8.79228×
-19
10
xFe[t,x]+4.44529×
-18
10
xMg[t,x]
.
∇
{x}
xMg[t,x]+
(1,0)
xCa
[t,x]0,{DirichletCondition[xMg[t,x]0.72996,x0.00188396],DirichletCondition[xFe[t,x]0.20652,x0.00188396],DirichletCondition[xCa[t,x]0.05635,x0.00188396]},{DirichletCondition[xMg[t,x]0.72449,x0.0038575],DirichletCondition[xFe[t,x]0.21915,x0.0038575],DirichletCondition[xCa[t,x]0.0491,x0.0038575]},xMg[0,x]
0.72449
x≥0.0033575
0.72996
x<0.0033575
0
True
,xFe[0,x]
0.21915
x≥0.0033575
0.20652
x<0.0033575
0
True
,xCa[0,x]
0.0491
x≥0.0033575
0.05635
x<0.0033575
0
True
,{xMg,xFe,xCa},{t,0,30000000000},{x,0.00188396,0.0038575}
In[]:=
(*equationsandconditionsforMethodoflines,TensorProductGrid*)
In[]:=
(*eqnTPG-threePDEs.IthasformasIwrote*)
In[]:=
eqnTPG=TableD[Xcomp[t,x][[i]],t]+
3
∑
j=1
(D[-Dij[[i,j]],x]*D[Xcomp[t,x][[j]],x])+
3
∑
j=1
((-Dij[[i,j]])*(D[Xcomp[t,x][[j]],{x,2}]))0,{i,1,3};​​(*bound1TPGandbound2TPG-boundaryconditionswhichsaythatmolefractionsatxbeginandxendareconstant.Forexample,Xcomp[[1]]=XMg[t,x].Wereplacexbyxbeginorxendforgettingboundaryconditionsattheendofinterval.*)​​bound1TPG=Table[(Xcomp[t,x][[i]]/.x->xbegin)list[[i]],{i,1,3}];​​bound2TPG=Table[(Xcomp[t,x][[i]]/.x->xend)list[[i+3]],{i,1,3}];​​(*initialconditionforPDEsasdescribedincommendforlist*)​​initTPG=Table[(Xcomp[t,x][[i]]/.t0)==Piecewise[{{list[[i+3]],x≥xrim},{list[[i]],x<xrim}}],{i,1,3}];​​(*Numericalsolution*)​​z2=NDSolve[{eqnTPG,bound1TPG,bound2TPG,initTPG},{xMg,xFe,xCa},{t,0,30000000000},{x,xbegin,xend}]
Out[]=
xMgInterpolatingFunction
Domain: 0.,3.×
10
10
,{0.00188,0.00386}
Output: scalar
,xFeInterpolatingFunction
Domain: 0.,3.×
10
10
,{0.00188,0.00386}
Output: scalar
,xCaInterpolatingFunction
Domain: 0.,3.×
10
10
,{0.00188,0.00386}
Output: scalar

In[]:=
(*PlotInterpolatingFunctionforXFe*)​​Plot3D[xFe[t,x]/.z2[[1,2]],{t,0,30000000000},{x,xbegin,xend}]
Out[]=
In[]:=
(*BelowIcheckthatz2useMethodofLineandTensorProductGrid*)​​z2Function=z2[[1,2,2]]​​
Out[]=
InterpolatingFunction
Domain: 0.,3.×
10
10
,{0.00188,0.00386}
Output: scalar

In[]:=
z2Function["ElementMesh"]
Out[]=
None
In[]:=
{state}=NDSolve`ProcessEquations[{eqnTPG,bound1TPG,bound2TPG,initTPG},{XMg,XFe,XCa},{t,0,300000000},{x,xbegin,xend}]
Out[]=
{NDSolve`StateData[<0.>]}
In[]:=
Dij[[1,1]]
Out[]=
1.11983×
-34
10
-1.08118×
-34
10
xCa[t,x]-1.14614×
-35
10
xFe[t,x]
8.59047×
-18
10
-8.29394×
-18
10
xCa[t,x]-8.79228×
-19
10
xFe[t,x]+4.44529×
-18
10
xMg[t,x]