Group Abstract

Message Boards

WOLFRAM COMMUNITY

8.2K Views

7 Replies

2 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Earth Science Mathematics Physics Graphics and Visualization Wolfram Language Numerical Computation

Numerical solution of multicomponent difussion in mineral

Konstantin Solovev

Konstantin Solovev, Sobolev Institute of Geology and Mineralogy

Posted 5 years ago

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} = NDSolve`ProcessEquations (={NDSolve`StateData["<" 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. Attachments: Diffusion.nb

POSTED BY: Konstantin Solovev

7 Replies

Sort By:

Konstantin Solovev

Konstantin Solovev, Sobolev Institute of Geology and Mineralogy

Posted 5 years ago

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.

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}];

 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.

POSTED BY: Konstantin Solovev

Konstantin Solovev

Konstantin Solovev, Sobolev Institute of Geology and Mineralogy

Posted 5 years ago

POSTED BY: Konstantin Solovev

Konstantin Solovev

Konstantin Solovev, Sobolev Institute of Geology and Mineralogy

Posted 5 years ago

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.