2025.1.25 Experiments¶

In [92]:
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)
In [93]:
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("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)

    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("")
        return pressure_line, velocity_line, viscosity_line, time_text, pl_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])

        return pressure_line, velocity_line, sound_speed_line, viscosity_line, time_text, pl_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.

In [94]:
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.24 * u.MPa
P_l = 1.4 * 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.7712765957446808 m / s

The predicted velocity is 25% less than the velocity we ultimately measure numerically.

In [95]:
# Still in the works 
HTML(animate_conduit_pressure("steady_state_flow", 100).to_html5_video())
Out[95]:
Your browser does not support the video tag.

Visocosity¶

At the moment the viscosity evolves with time and appears to drop close to zero starting from the x1 boundary.

I am currently defining the boundary condition as:

    'x1': {
        #'BCType': "SlipWall",
        'BCType': 'VelocityInlet1D',   # Inlet boundary condition
        'u': 1,
        'yWt': 0,
        'yC': 0,
        'trace_arho': 0,
    },

I wonder if something about the Inlet boundary condition is allowing some substance with a very low viscosity into the system. That said, it appears to me that for the material properties, the inlet boundary condition uses the same K and $\rho$ value that are defined in the physics section of the simulation.

Water content¶

I was unable to reduce "arhoWv": 0.008, # 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:

  File "/Users/paxton/git/quail_volcano/src/numerics/limiting/positivitypreserving.py", line 420, in limit_solution
    raise errors.NotPhysicalError
errors.NotPhysicalError
In [97]:
plot_pressure_vs_density("steady_state_flow")
No description has been provided for this image