Group Abstract

Message Boards

WOLFRAM COMMUNITY

8K 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:

Alan SAILLET

Alan SAILLET, CERN

Posted 5 years ago

POSTED BY: Alan SAILLET

Konstantin Solovev

Konstantin Solovev, Sobolev Institute of Geology and Mineralogy

Posted 5 years ago

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

POSTED BY: Konstantin Solovev

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

Alan SAILLET

Alan SAILLET, CERN

Posted 5 years ago

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

POSTED BY: Alan SAILLET

Konstantin Solovev

Konstantin Solovev, Sobolev Institute of Geology and Mineralogy

Posted 5 years ago

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.