2025.1.27 Experiments

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 solver.DG as DG
import solver.ADERDG as ADERDG
import solver.tools as solver_tools

import time
import multiprocessing as mp  
from multidomain import Domain, Observer

os.chdir(target_dir)

from IPython.display import HTML
from matplotlib import animation


VISCOSITY_FRICTION_SOURCE_IDX = 0

def initial_vars(folder):
    solver = readwritedatafiles.read_data_file(f"{folder}/test_output_{0}.pkl") 
    rho = solver.state_coeffs[..., 0:3].sum(axis=-1, keepdims=True)
    pressure = solver.physics.compute_additional_variable("Pressure", solver.state_coeffs, True)  

    print(f"The density at t=0, x=0 is: {rho.ravel()[0]} while the pressure: {pressure.ravel()[0]/1e6} MPa")

def plot_pressure_vs_density(folder, frame=100):
    solver = readwritedatafiles.read_data_file(f"{folder}/test_output_{frame}.pkl") 
    rho = solver.state_coeffs[..., 0:3].sum(axis=-1, keepdims=True)
    pressure = solver.physics.compute_additional_variable("Pressure", solver.state_coeffs, True) 



    plt.plot(rho.ravel(), pressure.ravel()/1e6)
    plt.xlabel("density [kg/m^3]")
    plt.ylabel("pressure [MPa]")
    plt.show()

def animate_conduit_pressure(folder, iterations=100):
    fig = plt.figure(figsize=(8,11))
    ax = fig.add_subplot(411,autoscale_on=False,\
                            xlim=(0,-1000),ylim=(0,3))
    
    ax2 = fig.add_subplot(412,autoscale_on=False,\
                            xlim=(0,-1000),ylim=(-0.5,1.5))
    
    ax3 = fig.add_subplot(413, autoscale_on=False,\
                            xlim=(0,-1000), ylim=(0,300))


    ax4 = fig.add_subplot(414, autoscale_on=False,\
                            xlim=(0,-1000), ylim=(0,1))

    pressure_line,  = ax.plot([], [], color="blue", label="pressure")
    velocity_line, = ax2.plot([], [], color="red", label="velocity")
    sound_speed_line, = ax3.plot([], [], color="green", label="speed of sound")
    viscosity_line, = ax4.plot([], [], color="orange", label="viscosity")

    ax4.set_xlabel("Depth [m]")

    ax.set_ylabel("Pressure [MPa]")
    ax2.set_ylabel("Velocity [m/s]")
    ax3.set_ylabel("Speed of sound [m/s]")
    ax4.set_ylabel("Effective viscosity [MPa * s]")

    time_template = 'time = %.2f [s]'
    time_text = ax.text(0.7,0.9,'',transform=ax.transAxes)

    pl_template = 'P_L = %2f [M Pa]'
    pl_text = ax.text(0.7, 0.84, "", transform=ax.transAxes)

    velocity_template = 'V = %2f [m/s]'
    velocity_text = ax2.text(0.7, 0.9, "", transform=ax2.transAxes)

    def init():
        pressure_line.set_data([], [])
        velocity_line.set_data([], [])
        sound_speed_line.set_data([], [])
        viscosity_line.set_data([], [])
        time_text.set_text("")
        pl_text.set_text("")
        velocity_text.set_text("")
        return pressure_line, velocity_line, viscosity_line, time_text, pl_text, velocity_text

    def animate(i):
        solver = readwritedatafiles.read_data_file(f"{folder}/test_output_{i}.pkl")
        flag_non_physical = True
        p = solver.physics.compute_additional_variable("Pressure", solver.state_coeffs, flag_non_physical)
        v = solver.physics.compute_additional_variable("Velocity", solver.state_coeffs, flag_non_physical)
        sound_speed = solver.physics.compute_additional_variable("SoundSpeed", solver.state_coeffs, flag_non_physical)

        fsource = solver.physics.source_terms[VISCOSITY_FRICTION_SOURCE_IDX]
        viscosity = fsource.compute_viscosity(solver.state_coeffs, solver.physics)

        # Get the position of of each nodal points (location corresponding to each entry of pDensityX)
        nodal_pts = solver.basis.get_nodes(solver.order)
        # Allocate [ne] x [nb, ndims]
        x = np.empty((solver.mesh.num_elems,) + nodal_pts.shape)
        for elem_ID in range(solver.mesh.num_elems):
            # Fill coordinates in physical space
            x[elem_ID] = mesh_tools.ref_to_phys(solver.mesh, elem_ID, nodal_pts)
        
        pressure_line.set_data(x.ravel(), p.ravel()/1e6)
        velocity_line.set_data(x.ravel(), v.ravel())
        sound_speed_line.set_data(x.ravel(), sound_speed.ravel())
        viscosity_line.set_data(x.ravel(), viscosity.ravel()/1e6)
        time_text.set_text(time_template % solver.time)
        pl_text.set_text(pl_template % (p.ravel()/1e6)[0])
        velocity_text.set_text(velocity_template % (v.ravel())[0])

        return pressure_line, velocity_line, sound_speed_line, viscosity_line, time_text, pl_text, velocity_text

    plt.close()
    return animation.FuncAnimation(fig, animate, np.arange(iterations), blit=False, init_func=init, interval=50)

