Elasticity solution for wedge-shaped problems
======================================

In this notebook an analytic solution for a wedge is constructed. Later on, we couple two wedge solutions to compute the displacement for a dipping fault with free surface.

Wedge solution
----------------------

Polar coordinates are given by
\begin{align}
    x &= r\cos(\theta) \\
    y &= r\sin(\theta) \\
    r &= \sqrt{x^2+y^2} \\
    \theta &= \arctan(y/x) \\
\end{align}

For a plane-strain problem, we assume there exists a biharmonic function $\phi$ such that [1]

\begin{align}
    2\mu u_r &= -\frac{\partial \phi}{\partial r} + (1-\nu)r\frac{\partial\psi}{\partial \theta} \\
    2\mu u_\theta &= -\frac{1}{r}\frac{\partial \phi}{\partial \theta} + (1-\nu)r^2\frac{\partial\psi}{\partial r} \\
    \sigma_{rr} &= \frac{1}{r}\frac{\partial \phi}{\partial r} +
    \frac{1}{r^2}\frac{\partial^2 \phi}{\partial \theta^2} \\
    \sigma_{\theta\theta} &= \frac{\partial^2 \phi}{\partial r^2} \\
    \sigma_{r\theta} &= -\frac{\partial }{\partial r} \left(
      \frac{1}{r}\frac{\partial \phi}{\partial\theta}
    \right) \\
\end{align}
where $\mu$ is the shear modulus, $\nu$ is Poisson's ratio, and $\phi$ and $\psi$ are related by
$$
\nabla^2\phi = \frac{\partial}{\partial r}\left(r\frac{\partial\psi}{\partial\theta}\right)
$$

Remark: Plane-stress is modelled in [1], therefore the formula are adjusted here for plane-strain.

[1] "Stress Singularites Resulting From Various Boundary Conditions in Angular Corners of Plates in Extension", M. L. Williams, Journal of Applied Mechanics, 19 (4), pp. 526-528.


```python
from sympy import *
init_printing()
```

Assume

$$
\phi = r^{\gamma+1} \left(A_1\sin((\gamma+1)\theta) + A_2\cos((\gamma+1)\theta) + A_3\sin((\gamma-1)\theta) + A_4\sin((\gamma-1)\theta)\right)
$$
and
$$
\psi = r^{\gamma-1}\left(-\frac{4A_3}{\gamma-1}\cos((\gamma-1)t) + \frac{4A_4}{\gamma-1}\sin((\gamma-1)t)\right)
$$


```python
r, t, x, y, gam = symbols('r theta x y gamma')
mu, nu, m, A0, A1, A2, A3, A4, S = symbols('mu nu m A_0 A_1 A_2 A_3 A_4 S')

phi = r**(gam+1) * (A1 * sin((gam+1)*t) + A2*cos((gam+1)*t) + A3*sin((gam-1)*t) + A4*cos((gam-1)*t))
psi = r**(gam-1) * (-4*A3/(gam-1)*cos((gam-1)*t) + 4*A4/(gam-1)*sin((gam-1)*t))

phi2 = phi.subs(gam, 2)
psi2 = psi.subs(gam, 2)
(phi2, psi2)
```

$\psi$ must be harmonic, i.e. $\nabla^2\psi=0$, and $\phi$ must be biharmonic, i.e. $\nabla^2(\nabla^2\phi=0)$:


```python
def laplace(f):
    return f.diff(r, 2) + f.diff(r)/r + f.diff(t, 2)/r**2
(laplace(psi).simplify(), laplace(laplace(phi)).simplify())
```

The relation of $\phi$ and $\psi$ must be satisfied:


```python
((r*psi.diff(t)).diff(r) - laplace(phi)).simplify()
```

Define displacement and stress:


