]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement the constant density sphere gravity model.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 15 Aug 2009 23:28:45 +0000 (23:28 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 15 Aug 2009 23:28:45 +0000 (23:28 +0000)
git-svn-id: https://svn.dealii.org/trunk@19276 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-32/doc/intro.dox
deal.II/examples/step-32/step-32.cc

index c7d1c09ffb7336df3ba193ca42c409f0dbcea23b..289c84f2b4db8fe0444fdbc27a0c139ac3382fe9 100644 (file)
@@ -428,21 +428,52 @@ SI system) even if we list them here in other units:
     \varphi(\mathbf x)
     =
     \int_{\text{earth}} 
-    -G \frac{\rho(\mathbf x)}{\|\mathbf x-\mathbf y\|} 
+    -G \frac{\rho(\mathbf y)}{\|\mathbf x-\mathbf y\|} 
     \ \text{d}y,
   @f]
   then $\mathbf g(\mathbf x) = -\nabla \varphi(\mathbf x)$. If we assume that
   the density $\rho$ is constant throughout the earth, we can produce an
-  analytical expression for the gravity vector:
-
-XXX
-  Of course, within this problem, we are only interested in the branch that
-  pertains to within the earth.
-
+  analytical expression for the gravity vector (don't try to integrate above
+  equation somehow -- it leads to elliptic integrals; the simpler way is to
+  notice that $-\Delta\varphi(\mathbf x) = 4\pi G \rho
+  \chi_{\text{earth}}(\mathbf x)$ and solving this
+  partial differential equation in all of ${\mathbb R}^3$ exploiting the
+  radial symmetry):
+  @f[
+    \mathbf g(\mathbf x) = 
+    \left\{
+      \begin{array}{ll}
+        -\frac{4}{3}\pi G \rho \|\mathbf x\| \frac{\mathbf x}{\|\mathbf x\|}
+       & \text{for} \ \|\mathbf x\|<R_1, \\
+        -\frac{4}{3}\pi G \rho R^3 \frac{1}{\|\mathbf x\|^2} 
+        \frac{\mathbf x}{\|\mathbf x\|}
+       & \text{for} \ \|\mathbf x\|\ge R_1
+      \end{array}
+    \right.
+  @f]
+  The factor $-\frac{\mathbf x}{\|\mathbf x\|}$ is the unit vector pointing
+  radially inward. Of course, within this problem, we are only interested in
+  the branch that pertains to within the earth, i.e. $\|\mathbf x\|<R_1$. In
+  the program, we therefore only consider the expression
+  @f[
+    \mathbf g(\mathbf x) = 
+        -\frac{4}{3}\pi G \rho \|\mathbf x\| \frac{\mathbf x}{\|\mathbf x\|}
+        =
+        -\frac{4}{3}\pi G \rho \mathbf x
+       =
+       - 9.81 \frac{\mathbf x}{R_1} \frac{\text{m}}{\text{s}^2},
+  @f]
+  where we can infer the last expression because we know Earth's gravity at
+  the surface (where $\|x\|=R_1$).
 
   There are two problems with this, however: (i) The Earth is not homogenous,
   i.e. the density $\rho$ depends on $\mathbf x$; in fact it is not even a
-  function that only depends on the radius $r$. (ii) The density, and by
+  function that only depends on the radius $r$. In reality, gravity therefore
+  does not always decrease as we get deeper: because the earth core is so much
+  denser than the mantle, gravity actually peaks at around $10.7
+  \frac{\text{m}}{\text{s}^2}$ at the core mantle boundary (see <a
+  target="_top" href="http://en.wikipedia.org/wiki/Earth's_gravity">this
+  article</a>). (ii) The density, and by
   consequence the gravity vector, is not even constant in time: after all, the
   problem we want to solve is the time dependent upwelling of hot, less dense
   material and the downwelling of cold dense material. This leads to a gravity
index 6a7b2e25a43739cb4abd1432975d3e81da84e87f..58917420082d28902c6425a64a457051aceb53e1 100644 (file)
@@ -99,7 +99,7 @@ namespace EquationData
   template <int dim>
   Tensor<1,dim> gravity_vector (const Point<dim> &p)
   {
-    return -9.81 * p / p.norm_square();
+    return -9.81 * p / R1;
   }
 
   

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.