From 38957a12fb3610d38d08b8d6f4b72a8a8f37c7ce Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 17 Aug 2009 02:46:52 +0000 Subject: [PATCH] Finish writing the testcase part. git-svn-id: https://svn.dealii.org/trunk@19281 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-32/doc/intro.dox | 106 +++++++++++++++++++++++-- deal.II/examples/step-32/step-32.cc | 22 +++-- 2 files changed, 113 insertions(+), 15 deletions(-) diff --git a/deal.II/examples/step-32/doc/intro.dox b/deal.II/examples/step-32/doc/intro.dox index c8a8ef5b2d..240a1ad131 100644 --- a/deal.II/examples/step-32/doc/intro.dox +++ b/deal.II/examples/step-32/doc/intro.dox @@ -329,7 +329,23 @@ The setup for this program is mildly reminiscent of the problem we wanted to solve in the first place (see the introduction of @ref step_31 "step-31"): convection in the earth mantle. As a consequence, we choose the following data, all of which appears in the program in units of meters and seconds (the -SI system) even if we list them here in other units: +SI system) even if we list them here in other units. + +As a reminder, the equations we want to solve are these: +@f{eqnarray*} + -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& + -\rho \; \beta \; T \mathbf{g}, + \\ + \nabla \cdot {\mathbf u} &=& 0, + \\ + \frac{\partial T}{\partial t} + + + {\mathbf u} \cdot \nabla T + - + \nabla \cdot \kappa \nabla T &=& \gamma, +@f} +augmented by boundary and initial conditions. We then have to choose data for +the following quantities: - -meaningful eta, kappa -do things in meters, second, not km -temperature boundary conditions at inner ring +
  • The density of the earth mantle varies spatially, but not by very + much. $\rho=3300 \frac{\text{kg}}{\text{m}^3}$ is a relatively good average + value. + +
  • The thermal expansion coefficient $\beta$ also varies with depth + (through its dependence on temperature and pressure). Close to the surface, + it appears to be on the order of $\beta=45\cdot 10^{-6} \frac 1{\text{K}}$, + whereas at the core mantle boundary, it may be closer to $\beta=10\cdot + 10^{-6} \frac 1{\text{K}}$. As a reasonable value, let us choose + $\beta=2\cdot 10^{-5} \frac 1{\text{K}}$. + +
  • The second to last parameter we need to specify is the viscosity + $\eta$. This is a tough one, because rocks at the temperatures and pressure + typical for the earth mantle flow so slowly that the viscosity can not be + determined accurately in the laboratory. So how do we know about the + viscosity of the mantle? The most commonly used route is to consider that + during and after ice ages, ice shields form and disappear on time scales + that are shorter than the time scale of flow in the mantle. As a + consequence, continents slowly sink into the earth mantle under the added + weight of an ice shield, and they rise up again slowly after the ice shield + has disappeared again (this is called postglacial + rebound). By measuring the speed of this rebound, we can infer the + viscosity of the material that flows into the area vacated under the + rebounding continental plates. + + Using this technique, values around $\eta=10^{21} \text{Pa\; s} + = 10^{21} \frac{\text{N\; s}}{\text{m}^2} + = 10^{21} \frac{\text{kg}}{\text{m\; s}}$ have been found as the most + likely, though the error bar on this is at least one order of magnitude. + + While we will use this value, we again have to caution that there are many + physical reasons to assume that this is not the correct value. First, it + should really be made dependent on temperature: hotter material is most + likely to be less viscous than colder material. In reality, however, the + situation is even more complex. Most rocks in the mantle undergo phase + changes as temperature and pressure change: depending on temperature and + pressure, different crystal configurations are thermodynamically favored + over others, even if the chemical composition of the mantle were + homogenous. For example, the common mantle material MgSiO3 exists + in its perovskite + structure throughout most of the mantle, but in the lower mantle the + same substance is stable only as post-perovskite. Clearly, + to compute realistic viscosities, we would not only need to know the exact + chemical composition of the mantle and the viscosities of all materials, but + we would also have to compute the thermodynamically most stable + configurations for all materials at each quadrature point. This is at the + time of writing this program not a feasible suggestion. + +
  • Our last material parameter is the thermal diffusivity $\kappa$, which + is defined as $\kappa=\frac{k}{\rho c_p}$ where $k$ is the thermal + conductivity, $\rho$ the density, and $c_p$ the specific heat. For + this, the literature indicates that it increases from around $0.7$ in the + upper mantle to around $1.7 \frac{\text{mm}^2}{\text{s}}$ in the lower + mantle, though the exact value + is not really all that important: heat transport through convection is + several orders of magnitude more important than through thermal + conduction. It may be of interest to know that perovskite, the most abundant + material in the earth mantle, appears to become transparent at pressures + above around 120 GPa (see, for example, J. Badro et al., Science 305, + 383-386 (2004)); in the lower mantle, it may therefore be that heat + transport through radiative transfer is more efficient than through thermal + conduction. + + In view of these considerations, let us choose + $\kappa=1 \frac{\text{mm}^2}{\text{s}} =10^{-6} \frac{\text{m}^2}{\text{s}}$ + for the purpose of this program. + +All of these pieces of equation data are defined in the program in the +EquationData namespace. diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index 5891742008..7f524782fb 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -80,14 +80,17 @@ using namespace dealii; // @sect3{Equation data} - // This program is mainly an extension of - // step-31 to operate in %parallel, so the - // equation data remains the same. + // In the following namespace, we define the + // various pieces of equation data. All of + // these are exhaustively discussed in the + // description of the testcase in the + // introduction: namespace EquationData { - const double eta = 1; - const double kappa = 1e-6; - const double Rayleigh_number = 10; + const double eta = 1e21; + const double kappa = 1e-6; + const double density = 3300; + const double beta = 2e-5; const double R0 = 6371000.-2890000.; const double R1 = 6371000.- 35000.; @@ -2041,8 +2044,11 @@ local_assemble_stokes_system (const typename DoFHandler::active_cell_iterat .quadrature_point(q)); for (unsigned int i=0; i