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).

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.

1.7 Modify frictional force to only resist movement¶

As Eric pointed out, at the moment the frictional force is always oriented in the same direction and is always the maximum value.

In really two things are true about friction:

  1. It only resits motion
  2. It never exceeds the pushing force.

To resolve this issue Eric had the idea to model tau as a sigmoid. In effect, we will take the maximal tau we calculated, and modify it with the function

$$\frac{2*\tau_{max}}{1 + e^{-k*u}} - \tau_{max}$$

This function will effectively smooth out tau around u=0 to avoid to case where the frictional force causes the velocity to actually change signs (which of course is not physical).

The downside of this method is that when u=0 we see slip in the system which is not exactly physical. Furthermore, to minimize that slip I use a large value of $k = 1e3$ which then causes me to need a very granular timestep to avoid weird oscillations.

1.8 Make sure get_jacobian for variable viscosity shear stress is correct.¶

In order to encourage get_jacobian to be called, I need to make two changes:

  1. Set the source_treatment to implicit for the frictional viscosity shear stress term.
  2. Set the timestepper type to Strang

Looking at the documentation for the timestepper, I came across this comment:

		Explicit Schemes:
		-----------------
		- Forward Euler (FE)
		- 4th-order Runge Kutta (RK4)
		- Low storage 4th-order Runge Kutta (LSRK4)
		- Strong-stability preserving 3rd-order Runge Kutta (SSPRK3)
		- Arbitrary DERivatives in space and time (ADER)
			-> used in tandem with ADERDG solver

		Operator Splitting Type Schemes:
		--------------------------------
		- Strang Splitting (Strang)
		- Simpler Splitting (Simpler)

In the below example, the viscous term is interpreted implicitly. Some of the challenges that introduces is that we can no longer easily get the viscous drag term.

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

HTML(ani.to_html5_video())

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

For example, we can see below that we only see the viscous drag term when we first plot the source terms.

In [16]:
folder = "slip_variable_addition"
file_prefix = "short_plug_friction_resistive_force"

solver = readwritedatafiles.read_data_file(f"{folder}/{file_prefix}_{0}.pkl")

print(solver.physics.source_terms)


solver = readwritedatafiles.read_data_file(f"{folder}/{file_prefix}_{1}.pkl")

print(solver.physics.source_terms)

[<physics.multiphasevpT.functions.SlipSource object at 0x12664e440>, <physics.multiphasevpT.functions.FrictionVolSlip object at 0x125fb4050>, <physics.multiphasevpT.functions.FrictionVolFracVariableMu object at 0x126688350>]
[<physics.multiphasevpT.functions.SlipSource object at 0x12715ea80>, <physics.multiphasevpT.functions.FrictionVolSlip object at 0x12731fcf0>]

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.