Modeling the Horizontal Buckling Force of a Lithospheric Plate Using Ch21.4 Techniques

In this project I will be modeling two cases:

• Case 1 : Homogeneous Equation using methods similar to HW7

• Case 2 : Homogeneous Equation using matrix determinants

\[ \begin{aligned} y''+ Dy &=0\\ y(0) &=0\\ y(L) &=0\\ D &= 3.3 \times10^{20} \;\textrm{Nm} \end{aligned} \]

For correctly modeling the buckling of a lithospheric plate, I am working under the following assumptions:

1. The density of the rock is \( 2360\frac{kg}{m^3} \) at a depth of \( 5000m \) below mean sea level

2. \[ \begin{align*} P& = 120 \; \textrm{MPa}\\ D& = 3.3 \times10^{20} \;\textrm{Nm} \end{align*} \]

3. I have assumed the length of the plate is 20 meters long and the ends of the plate are fixed

I have sourced these above assumptions for 1 and 2 from Otsubo, M., Katayama, I., Miyakawa, A. et al., 2020.

Side note: \( 1\textrm{Pa}=\frac{1\textrm{N}}{\textrm{m}^2} \)

Using the IVP from the introduction slide, a difference equation found by a Taylor expansions will numerically solve our IVP. The Taylor expansion is shown below:

\[ y^{''}=\frac{1}{h^2}[y(x+h)-2y(x)+y(x-h)] \]

Simplifying the results and applying it to the case that models the buckling of a plate:

\[ \begin{aligned} 0& = (\frac{1}{h^2}[y(x+h)-2y(x)+y(x-h)])*h^2\\ 0& = [y_{i+1}-2y_i+y_{i-1}]+h^2Dy_i\\ 0& = -y_{i+1}+y_i(2-Dh^2)-y_{i-1}\\ a& = 2-Dh^2 \end{aligned} \]

\[ \begin{array} {} i&=1:& u_0 & +& au_1 & -& u_2 & +& 0u_3 & + & 0u_4 +& 0u_5 &= 0\\ i&=2:& u_0 & -& u_1 & +& au_2 & -& u_3 & + & 0u_4 +& 0u_5 &= 0\\ i&=3:& u_0 & +& 0u_1 & -& u_2 & +& au_3 & - & u_4 +& 0u_5 &= 0\\ i&=4:& u_0 & +& 0u_1 & -& 0u_2 & -& u_3 & + & au_4 -& u_5 &= 0\\ \end{array} {} \]

\[ \begin{array} {} i&=1:& 0 & +& au_1 & -& u_2 & +& 0u_3 & + & 0u_4 +& 0 &= 0\\ i&=2:& 0 & -& u_1 & +& au_2 & -& u_3 & + & 0u_4 +& 0 &= 0\\ i&=3:& 0 & +& 0u_1 & -& u_2 & +& au_3 & - & u_4 +& 0 &= 0\\ i&=4:& 0 & +& 0u_1 & -& 0u_2 & -& u_3 & + & au_4 -& 0 &= 0\\ \end{array} {} \]

\[ A = \begin{bmatrix} a & -1 & 0 & 0 \\ -1 & a & -1 & 0 \\ 0 & -1 & a & -1 \\ 0 & 0 & -1 & -a \end{bmatrix} \, \begin{bmatrix} u_1 \\ u_2 \\ u_3 \\ u_4 \end{bmatrix} \ = \begin{bmatrix} 0\\ 0\\ 0\\ 0 \end{bmatrix} \]

disp('Enter in A and b:')
D=120;          #plate rigidity
h=10-2/2;       #step size
a=2-D*(h^2);
A = [a,-1,0,0;-1,a,-1,0;0,-1,a,-1;0,0,-1,-a]  #Initialize matrix A
b=[0;0;0;0]                #Initialize vector b
disp('Use octave to find inverse G of A, using G=inv(A) command:')
G=inv(A)                           #Finds the inverse of matrix A
disp('Solve Au=b by computing u=G*b:')
u=G*b

A =

  -9718     -1      0      0
     -1  -9718     -1      0
      0     -1  -9718     -1
      0      0     -1   9718

b =

   0
   0
   0
   0

u =

   0
   0
   0
   0

The graph below interpolates the displacement at each of the four mesh points within the lithospheric plate.

N = 10;      #Nodes
L = 20;       #Plate Length
M = N+2;      #M = number of boundary nodes plus interior nodes
h = L/(N+1);  #step size along the plate 
h2= h^2;      #from notes, needed for calculation
p = pi/L;     #First eigen value for analytical solution
a = 2-((h^2)*(p^2));

y = zeros(M,1);     #initialize y-values
t = linspace(0,L,M);  #initialize t-values

y(1)=sin(p*t(1));
y(2)=sin(p*t(2));

for k= L:M-1;         #for-loop for critical buckling force
  y(k)=a*y(k-1)-y(k-2);
end

A= zeros(N,N);        #Set up matrix
A(1,1)=a;
A(1,2)=-1;
for i=2:N-1
  A(i,i-1:i+1)=[-1,a,-1];     #For-loop for tri-diagonal matrix
end 
A(N,N-1)=-1;
A(N,N)=a;
A;
det = det(A)

eigenvalues = eig(A);
e1 = eigenvalues(1)

#To find eigenvector u, solve (A-cI)u = 0
I = eye(N);
B = A-e1*1
C = rref(B)

C(N,N) = 1;
d = [zeros(1,N-1), sin(p*t(M-1))];
u=C\d;

figure

subplot(1,2,1)
plot(t,y)

subplot(1,2,2)
t1=linspace(0,L,N);
plot(t1,u)

print("Case2Plot.png")

Displays the plot that interpolates the displacement of a lithospheric plate. The first eigen value is the critical buckling force.

det = -0.038122
e1 = -5.5293e-04

• Case 1 is displaying the correct values for when our \( \mathbf{b} \) has all zero values

• Case 2 is displaying a correct displacement graph and the determinant values (not equaling zero) were accounted for using an eigenvector

• Increasing the nodes in Case 2 models the displacement graph appropriately, the curve is sinusoidal

• Modeling lithospheric flexure without fixed ends

• Apply actual field data from Grand Junction into this numerical method to test accuracy

• Add more comments to code

• HW7: Ch21.4

• Class notes: Ch21.4

• Textbook: Ch21.4

• Otsubo, M., Katayama, I., Miyakawa, A. et al. Inelastic behavior and mechanical strength of the shallow upper crust controlled by layer-parallel slip in the high-strain zone of the Niigata region, Japan. Earth Planets Space 72, 30 (2020). https://doi.org/10.1186/s40623-020-01154-w

• Schubert, G., Turcotte, D., 2002, Geodynamics: Cambridge University Press, p. 211-212