Set Up
I am attempting to make a 2D model of the diffusion of ions in solution due to concentration gradients and an electric field.
The equation that needs to be solved to describe the diffusion is
∂C/∂t=-∇∙(-D∇C-zμC∇ϕ)
C: Concentration of species in a given space
D:Diffusion Coefficient
z:Charge of a species molecule
μ:Mobility of a species molecule
ϕ:Electric potential at a given point
-
The code is as follows:
All constants are set equal to "1"
diffusionC = 1;
charge = 1;
mobility = 1;
The equation that describes the electric field is defined
\[CapitalPhi][t_, x_, y_] := -2.5*x + 5;
The gradients inside the parenthesis are defined
diffGrad = Grad[u[t, x, y], {x, y}];
elecGrad = Grad[\[CapitalPhi][t, x, y], {x, y}];
The terms inside the parenthesis are defined
equationP1 = -diffusionC*(diffGrad[[1]] + diffGrad[[2]]);
equationP2 = charge*mobility*u[t, x, y]*elecGrad[[1]];
The final gradient outside the parenthesis is defined
fullGrad = Grad[equationP1 + equationP2, {x, y}]
The full equation is defined with flux boundary conditions for the top, bottom, and right sides
equation =
D[u[t, x, y], t] + fullGrad[[1]] ==
NeumannValue[0, y == 0 || y == 2] + NeumannValue[3, x == 2]
The area over which the equation is to be solved is defined
R = ImplicitRegion[0 <= x <= 2 && 0 <= y <= 2, {x, y}];
The initial condition and boundary condition for the left side are defined
ic = u[0, x, y] == 1
bc = u[t, 0, y] == 8 - 7*Exp[-1000*t]
The equation is solved
sol = NDSolve[{equation, ic, bc}, {u, \[CapitalPhi]}, {t, 0,
10}, {x, y} \[Element] R]
The solution is animated
ListAnimate@
Table[Plot3D[u[t, x, y] /. sol, {x, 0, 2}, {y, 0, 2},
PlotRange -> All, AxesLabel -> {x, y, u}], {t, 0, 10, .025}]
ListAnimate@
Table[Plot3D[\[CapitalPhi][t, x, y] /. sol, {x, 0, 2}, {y, 0, 2},
PlotRange -> All, AxesLabel -> {x, y, u}], {t, 0, 10, .025}]
Problem
Upon running the line to solve the equation, the following errors are presented.
NDSolve::femcscd: The PDE is convection dominated and the result may not be
stable. Adding artificial diffusion may help. >>
NDSolve::ndsz: At t == 0.06490840231453537`, step size is effectively zero;
singularity or stiff system suspected. >>
The solution given is:
which only solves up until the time given in the error rather than to "10" as defined in the code.
I request a solution to resolve the errors and get the model to solve properly for the desired time range. If additional information is required, please do let me know. Thank you for your time and consideration.