```python
def quants(phi, psi):
    ur = (-phi.diff(r) + (1-nu)*r*psi.diff(t)) / (2 * mu)
    ut = (-phi.diff(t)/r + (1-nu)*r**2*psi.diff(r)) / (2 * mu)
    srr = (phi.diff(r) / r + phi.diff(t, 2) / r**2)
    stt = phi.diff(r, 2)
    srt = -(phi.diff(t)/r).diff(r)
    return (ur.simplify(), ut.simplify(), srr.simplify(), stt.simplify(), srt.simplify())
ur, ut, srr, stt, srt = quants(phi, psi)
```

Check that stress-displacement relation is satisfied:


```python
err = ur.diff(r)
ett = ut.diff(t)/r + ur/r
ert = (ur.diff(t)/r + ut.diff(r) - ut/r) / 2
((err - ((1-nu)*srr - nu*stt)/(2*mu)).simplify(),
 (ett - ((1-nu)*stt - nu*srr)/(2*mu)).simplify(),
 (ert - srt/(2*mu)).simplify())
```

The planes $\theta=0$ and $\theta=-\pi/3$ form an infinite wedge.
We set $\gamma=2$ and recover the coefficients $A_1,A_2,A_3,A_4$ which satisfy
\begin{align}
    u_\theta(r,-\pi/3) &= 0 \\
    u_r(r,-\pi/3) &= Sr^2 \\
    \sigma_{r\theta}(r,0) &= 0 \\
    \sigma_{\theta\theta}(r,0) &= 0
\end{align}


```python
w1 = 0
w2 = -pi/3
mygam = 2
solution = linsolve([ut.subs(gam, mygam).subs(t, w2)/r**2,
                     ur.subs(gam, mygam).subs(t, w2)/r**2 - S,
                     srt.subs(gam, mygam).subs(t, w1)/r,
                     stt.subs(gam, mygam).subs(t, w1)/r],
         A1, A2, A3, A4)
a1, a2, a3, a4 = solution.args[0]
a1, a2, a3, a4 = a1.simplify(), a2.simplify(), a3.simplify(), a4.simplify()
(a1, a2, a3, a4)
```

We check that the boundary conditions are satisfied:


```python
srr_s = srr.subs(gam, mygam).subs(A1, a1).subs(A2, a2).subs(A3, a3).subs(A4, a4).simplify()
stt_s = stt.subs(gam, mygam).subs(A1, a1).subs(A2, a2).subs(A3, a3).subs(A4, a4).simplify()
srt_s = srt.subs(gam, mygam).subs(A1, a1).subs(A2, a2).subs(A3, a3).subs(A4, a4).simplify()
(stt_s.subs(t, w1), srt_s.subs(t, w1))
```


```python
ur_s = ur.subs(gam, mygam).subs(A1, a1).subs(A2, a2).subs(A3, a3).subs(A4, a4).simplify()
ut_s = ut.subs(gam, mygam).subs(A1, a1).subs(A2, a2).subs(A3, a3).subs(A4, a4).simplify()
(ur_s.subs(t, w2).simplify(), ut_s.subs(t, w2).simplify())
```

We convert the solution into Lua code:


```python
def lua_code(expr):
    return rcode(expr).replace('cos', 'math.cos').replace('sin', 'math.sin').replace('sqrt', 'math.sqrt').replace('pi', 'math.pi')
print("local A_1 =", lua_code(a1))
print("local A_2 =", lua_code(a2))
print("local A_3 =", lua_code(a3))
print("local A_4 =", lua_code(a4))
print("local ur =", lua_code(ur.subs(gam, mygam).simplify()))
print("local ut =", lua_code(ut.subs(gam, mygam).simplify()))
```

    local A_1 = (1.0/3.0)*math.sqrt(3)*S*mu*(4*nu - 5)/(-16*nu^2 + 20*nu + 5)
    local A_2 = S*mu*(7 - 4*nu)/(-16*nu^2 + 20*nu + 5)
    local A_3 = math.sqrt(3)*S*mu*(5 - 4*nu)/(-16*nu^2 + 20*nu + 5)
    local A_4 = S*mu*(4*nu - 7)/(-16*nu^2 + 20*nu + 5)
    local ur = (1.0/2.0)*r^2*(-3*A_1*math.sin(3*theta) - 3*A_2*math.cos(3*theta) - 4*A_3*nu*math.sin(theta) + A_3*math.sin(theta) - 4*A_4*nu*math.cos(theta) + A_4*math.cos(theta))/mu
    local ut = (1.0/2.0)*r^2*(-3*A_1*math.cos(3*theta) + 3*A_2*math.sin(3*theta) - A_3*math.cos(theta) + A_4*math.sin(theta) + (4*nu - 4)*(A_3*math.cos(theta) - A_4*math.sin(theta)))/mu


