Skip to content

generate.lithostatic

Lithostatic stress field simulation.

This module provides functions for simulating the photoelastic response of a material under lithostatic stress conditions (stress increasing linearly with depth).

lithostatic

Functions

lithostatic_stress_cartesian(X, Y, rho, g=9.81, K0=0.5)

Lithostatic stress field (increasing with depth).

Parameters:

Name Type Description Default
X array - like

X coordinates (horizontal position).

 required
Y array - like

Y coordinates (depth, positive downward from surface).

 required
rho float

Density (kg/m^3).

 required
g float

Gravity (m/s^2).

 9.81
K0 float

Coefficient of lateral earth pressure (sigma_xx / sigma_yy).

 0.5

Returns:

Name Type Description
sigma_xx array - like

Normal stress in x direction (Pa).
sigma_yy array - like

Normal stress in y direction (Pa).
tau_xy array - like

Shear stress (Pa).
Source code in photoelastimetry/generate/lithostatic.py 
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
def lithostatic_stress_cartesian(X, Y, rho, g=9.81, K0=0.5):
    """
    Lithostatic stress field (increasing with depth).

    Parameters
    ----------
    X : array-like
        X coordinates (horizontal position).
    Y : array-like
        Y coordinates (depth, positive downward from surface).
    rho : float
        Density (kg/m^3).
    g : float
        Gravity (m/s^2).
    K0 : float
        Coefficient of lateral earth pressure (sigma_xx / sigma_yy).

    Returns
    -------
    sigma_xx : array-like
        Normal stress in x direction (Pa).
    sigma_yy : array-like
        Normal stress in y direction (Pa).
    tau_xy : array-like
        Shear stress (Pa).
    """
    # Vertical stress increases linearly with depth
    sigma_yy = rho * g * Y

    # Ensure no negative depth (if Y coordinate system differs) produces negative stress
    # Assuming soil/rock/granular material doesn't support tension this way usually
    sigma_yy = np.maximum(0, sigma_yy)

    # Horizontal stress is a fraction of vertical stress
    sigma_xx = K0 * sigma_yy

    # No shear stress in simple lithostatic case
    tau_xy = np.zeros_like(X)

    return sigma_xx, sigma_yy, tau_xy

generate_synthetic_lithostatic(X, Y, rho, g, K0, S_i_hat, wavelengths_nm, thickness, C)

Generate synthetic lithostatic stress data.

Source code in photoelastimetry/generate/lithostatic.py 
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
def generate_synthetic_lithostatic(X, Y, rho, g, K0, S_i_hat, wavelengths_nm, thickness, C):
    """
    Generate synthetic lithostatic stress data.
    """

    # Get stress components
    sigma_xx, sigma_yy, tau_xy = lithostatic_stress_cartesian(X, Y, rho, g, K0)

    # Principal stress difference and angle
    principal_diff = np.abs(sigma_xx - sigma_yy)
    theta_p = 0.5 * np.arctan2(2 * tau_xy, sigma_xx - sigma_yy)

    height, width = sigma_xx.shape
    n_colors = len(wavelengths_nm)

    synthetic_images = np.zeros((height, width, n_colors, 4))

    for i, lambda_light in enumerate(wavelengths_nm):
        # Stress-optic coefficient for this wavelength
        c_lambda = C[i]

        # Simulate polarimetry
        I0, I45, I90, I135 = simulate_four_step_polarimetry(
            sigma_xx, sigma_yy, tau_xy, c_lambda, 1.0, thickness, lambda_light, S_i_hat  # nu (solid fraction)
        )

        # Stack into our array
        synthetic_images[:, :, i, 0] = I0
        synthetic_images[:, :, i, 1] = I45
        synthetic_images[:, :, i, 2] = I90
        synthetic_images[:, :, i, 3] = I135

    return synthetic_images, principal_diff, theta_p, sigma_xx, sigma_yy, tau_xy