Week of 2025.2.23¶

Last week we finished setting the simulation up to have a numeric steady state solution. Let's recap that here:

1.1 Steady state solution reviewed¶

The force balance for the plug can be written like so

$$ M \ddot{s} = A (p0 + \Delta p(s(t))) - P_{atm} A - 2 \pi R L \tau(s(t)) $$

We are basically just stating that the force exerted on the plug is the delta in pressure multiplied by the cross section surface area, subtracted by the shear stress multiplied by the surface area of the conduit.

Equilibrium is defined as

$$ \begin{align} 0 = P_0 A - P_{atm}A - 2 \pi R L \tau(s)\\ P_0 = p_{atm} + \frac{2 \pi R L \tau_p}{A} \\ P_0 = p_{atm} + \frac{2 L \tau_p}{R} \end{align} $$

Solving that we get

$$ P_{-500} = p_{atm} + \frac{2 L \tau_p}{R} \\ P_{-500} = 1e5 [Pa] + \frac{2 * 500 [m] * 2e5 [Pa]}{50 * [m]} \\ P_{-500} = 4.1e6 [Pa] $$

1.2 Imports¶

In [1]:
import os
import numpy as np
import matplotlib.pyplot as plt

# Modify base path for depending on your file structure.
BASE_PATH = "/Users/paxton/git"

# Specify path where .pkl files are located
target_dir = f"{BASE_PATH}/quail_volcano/scenarios/simple_1D_test"
# Specify path for Quail source code
source_dir = f"{BASE_PATH}/quail_volcano/src"
# Change to working directory
os.chdir(target_dir)


# Import quail modules
os.chdir(source_dir)

import meshing.tools as mesh_tools

import numerics.helpers.helpers as helpers
import numerics.timestepping.tools as stepper_tools

import physics.zerodimensional.zerodimensional as zerod
import physics.euler.euler as euler
import physics.navierstokes.navierstokes as navierstokes
import physics.scalar.scalar as scalar
import physics.chemistry.chemistry as chemistry
import physics.multiphasevpT.multiphasevpT as multiphasevpT

import processing.readwritedatafiles as readwritedatafiles
import processing.post as post
import processing.plot as plot
import processing.animate as animate

import solver.DG as DG
import solver.ADERDG as ADERDG
import solver.tools as solver_tools

import time
from IPython.display import HTML
import multiprocessing as mp  
from multidomain import Domain, Observer

os.chdir(target_dir)

1.3 Updated geophysics¶

  • Reduce plug height to 100m
  • Reduce radius to 10m
  • Reduce bulk modulus by factor of 10 to 1e9

$$ \begin{align} P_0 &= p_{atm} + \frac{2 L \tau_p}{R} \\ P_0 &= 1e5 [Pa] + \frac{2 * 100 [m] * 2e5 [Pa]}{10m} \\ P_0 &= 4.1e6 [Pa] \end{align} $$

We see that those changes did not affect the steady state because the length of the plug and radius were reduced proportionally.

Problem: Because the boundary between melt and plug is no longer symetric, the diffusion error causes the boundary to shift more than we want.

Solution:

  • Let's start using the 'slip' variable in conjunction with the initial boundary to track the plug/melt boundary.

Note:

  • I currently do not allow a negative velocity to act as a source term on slip. Instead I return np.maxmimum(rho*u, 0).
In [3]:
ani = animate.animate_conduit_pressure("slip_variable_addition", iterations=100, file_prefix="split_domain_short_plug")

HTML(ani.to_html5_video())

/Users/paxton/git/quail_volcano/scenarios/simple_1D_test
Out[3]:
Your browser does not support the video tag.

1.4 Stability analysis for dynamic system¶

$$ \begin{align} M \ddot{s} = A \Delta p(s(t)) + 2 \pi R L (\tau_p - \tau(s)) \\ M \ddot{s} = \frac{-A^2s(t)}{S} + 2 \pi R L (\tau_p - \tau(s)) \end{align} $$