Dipping fault
------------------

We now couple two wedge solutions to impose a displacement in a dipping fault with a free-surface condition, as illustrated in the following figure:



where $[\![ u ]\!] = u_A - u_B$. The A-side and B-side of the fault are separated by the fault plane $\theta=\omega$. Also note that the origin, marked with a circle, is at the intesection of the fault and the free surface.

We introduce potentials for the A-side and the B-side of the fault:


```python
B1, B2, B3, B4 = symbols('B_1 B_2 B_3 B_4')

phiA = phi.simplify()
psiA = psi.simplify()
phiB = phiA.subs(A1, B1).subs(A2, B2).subs(A3, B3).subs(A4, B4).simplify()
psiB = psiA.subs(A1, B1).subs(A2, B2).subs(A3, B3).subs(A4, B4).simplify()
(phiA, psiA, phiB, psiB)
```


```python
urA, utA, srrA, sttA, srtA = quants(phiA, psiA)
urB, utB, srrB, sttB, srtB = quants(phiB, psiB)
```

We remove the dependency on $r$ here to solve for the A and B coefficients:


```python
sttA_e = (sttA*r**(1-gam)/gam).simplify()
srtA_e = (srtA*r**(1-gam)/gam).simplify()
sttB_e = (sttB*r**(1-gam)/gam).simplify()
srtB_e = (srtB*r**(1-gam)/gam).simplify()
(sttA_e, srtA_e, sttB_e, srtB_e)
```


```python
def simp_A(expr):
    return factor(expr, *parsA).simplify()
def simp_B(expr):
    return factor(expr, *parsB).simplify()

parsA = (A1, A2, A3, A4)
urA_e = simp_A(urA*r**(-gam)*2*mu)
utA_e = simp_A(utA*r**(-gam)*2*mu)
parsB = (B1, B2, B3, B4)
urB_e = simp_B(urB*r**(-gam)*2*mu)
utB_e = simp_B(utB*r**(-gam)*2*mu)
(urA_e, utA_e, urB_e, utB_e)
```

We first compute coefficients which ensure the free-surface boundary condition. These calculations are independent for the A-side and B-side. Note that the A-wedge is given by planes $\theta=\pi$ and $\theta=\omega$, and the B-wedge is given by planes $\theta=0$ and $\theta=\omega$.


```python
part_solution_A = solve([sttA_e.subs(t, pi),
                         srtA_e.subs(t, pi)],
                         *parsA)
part_solution_A = {k: v.simplify() for k, v in part_solution_A.items()}
part_solution_B = solve([sttB_e.subs(t, 0),
                         srtB_e.subs(t, 0)],
                         *parsB)
part_solution_B = {k: v.simplify() for k, v in part_solution_B.items()}
(part_solution_A, part_solution_B)
```

We can now express $A_1,A_2,B_1,B_2$ in terms of $A_3,A_4,B_3,B_3$, which we substitute in the following:


```python
sttA_e2 = simp_A(sttA_e.subs(part_solution_A.items()))
srtA_e2 = simp_A(srtA_e.subs(part_solution_A.items()))
sttB_e2 = simp_B(sttB_e.subs(part_solution_B.items()))
srtB_e2 = simp_B(srtB_e.subs(part_solution_B.items()))
(sttA_e2, srtA_e2, sttB_e2, srtB_e2)
```


