Static or dynamic pressure?

2.4 Static or dynamic pressure?

One could reformulate equation (1) somewhat. To this end, let us say that we would want to represent the pressure $p$ as the sum of two parts that we will call static and dynamic, $p = p_{s} + p_{d}$ . If we assume that $p_{s}$ is already given, then we can replace (1) by

\begin{matrix} - \nabla \cdot 2 η \nabla u + \nabla p_{d} = ρ g - \nabla p_{s} . \end{matrix}

One typically chooses $p_{s}$ as the pressure one would get if the whole medium were at rest – i.e., as the hydrostatic pressure. This pressure can be computed noting that (1) reduces to

\begin{matrix} \nabla p_{s} = ρ (p_{s}, T_{s}, x) g = \bar{ρ} g \end{matrix}

in the absence of any motion where $T_{s}$ is some static temperature field (see also Section 2.6). This, our rewritten version of (1) would look like this:

\begin{matrix} - \nabla \cdot 2 η \nabla u + \nabla p_{d} = [ρ (p, T, x) - ρ (p_{s}, T_{s}, x)] g . \end{matrix}

In this formulation, it is clear that the quantity that drives the fluid flow is in fact the buoyancy caused by the variation of densities, not the density itself.

This reformulation has a number of advantages and disadvantages:

One can notice that in many realistic cases, the dynamic component $p_{d}$ of the pressure is orders of magnitude smaller than the static component $p_{s}$ . For example, in the earth, the two are separated by around 6 orders of magnitude at the bottom of the earth mantle. Consequently, if one wants to solve the linear system that arises from discretization of the original equations, one has to solve it a significant degree of accuracy (6–7 digits) to get the dynamic part of the pressure correct to even one digit. This entails a very significant numerical effort, and one that is not necessary if we can split the pressure in a way so that the pre-computed static pressure $p_{s}$ (or, rather, the density using the static pressure and temperature from which $p_{s}$ results) absorbs the dominant part and one only has to compute the remaining, dynamic pressure to 2 or 3 digits of accuracy, rather than the corresponding 7–8 for the total pressure.
On the other hand, the pressure $p_{d}$ one computes this way is not immediately comparable to quantities that we use to look up pressure-dependent quantities such as the density. Rather, one needs to first find the static pressure as well (see Section 2.6) and add the two together before they can be used to look up material properties or to compare them with experimental results. Consequently, if the pressure a program outputs (either for visualization, or in the internal interfaces to parts of the code where users can implement pressure- and temperature-dependent material properties) is only the dynamic component, then all of the consumers of this information need to convert it into the total pressure when comparing with physical experiments. Since any code implementing realistic material models has a great many of these places, there is a large potential for inadvertent errors and bugs.
Finally, the definition of a reference density $ρ (p_{s}, T_{s}, x)$ derived from static pressures and temperatures is only simple if we have incompressible models and under the assumption that the temperature-induced density variations are small compared to the overall density. In this case, we can choose $ρ (p_{s}, T_{s}, x) = ρ_{0}$ with a constant reference density $ρ_{0}$ . On the other hand, for more complicated models, it is not a priori clear which density to choose since we first need to compute static pressures and temperatures – quantities that satisfy equations that introduce boundary layers, may include phase changes releasing latent heat, and where the density may have discontinuities at certain depths, see Section 2.6.
Thus, if we compute adiabatic pressures and temperatures ${\bar{p}}_{s}, {\bar{T}}_{s}$ under the assumption of a thermal boundary layer worth 900 Kelvin at the top, and we get a corresponding density profile $\bar{ρ} = ρ ({\bar{p}}_{s}, {\bar{T}}_{s}, x)$ , but after running for a few million years the temperature turns out to be so that the top boundary layer has a jump of only 800 Kelvin with corresponding adiabatic pressures and temperatures ${\hat{p}}_{s}, {\hat{T}}_{s}$ , then a more appropriate density profile would be $\hat{ρ} = ρ ({\hat{p}}_{s}, {\hat{T}}_{s}, x)$ .
The problem is that it may well be that the erroneously computed density profile $\hat{ρ}$ does not lead to a separation where $| p_{d} | « | p_{s} |$ because, especially if the material undergoes phase changes, there will be entire areas of the computational domain in which $| ρ - {\hat{ρ}}_{s} | « | ρ |$ but $| ρ - {\bar{ρ}}_{s} | «̸ | ρ |$ . Consequently the benefits of lesser requirements on the iterative linear solver would not be realized.

We do note that most of the codes available today and that we are aware of split the pressure into static and dynamic parts nevertheless, either internally or require the user to specify the density profile as the difference between the true and the hydrostatic density. This may, in part, be due to the fact that historically most codes were written to solve problems in which the medium was considered incompressible, i.e., where the definition of a static density was simple.

On the other hand, we intend ASPECT to be a code that can solve more general models for which this definition is not as simple. As a consequence, we have chosen to solve the equations as stated originally – i.e., we solve for the full pressure rather than just its dynamic component. With most traditional methods, this would lead to a catastrophic loss of accuracy in the dynamic pressure since it is many orders of magnitude smaller than the total pressure at the bottom of the earth mantle. We avoid this problem in ASPECT by using a cleverly chosen iterative solver that ensures that the full pressure we compute is accurate enough so that the dynamic pressure can be extracted from it with the same accuracy one would get if one were to solve for only the dynamic component. The methods that ensure this are described in detail in [KHB12] and in particular in the appendix of that paper.

Note: By default, ASPECT uses the full pressure in the equations, and only prescribing density deviations from a reference state on the right-hand side of (1) would lead to negative densities in the energy equation (3). However, when using one of the approximations described in Section 2.10, the energy balance uses the reference density

\bar{ρ}

instead of the full density, which makes it possible to formulate the Stokes system in terms of the dynamic instead of the full pressure. In order to do this, one would have to use a material model (see Section 6.3.1) in which the density is in fact a density variation, and then the pressure solution variable would only be the dynamic pressure.

[next] [prev] [prev-tail] [front] [up]