1. Steady state flow

• No water, no air, no crystals
• constant viscosity.
• g = 0

Conservation of momentum

$$ \rho (\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial y}) = - \frac{\partial p}{\partial y} - \rho g - \frac{2 \tau}{R} $$

The term on the left is the inertial term, and the terms on the right are pressure gradient, weight, and drag. For the sake of our rough analytical solution, we assume that the fluid is largely incompressible and as a result, we expect the steady state of have a constant velocity with respect to both time and space.

$$ 0 = - \frac{\partial p}{\partial x} - \frac{2 \tau}{R} $$

For laminar flows, $\tau = \frac{4 \mu u}{R}$ where $\mu$ is viscosity. For our steady state flow, pressure changes with respect to x linearly. $\frac{\partial p}{\partial t} = \frac{P_L - P_0}{L}$. Plugging those quantities into the above equation, we get.

$$ \begin{align} 0 &= \frac{P_L - P_0}{L} - \frac{8\mu u}{R^2} \\ u &= \frac{P_L - P_0}{L} \frac{R^2}{8 \mu} \end{align} $$

It is worth noting that because the fluid is compressible, in reality $p = p(\rho)$ and specifically pressure increases slightly with density. However, our analytical solution above should still be approximately correct given that our fluid is not that compressible.

from astropy import units as u

# the viscosity is calculated from the numerical simulation. 
mu = 4.7 * 1e5 * u.Pa * u.s
R = 50 * u.m
L = 1000 * u.m
P_0 = 0.5 * u.MPa
P_l = 0.994 * u.MPa

u_steady_state = ((P_l - P_0)/L) * (R**2 / (8*mu))

print(f"The predicted velocity of the steady state solution: {u_steady_state.to(u.m/u.s)}")

The predicted velocity of the steady state solution: 0.3284574468085106 m / s

# Still in the works 
HTML(animate_conduit_pressure("steady_state_flow", 100).to_html5_video())

Water content

I was unable to reduce "arhoWv": 0.01, # Mass water in fluid state per mixture volume without running into challenges with the numerical stability.

If I decrease the water fraction down to 0.001 I get a non-physical error:

Traceback (most recent call last):
  File "/Users/paxton/git/quail_volcano/src/quail", line 741, in <module>
    main(sys.argv[1:])
    ~~~~^^^^^^^^^^^^^^
  File "/Users/paxton/git/quail_volcano/src/quail", line 688, in main
    domain_list[0].solve()
    ~~~~~~~~~~~~~~~~~~~~^^
  File "/Users/paxton/git/quail_volcano/src/multidomain.py", line 105, in solve
    self.solver.solve()
    ~~~~~~~~~~~~~~~~~^^
  File "/Users/paxton/git/quail_volcano/src/solver/base.py", line 671, in solve
    res = stepper.take_time_step(self)
  File "/Users/paxton/git/quail_volcano/src/numerics/timestepping/stepper.py", line 249, in take_time_step
    solver.apply_limiter(U)
    ~~~~~~~~~~~~~~~~~~~~^^^
  File "/Users/paxton/git/quail_volcano/src/solver/base.py", line 571, in apply_limiter
    limiter.limit_solution(self, U)
    ~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^
  File "/Users/paxton/git/quail_volcano/src/numerics/limiting/positivitypreserving.py", line 420, in limit_solution
    raise errors.NotPhysicalError
errors.NotPhysicalError

arhoWt is currently set to 0, but modifying that value does not appear to affect my ability to decrease arhoWv.

plot_pressure_vs_density("steady_state_flow")