```python
utA_e2 = simp_A(utA_e.subs(part_solution_A.items()))
urA_e2 = simp_A(urA_e.subs(part_solution_A.items()))
urB_e2 = simp_B(urB_e.subs(part_solution_B.items()))
utB_e2 = simp_B(utB_e.subs(part_solution_B.items()))
(urA_e2, utA_e2, urB_e2, utB_e2)
```

We now find the remaining coefficients for the imposed displacement on the fault plane $\omega = -\pi/3$. Note that we require the stress tensor to be continuous.


```python
w = -pi/3
equations = [(sttA_e2 - sttB_e2).subs(t, w).simplify(),
             (srtA_e2 - srtB_e2).subs(t, w).simplify(),
             (utA_e2 - utB_e2).subs(t, w).simplify(),
             (urA_e2 - urB_e2 - S).subs(t, w).simplify()]
equations
```


```python
part_equations_0 = [factor(equations[0] + equations[3], *parsA, *parsB).simplify(),
                    factor(equations[1] - equations[2], *parsA, *parsB).simplify()]
part_equations_0
```


```python
part_solution_0 = solve(part_equations_0, A3, A4, B3, B4)
part_solution_0 = {k: v.simplify() for k, v in part_solution_0.items()}
part_solution_0
```


```python
part_equations_1 = [factor(equations[0].subs(part_solution_0.items()).simplify(), *parsB).simplify(),
                    factor(equations[1].subs(part_solution_0.items()).simplify(), *parsB).simplify()]
```


```python
solution = solve(part_equations_1, B3, B4)
solution = {k: v.simplify() for k, v in solution.items()}
solution
```

Note that $\sin(\pi\gamma)$ appears in the denominator above.
Thus, $\gamma$ must not be an integer.

We substitute the parameters into the quantities for $\gamma = 0.5$:


```python
mygam = 0.5
def subs_solution_A(expr):
    return expr.subs(part_solution_A.items()).subs(part_solution_0.items()).subs(solution.items()).subs(gam, mygam)
def subs_solution_B(expr):
    return expr.subs(part_solution_B.items()).subs(part_solution_0.items()).subs(solution.items()).subs(gam, mygam)
```


```python
srrA_s = subs_solution_A(srrA).simplify()
sttA_s = subs_solution_A(sttA).simplify()
srtA_s = subs_solution_A(srtA).simplify()
urA_s = subs_solution_A(urA).simplify()
utA_s = subs_solution_A(utA).simplify()

srrB_s = subs_solution_B(srrB).simplify()
sttB_s = subs_solution_B(sttB).simplify()
srtB_s = subs_solution_B(srtB).simplify()
urB_s = subs_solution_B(urB).simplify()
utB_s = subs_solution_B(utB).simplify()
```

In the following we test whether the boundary conditions are satisifed.


```python
(sttA_s.subs(t, pi).simplify(),
 srtA_s.subs(t, pi).simplify(),
 sttB_s.subs(t, 0).simplify(),
 srtB_s.subs(t, 0).simplify())
```


```python
((sttA_s.subs(t, w)-sttB_s.subs(t, w)).simplify().evalf(),
 (srtA_s.subs(t, w)-srtB_s.subs(t, w)).simplify().evalf(),
 (utA_s.subs(t, w)-utB_s.subs(t, w)).simplify().evalf())
```


```python
((urA_s.subs(t, w) - urB_s.subs(t, w) - S).simplify() + S).evalf().simplify()
```

Lastly, we generate Lua code for the A-side and the B-side:


```python
def lua_code_cse(var_expr_pairs):
    vs, es = zip(*var_expr_pairs)
    e_cse = cse(es)
    for k, v in e_cse[0]:
        print("local {} = {}".format(k, lua_code(v)))
    for k, v in zip(vs, e_cse[1]):
        print("local {} = {}".format(k, lua_code(v)))
```

$B_3$ and $B_4$


