5
$\begingroup$

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.

$\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
2
$\begingroup$

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.

$\endgroup$
  • $\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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.