Let's consider the two primary competing forces for stability. First there is the shear stress drag as a result of slip. That can be express as:

$$ \begin{align} F_d &= 2 \pi R L (\tau_p - \tau(s)) \\ \end{align} $$

Note that $\tau(s) = \tau_r - (\tau_r - \tau_p) \exp{-s(t)/D_c}$. Also note that with a first order taylor expansion, $\exp(-s(t)/D_c) = 1 - s(t)/D_c$. So a first order estimate of $\tau(s)$ is: $\tau(s) = \tau_p + (\tau_r - \tau_p)(\frac{s(t)}{D_c})$.

Applying that approximation we get:

$$ \begin{align} F_d &= 2 \pi R L (\tau_p - \tau_r)(\frac{s(t)}{D_c}) \\ \end{align} $$

Now, let's consider the force applied by the pressure on the left hand side of the domain

$F_p = \frac{-A^2 s(t)}{S}$ where $S = \frac{A L_{magma}}{K}$ so that becomes:

$$F_p = \frac{-A K s(t)}{L_{magma}}$$

Now, to think about stability lets call the system unstable when:

$$ \frac{\frac{d F_d}{d s}}{\frac{d F_p}{d S}} < 1$$

The derivatives are easy to take:

$$ \begin{align} \frac{d F_d}{d s} &= 2 \pi R L_{plug} (\frac{\tau_p - \tau_r}{D_c}) \\ \frac{d F_p}{d s} &= \frac{A K}{L_{magma}} \end{align} $$

$$ \begin{align} \frac{2 \pi R L_{plug} (\frac{\tau_p - \tau_r}{D_c})}{\frac{A K}{L_{magma}}} > 1 \\ \frac{2 \pi R L_{plug} L_{magma} (\tau_p - \tau_r)}{A K D_c} > 1 \\ \end{align} $$

So we conclude that the system is become more unstable when slip is the range $(0, D_c)$ and

$$ D_c < \frac{2 L_{magma} L_{plug} (\tau_p - \tau_r)}{R K} $$

Plugging in our numbers we get:

$$ D_c < \frac{2 * 100m[m] * 900 [m] * 2e5 [Pa]}{10 [m] * 1e9 [Pa]} $$

$$ D_c < 3.6m $$

What does this tell us? It says something like: once we start slipping, if $D_c < 3.6$ we will continue to slip at least 3.6m until slip > D_c.

Note we use a first order approximation of the exponent function for tau, which may affect these results slightly.

1.5 Numerical example with unstable slip.¶

Let's slightly increase the pressure in the melt and verify that once we start slipping we continue for 3.5 [m] when D_c = 3.5m.

  • Increase pressure in the melt to 10 [Pa] from 4.1
  • Reduce $D_c$ from 10 [m] to 3.5 [m].
In [2]:
ani = animate.animate_conduit_pressure("slip_variable_addition", iterations=100, file_prefix="short_plug_dynamic")

HTML(ani.to_html5_video())

/Users/paxton/git/quail_volcano/scenarios/simple_1D_test
Out[2]:
Your browser does not support the video tag.

1.6 Evolving from stable state with VelocityInlet1D¶

Next, let's experiment with starting at steady state and slowly increasing the pressure from the left with a velocity inlet. This is probably not a very physical example (the melt pressure increases very quickly) but still fun to see that we can evolve from steady state with a velocity inlet condition.

In [3]:
ani = animate.animate_conduit_pressure("slip_variable_addition", iterations=100, file_prefix="short_plug_evolve_with_inlet")

HTML(ani.to_html5_video())

/Users/paxton/git/quail_volcano/scenarios/simple_1D_test
Out[3]:
Your browser does not support the video tag.

Open Questions¶

  1. In theory the slip drag should only ever apply resistive force. It appears in the current simulation that is actually applies negative pushing force. That is a problem. I somewhat mitigated it by only allowing the source term of rho * u to be positive. All the same the velocity graph looks wrong to me.
  2. Verify that my new model for tracking the boundary holds up when the boundary is moving.