```python
lua_code_cse(solution.items())
```

    local x0 = math.pi*((1.0/3.0)*gamma + 1.0/3.0)
    local x1 = math.sin(x0)
    local x2 = 2*x1
    local x3 = math.pi*((7.0/3.0)*gamma + 1.0/3.0)
    local x4 = math.sin(x3)
    local x5 = 2*x4
    local x6 = 3*gamma
    local x7 = math.cos(x0)
    local x8 = math.sqrt(3)*gamma
    local x9 = math.cos(x3)
    local x10 = math.pi*gamma
    local x11 = (8.0/3.0)*x10
    local x12 = math.sin(x11)
    local x13 = 2*x7
    local x14 = 2*x9
    local x15 = math.cos(x11)
    local x16 = (1.0/32.0)*S/((nu - 1)*math.sin(x10)^2)
    local B_3 = x16*(x1*x6 + x12*x13 - x12*x14 - x15*x2 + x15*x5 + x2 - x4*x6 - x5 + x7*x8 - x8*x9)
    local B_4 = x16*(-x1*x8 - x12*x2 + x12*x5 - x13*x15 + x13 + x14*x15 - x14 + x4*x8 + x6*x7 - x6*x9)


$A_3$ and $A_4$


```python
lua_code_cse(part_solution_0.items())
```

    local x0 = nu - 1
    local x1 = 1.0/x0
    local x2 = math.pi*((1.0/3.0)*gamma + 1.0/6.0)
    local x3 = (1.0/4.0)*S
    local A_3 = x1*(B_3*x0 - x3*math.cos(x2))
    local A_4 = x1*(B_4*x0 - x3*math.sin(x2))


$A_1$ and $A_2$


```python
lua_code_cse(part_solution_A.items())
```

    local x0 = 1.0/(gamma + 1)
    local x1 = 2*math.pi*gamma
    local x2 = math.cos(x1)
    local x3 = math.sin(x1)
    local A_1 = x0*(-A_3*gamma + A_3*x2 - A_4*x3)
    local A_2 = -x0*(A_3*x3 + A_4*gamma + A_4*x2)


$B_1$ and $B_2$


```python
lua_code_cse(part_solution_B.items())
```

    local B_2 = -B_4
    local B_1 = B_3*(1 - gamma)/(gamma + 1)



```python
urA
```

A-side: $u_r$ and $u_\theta$


```python
lua_code_cse([('urA', urA), ('utA', utA)])
```

    local x0 = gamma - 1
    local x1 = theta*x0
    local x2 = math.sin(x1)
    local x3 = math.cos(x1)
    local x4 = A_3*x2 + A_4*x3
    local x5 = gamma + 1
    local x6 = theta*x5
    local x7 = math.sin(x6)
    local x8 = math.cos(x6)
    local x9 = (1.0/2.0)*r^gamma/mu
    local x10 = A_4*x2
    local x11 = A_3*x3
    local urA = x9*(4*x4*(1 - nu) - x5*(A_1*x7 + A_2*x8 + x4))
    local utA = x9*(-A_1*x5*x8 + A_2*x5*x7 + x0*x10 - x0*x11 + (4*nu - 4)*(-x10 + x11))


B-side: $u_r$ and $u_\theta$


```python
lua_code_cse([('urB', urB), ('utB', utB)])
```

    local x0 = gamma - 1
    local x1 = theta*x0
    local x2 = math.sin(x1)
    local x3 = math.cos(x1)
    local x4 = B_3*x2 + B_4*x3
    local x5 = gamma + 1
    local x6 = theta*x5
    local x7 = math.sin(x6)
    local x8 = math.cos(x6)
    local x9 = (1.0/2.0)*r^gamma/mu
    local x10 = B_4*x2
    local x11 = B_3*x3
    local urB = x9*(4*x4*(1 - nu) - x5*(B_1*x7 + B_2*x8 + x4))
    local utB = x9*(-B_1*x5*x8 + B_2*x5*x7 + x0*x10 - x0*x11 + (4*nu - 4)*(-x10 + x11))


A-side: $u_{r,r}, u_{r,\theta}, u_{\theta,r}, u_{\theta,\theta}$


