I'm trying to solve the system of equations that relates to blood flows in arteries i.e. flow in elastic tubes.
The system looks as follows
$$\frac{\partial A}{\partial t}+\frac{\partial\left(Au\right)}{\partial x}=0$$
$$ \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+\frac{1}{\rho}\frac{\partial p}{\partial x}=0$$
The first equation is for volume conservation. Fluid velocity is $u$. Pressure $p$ and area $A$ are related directly by a "tube law" so we don't solve for $p$ and the only quantities that are solved for every timestep are $A$ and $u$.
The tube law has a dissipative term $\partial A/\partial t$ so through ${\partial p}/{\partial x}$ and the continuity equation, we end up with a $\partial^2 u /\partial x^2$ term so there's some nice dissipation for the momentum equation.
Now I wanted to quickly solve this using MATLAB's PDEPE and realized that I needed a second spatial derivative for the volume conservation equation as well so I modified it to be
$$\frac{\partial A}{\partial t}+\frac{\partial\left(Au\right)}{\partial x} + \epsilon\frac{\partial^2 A}{\partial x^2}=0$$ where $\epsilon$ is really small.
I ran a couple of cases and am seeing nice physically sensible tube deformations and flows. I compute the residual for the original unadulterated volume conservation equation and the values seem acceptable.
But from a theoretical standpoint, I wanted to know how egregious this error is or if it is an acceptable hack for solving the system? I have converted a PDE that had first order spatial derivatives into one that has second order derivatives. Is the modified system completely different or is it the same as the original as $\epsilon \to 0$? Thanks!
Usually, addition of an artificial dissipation (AD) term will regularize a numerical solution at the cost of accuracy and it is common practice in computational methods. In some cases, it is even necessary to employ AD to obtain a solution for a problem, but this may depend on the discretization scheme you are using.
For example, at transonic flows with shock waves a central difference FV scheme that employs higher order scalar dissipation at smooth regions and 1st order scalar dissipation near shock-waves is very popular in aeronautics (Jameson-Schmid-Turkel scheme) and important to avoid the odd-even decoupling of the central scheme. Other, more advanced options are to use anisotropic scaling of the dissipation coefficient with respect to local flow properties. For instance, matrix dissipation schemes are used in flows with boundary layers where the AD terms are scaled appropriately to obtain sharper boundary layer gradients.
So, informally speaking, you need $\epsilon$ small enough for accuracy but large enough for stability and this $\epsilon$ may strongly depend on local flow gradients/Jacobians and is problem depended.
What you are asking I think is a bit general, but there is the theory of viscosity solutions and the method of vanishing viscosity which is the study of how some estimates of your modified (approximate) solution $(u^{\epsilon}, A^{\epsilon})$ are affected as $\epsilon \to 0$. You may want to check Evan's PDE book Chapter 10 and its references.