I am trying to solve numerically next PDE system: $$\frac{\partial c}{\partial t}=\epsilon\frac{\partial}{\partial x}(\frac{\partial c}{\partial x}+\rho\frac{\partial \varphi}{\partial x}+\frac{vc}{1-vc}\frac{\partial c}{\partial x})$$ $$\frac{\partial \rho}{\partial t}=\epsilon\frac{\partial}{\partial x}(\frac{\partial \rho}{\partial x}+c\frac{\partial \varphi}{\partial x}+\frac{v\rho}{1-vc}\frac{\partial c}{\partial x})$$ $$-\epsilon^2\frac{\partial \varphi^2}{\partial x^2}=\rho$$ with completely blocking boundary conditions at x= ±1; $v$ and $\epsilon$ are constants. However, at some point of calculation the oscillations near the boundaries occur due to the too big gradients near the boundaries. I have found some information that there are several techniques for handling numerical instabilities. (Petrov-Galerkin (SUPG) or Galerkin least-squares (GLS) method). Unfortunately all of them usually described for the case of simple Convection-Diusion equation. Can anybody help me how to adjust some of those methods (or maybe another method) for my case? Also I need to add, that problem is not in the fact that $(1-vc)$ may tended to zero in some points.
$\begingroup$
$\endgroup$
-
$\begingroup$ How do you solve the non-linearity? $\endgroup$ – Dr_Sam Jun 24 '13 at 7:31
-
$\begingroup$ I don't quite understand your question. I can say only, that I'm using special software for solving PDE, not my own C++ code. $\endgroup$ – oleksandrv Jun 24 '13 at 11:57
-
$\begingroup$ You have a non-linear system of PDEs, so I wanted to know how you linearize it to be able to solve using finite elements. $\endgroup$ – Dr_Sam Jun 24 '13 at 12:19
-
$\begingroup$ Newton's method. $\endgroup$ – oleksandrv Jun 24 '13 at 13:44
$\begingroup$
$\endgroup$
If you are using a finite element method to solve this problem, why not put a few grid points in the boundary layer? 1-D implementations are usually pretty trivial to create that support meshes of varying cell size.
-
$\begingroup$ It is highly nonlinear task, therefore it requires very dense grid which in turn leads to the quite long computation time. I have already tried to do it (increase a density of the grid near the boundaries) within the computational limits of my computer (I mean amount of RAM), however it did not give any positive results. $\endgroup$ – oleksandrv Jun 23 '13 at 20:42
-
$\begingroup$ Maybe if you plotted the troublesome solution for us, gave the values of $\nu$ and $\epsilon$, and gave the mathematical form of the boundary conditions, we could help out a little more. $\endgroup$ – Bill Barth Jun 23 '13 at 23:19
-
$\begingroup$ Boundary conditions can be written in the next form: $$\frac{\partial c}{\partial x}+\rho\frac{\partial \varphi}{\partial x}+\frac{vc}{1-vc}\frac{\partial c}{\partial x}=0$$ $$\frac{\partial \rho}{\partial x}+c\frac{\partial \varphi}{\partial x}+\frac{v\rho}{1-vc}\frac{\partial c}{\partial x}=0$$ $$\varphi\pm\delta\epsilon\frac{\partial \varphi}{\partial x}=\pm v$$ where $\delta$ is a constant which measures the effective thickness of the surface insulating layer (according to the physics model). $\endgroup$ – oleksandrv Jun 24 '13 at 9:27
-
$\begingroup$ To be honest I'm using these equations in a little bit another form (which is a little harder for analysis, for my opinion) $$\frac{\partial cn}{\partial t}=\frac{\partial }{\partial x}(A*cn \frac{\partial }{\partial x}(-B*\phi+C*(\log{cn}-\log{(1-D*(cp+cn))}))$$ $$\frac{\partial cp}{\partial t}=\frac{\partial }{\partial x}(A*cp \frac{\partial }{\partial x}(B*\phi+C*(\log{cp}-\log{(1-D*(cp+cn))}))$$ $$\frac{\partial \phi^2}{\partial x^2}=D(cp-cn)$$ However, this system can be easely converted to previous form with simple change of variables. $\endgroup$ – oleksandrv Jun 24 '13 at 9:27
-
$\begingroup$ Boundary conditions: $$\mathbf{n}\cdot((A*cn \frac{\partial }{\partial x}(-B*\phi+C*(\log{cn}-\log{(1-D*(cp+cn))})))=0$$ $$\mathbf{n}\cdot((A*cp \frac{\partial }{\partial x}(B*\phi+C*(\log{cp}-\log{(1-D*(cp+cn))})))=0$$ $$\mathbf{n}\cdot(80*\frac{\partial \phi}{\partial x})=0$$ Values of parameters: A=7.26E-13 B=96485 C=2436 D=2.21E-4 Initial conditions: $cn$=$cp$=250; $\phi$=3 (smoothly increasing from 0 to 3) Oscillations: imgur.com/OY85XOE $\endgroup$ – oleksandrv Jun 24 '13 at 9:32