```python
lua_code_cse([('urA_r', urA.diff(r).simplify()),
              ('urA_t', urA.diff(t).simplify()),
              ('utA_r', utA.diff(r).simplify()),
              ('utA_t', utA.diff(t).simplify())])
```

    local x0 = 4*nu - 4
    local x1 = gamma - 1
    local x2 = theta*x1
    local x3 = math.sin(x2)
    local x4 = A_3*x3
    local x5 = math.cos(x2)
    local x6 = A_4*x5
    local x7 = x4 + x6
    local x8 = gamma + 1
    local x9 = theta*x8
    local x10 = math.sin(x9)
    local x11 = A_1*x10
    local x12 = math.cos(x9)
    local x13 = A_2*x12
    local x14 = (1.0/2.0)/mu
    local x15 = gamma*r^x1*x14
    local x16 = A_3*x5
    local x17 = A_4*x3
    local x18 = x16 - x17
    local x19 = 4*x1*(nu - 1)
    local x20 = A_1*x12*x8
    local x21 = x1*x16
    local x22 = A_2*x10*x8
    local x23 = x1*x17
    local x24 = r^gamma*x14
    local x25 = x8^2
    local x26 = x1^2
    local urA_r = -x15*(x0*x7 + x8*(x11 + x13 + x7))
    local urA_t = -x24*(x18*x19 + x8*(x20 + x21 - x22 - x23))
    local utA_r = x15*(x0*x18 - x20 - x21 + x22 + x23)
    local utA_t = x24*(x11*x25 + x13*x25 - x19*x7 + x26*x4 + x26*x6)


B-side: $u_{r,r}, u_{r,\theta}, u_{\theta,r}, u_{\theta,\theta}$


```python
lua_code_cse([('urB_r', urB.diff(r).simplify()),
              ('urB_t', urB.diff(t).simplify()),
              ('utB_r', utB.diff(r).simplify()),
              ('utB_t', utB.diff(t).simplify())])
```

    local x0 = 4*nu - 4
    local x1 = gamma - 1
    local x2 = theta*x1
    local x3 = math.sin(x2)
    local x4 = B_3*x3
    local x5 = math.cos(x2)
    local x6 = B_4*x5
    local x7 = x4 + x6
    local x8 = gamma + 1
    local x9 = theta*x8
    local x10 = math.sin(x9)
    local x11 = B_1*x10
    local x12 = math.cos(x9)
    local x13 = B_2*x12
    local x14 = (1.0/2.0)/mu
    local x15 = gamma*r^x1*x14
    local x16 = B_3*x5
    local x17 = B_4*x3
    local x18 = x16 - x17
    local x19 = 4*x1*(nu - 1)
    local x20 = B_1*x12*x8
    local x21 = x1*x16
    local x22 = B_2*x10*x8
    local x23 = x1*x17
    local x24 = r^gamma*x14
    local x25 = x8^2
    local x26 = x1^2
    local urB_r = -x15*(x0*x7 + x8*(x11 + x13 + x7))
    local urB_t = -x24*(x18*x19 + x8*(x20 + x21 - x22 - x23))
    local utB_r = x15*(x0*x18 - x20 - x21 + x22 + x23)
    local utB_t = x24*(x11*x25 + x13*x25 - x19*x7 + x26*x4 + x26*x6)



```python

```

f579a00c5d5d2f87edfb396fc236f9174a07b72d
ipynb
Jupyter Notebook
notebooks/wedge.ipynb
TEAR-ERC/tandem
2db44f324c3205e81d28633ede261ad179334aa3
BSD-3-Clause
2021-11-03T17:11:00.000Z
2022-02-16T07:51:01.000Z
notebooks/wedge.ipynb
TEAR-ERC/tandem
2db44f324c3205e81d28633ede261ad179334aa3
BSD-3-Clause
2020-05-18T14:51:13.000Z
2021-12-16T12:56:31.000Z
notebooks/wedge.ipynb
TEAR-ERC/tandem
2db44f324c3205e81d28633ede261ad179334aa3
BSD-3-Clause
2020-10-23T08:04:42.000Z
2021-11-15T12:23:59.000Z
Qwen/Qwen-72B
1. YES
2. YES