Numerical solution of multicomponent difussion in mineral

Questions
, , ,
1

Hi,
I am trying to get numerical solution of multicomponent diffusion in one dimension in mineral. I have system of transient PDEs with transient non-linear coefficient, boundary and intial conditions (see image with general equation and information below).

\frac{\partial}{\partial t}X_{i}=\frac{\partial}{\partial x}\sum_{j=1}^{n-1}D_{ij}\frac{\partial}{\partial x}X_{j}

where D_{ij} is

D_{ij}=D_{i}\delta_{ij}-\frac{D_{i}X_{i}Z_{j}Z_{i}(D_{j}-D_{n})}{\sum_{l=1}^nX_{l}Z_{l}^2D_{l}}

where Z is charge of element

There are four components with equal charge (XMg+XFe+XCa+XMn=1, n=4, and we consider Mn as dependent component), interval of diffusion (in meter, x) and time interval (in sec, t), intial concentration (at t=0, Piecewise function), assumption that concentrations of components at the ends of interval do not change (boundary condition) and self-diffusion coefficient (D_{i} or D_{j}).

So, XMg, XFe, and XCa depend on t and x, and D_{ij} depends on XMg[t,x], XFe[t,x], and XCa[t,x]. Depend variables occur as in denominator as in numerator in D_{ij}, so coefficients of system of 3 transient PDE is transient and non-linear.
I read that Wolfram Mathematica v 12.1 can solve such equations, so I tried to implement code for getting of numerical solutions.

In general, I got not bad results with using Method of Lines and TensorProductGrid (NDSolve chose it automatically). I checked it by “ElementMesh” (=None) and {state} = NDSolveProcessEquations (={NDSolveStateData[“<” 0. “>”]}).
But I want to get result with Method of Lines and FEM or with purely FEM. Here I have a problem.

When I re-write equations in the Coefficient Form (with Div, Grad, Inactive), condition in DirichletCondition form and try to evalute NDSolve, Wolfram Mathematica send me error message (see in attached notebook) “NDSolveValue::femper”; “PDE parsing error of Div … Inconsistent equation dimensions.”

What is going wrong? How should I correct my code? If Method of Line with TensorProductGrid works, at least Method of Line with FiniteElement will work too.

P.S.
I tried to prediscribe fuctions for two of three components (for exmaple, XCa[t,x] and XFe[t,x]) and NDSolve can get result for third component (hence, equations are correct), but it did not work with two or more dependent variables.

Open in Wolfram Cloud ↗

Diffusion.nb (282 KB)

2

Can you add comment lines between each line? Because just like this, it is hard to understand your code

3

I tried to clarify my code and I added comments for lines. Thank you for your response.

4

When I replaced

 eqnFEM = Table[D[Xcomp[t, x][[i]], t] + Inactive[Div][Sum[((-D[[i,j]]*IdentityMatrix[1]).Inactive[Grad][Xcomp[[j]],{x}],{j,1,3}],{x}]==0, {i, 1, 3}];
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}];

by

 eqnFEM = Table[D[Xcomp[t, x][[i]], t] + Sum[D[-D[[i,j]],x]*D[Xcomp[[j]],x],{j,1,3}], {i, 1, 3}]+Sum[-D[[i,j]]*D[Xcomp[[j]],{x,2}],{j,1,3}]=NeumannValue[0,x==xend], {i, 1, 3}];
     bound1FEM = 
       Table[DirichletCondition[Xcomp[t, x][[i]] == list[[i]], 
         x == xbegin], {i, 1, 3}];

I fixed the problem. I happy that it works, but I do not understand yet why it fixed the problem. Guide from Wolfram recommends first form, rather than second.

5

My first question (for my personal knowledge): Why do you want to solve it using FEM only?

If you do not specify <<Method → “FiniteElement”>> in NDSolveValue, how can you be sure it is choosing automatically the FEM method? Actually, By specifying <<Method → “FiniteElement”>>, the error is not the same.

As you said, using FEM, you cannot solve 2nd order derivatives. Since <<D[Sum(Dij*D[Xcomp[[j]],x],x]>> is directly a 2nd derivative (you can see that by showing its expression), I think that if you have some results, it has not been obtained using the FEM method

6

Also, in the results you obtain, were can clearly see that there is almost no diffusion over time. Is that what you expected?

7

Time for diffusion in example just 300000000 sec or about 10 years. It is extremly small for diffusion in munerals. If you increase time up to 100 000 000 000 sec or 3000 years, you will see the result. In the most geological problem diffusion run during more than 1 million years. In my notebook, I also implemented fitting of initial, boundary parameters and time for getting time of diffusion for my data. It return about 2500 years of diffusion. In more ussual geological problems, I would get millions years of diffusion.

8

About you first question: I just want to see difference betwen discretization by FEM and TPG.
About you second question: If you re-write example as I wrote in previous message and write

z1Function = z1[[2]]
z1Function["ElementMesh"]

you will see that it is not none. There is Mesh, so FEM is used. If you write

{state} = 
 NDSolve`ProcessEquations[{eqnFEM, bound1FEM, bound2FEM, 
   initFEM}, {XMg, XFe, XCa}, {t, 0, 300000000}, {x, xbegin, xend}]

you will see

{NDSolve`StateData["<" 0. ">"]}

so t is considered as temporal variable and, hence, NDSolve automatically choose Method of Lines and FEM as method of discretization. I try to specify NDSolve use purely FEM, but calculation it too long. I think when NDSolve see Dirichlet Condition and initial problem at t=0. It automatically choose Method of Lines and FEM. If you re-write boundary condition without DirchletCondition, you will see that NDSolve automatically choose Method Of Lines and TPG as method of discretization.