# Numerical solution of multicomponent difussion in mineral

Posted 9 months ago
1123 Views
|
7 Replies
|
2 Total Likes
|
 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 elementThere 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. Attachments:
7 Replies
Sort By:
Posted 9 months ago
 Can you add comment lines between each line? Because just like this, it is hard to understand your code
Posted 9 months ago
 I tried to clarify my code and I added comments for lines. Thank you for your response.
Posted 9 months 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.
Posted 9 months ago
 My first question (for my personal knowledge): Why do you want to solve it using FEM only?If you do not specify < "FiniteElement">> in NDSolveValue, how can you be sure it is choosing automatically the FEM method? Actually, By specifying < "FiniteElement">>, the error is not the same.As you said, using FEM, you cannot solve 2nd order derivatives. Since <> 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 9 months 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} = NDSolveProcessEquations[{eqnFEM, bound1FEM, bound2FEM, initFEM}, {XMg, XFe, XCa}, {t, 0, 300000000}, {x, xbegin, xend}] you will see {NDSolveStateData["<" 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.