# Circular hole probllem

In this notebook an analytic solution for a disk-shaped hole is constructed.

The domain defined as
$$
    \Omega = \{ x \in \mathbb{R^2} : 0 < a < |x| < A\}
$$
Let
$$
    \Gamma^r = \{ x \in \mathbb{R^2} : |x| = r\},
$$
then $\partial\Omega = \Gamma^a \cup \Gamma^A$.

## Poisson

Sought is the solution to the problem
$$
\begin{aligned}
    \nabla^2 u = 0, & \quad x \in \Omega, \\
    \frac{\partial u}{\partial r} = 0, & \quad x \in \Gamma^a, \\
    u = g, & \quad x \in \Gamma^A, \\
\end{aligned}
$$
where the Laplacian in polar coordinates is
$$
    \nabla^2 = \frac{\partial^2}{\partial r^2}
               + \frac{1}{r}\frac{\partial}{\partial r}
               + \frac{1}{r^2}\frac{\partial^2}{\partial \theta^2}
$$
and $g$ is chosen such that the solution is easy to construct :-)


```python
from sympy import *
from util import lua_code
```


```python
# Laplacian in polar coordinates
def laplacian(f):
    return f.diff(r, 2) + f.diff(r)/r + f.diff(t, 2)/r**2
```

Solution ansatz:


```python
r, t = symbols('r theta')
a, B = symbols('a B')
u = r**2 * sin(2*t) + B * sin(2*t) / r**2
u
```




$\displaystyle \frac{B \sin{\left(2 \theta \right)}}{r^{2}} + r^{2} \sin{\left(2 \theta \right)}$



Check Laplacian


```python
simplify(laplacian(u))
```




$\displaystyle 0$



Check inner circle boundary


```python
u.diff(r).subs(r, a)
```




$\displaystyle - \frac{2 B \sin{\left(2 \theta \right)}}{a^{3}} + 2 a \sin{\left(2 \theta \right)}$



Hence,


```python
u = u.subs(B, a**4)
u.diff(r).subs(r, a)
```




$\displaystyle 0$



We set $g=u$ on the outer circle boundary.

Gradient:


```python
u_x = simplify(u.diff(r) * cos(t) - u.diff(t) * sin(t) / r)
u_x
```




$\displaystyle - \frac{2 a^{4} \sin{\left(3 \theta \right)}}{r^{3}} + 2 r \sin{\left(\theta \right)}$




```python
u_y = simplify(u.diff(r) * sin(t) + u.diff(t) * cos(t) / r)
u_y
```




$\displaystyle \frac{2 a^{4} \cos{\left(3 \theta \right)}}{r^{3}} + 2 r \cos{\left(\theta \right)}$



## Elasticity

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, 1952, Journal of Applied Mechanics, 19 (4), pp. 526-528.

Sought is the solution to the problem
$$
\begin{aligned}
    \sigma_{rr} = 0, & \quad x \in \Gamma^a, \\
    \sigma_{r\theta} = 0, & \quad x \in \Gamma^a, \\
    u_r = g_r, & \quad x \in \Gamma^A, \\
    u_\theta = g_\theta, & \quad x \in \Gamma^A, \\
\end{aligned}
$$
for $\phi,\psi$ such that the equations of elasticity are satisfied.

Solution ansatz (from J.R. Barber, Elasticity, 2002):


```python
r, t = symbols('r theta')
a, S, mu, nu = symbols('a S mu nu')
A = S*a**2
B = -S*a**4/2
phi = -S * r**2 * sin(2*t) / 2 + A * sin(2*t) + B * sin(2*t) / r**2
phi
```




$\displaystyle - \frac{S a^{4} \sin{\left(2 \theta \right)}}{2 r^{2}} + S a^{2} \sin{\left(2 \theta \right)} - \frac{S r^{2} \sin{\left(2 \theta \right)}}{2}$



Check that phi is biharmonic.


```python
laplacian(simplify(laplacian(phi)))
```




$\displaystyle 0$



Compute stresses


```python
srr = simplify(phi.diff(r) / r + phi.diff(t, 2) / r**2)
stt = simplify(phi.diff(t, 2))
srt = simplify(-(phi.diff(t) / r).diff(r))
```

Check boundary condition on inner circle.


```python
(srr.subs(r, a), srt.subs(r, a))
```




    (0, 0)



Find $\psi$ from compatbility condition, disregarding integration constants.


```python
psi = integrate(integrate(simplify(laplacian(phi)), r) / r, t)
psi
```




$\displaystyle - \frac{2 S a^{2} \cos{\left(2 \theta \right)}}{r^{2}}$



Check compatbility condition.


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




$\displaystyle 0$



Compute displacement


```python
ur = simplify(-phi.diff(r) + (1 - nu) * r * psi.diff(t)) / (2*mu)
ut = simplify(-phi.diff(t) / r + (1 - nu) * r**2 * psi.diff(r)) / (2*mu)
```


```python
ur
```




$\displaystyle \frac{S \left(- a^{4} - 4 a^{2} r^{2} \left(\nu - 1\right) + r^{4}\right) \sin{\left(2 \theta \right)}}{2 \mu r^{3}}$




```python
lua_code(ur)
```




    '(1.0/2.0)*S*(-a^4 - 4*a^2*r^2*(nu - 1) + r^4)*math.sin(2*theta)/(mu*r^3)'




```python
ut
```




$\displaystyle \frac{S \left(a^{4} - 4 a^{2} \nu r^{2} + 2 a^{2} r^{2} + r^{4}\right) \cos{\left(2 \theta \right)}}{2 \mu r^{3}}$




```python
lua_code(ut)
```




    '(1.0/2.0)*S*(a^4 - 4*a^2*nu*r^2 + 2*a^2*r^2 + r^4)*math.cos(2*theta)/(mu*r^3)'



We set $g_r=u_r$ and $g_\theta = u_\theta$ on the outer circle boundary.

Compute Jacobian


```python
J = {}
J['ur_r'] = simplify(ur.diff(r))
J['ur_t'] = simplify(ur.diff(t))
J['ut_r'] = simplify(ut.diff(r))
J['ut_t'] = simplify(ut.diff(t))
```


```python
for key, value in J.items():
    print(f'local {key} = {lua_code(value)}')
```

    local ur_r = (1.0/2.0)*S*(3*a^4 + 4*a^2*nu*r^2 - 4*a^2*r^2 + r^4)*math.sin(2*theta)/(mu*r^4)
    local ur_t = -S*(a^4 + 4*a^2*r^2*(nu - 1) - r^4)*math.cos(2*theta)/(mu*r^3)
    local ut_r = (1.0/2.0)*S*(-3*a^4 + 4*a^2*nu*r^2 - 2*a^2*r^2 + r^4)*math.cos(2*theta)/(mu*r^4)
    local ut_t = S*(-a^4 + 4*a^2*nu*r^2 - 2*a^2*r^2 - r^4)*math.sin(2*theta)/(mu*r^3)



```python

```

5e90840637c255127fa4e5516890c9b5cddfe302
ipynb
Jupyter Notebook
notebooks/circular_hole.ipynb
TEAR-ERC/tandem
2db44f324c3205e81d28633ede261ad179334aa3
BSD-3-Clause
2021-11-03T17:11:00.000Z
2022-02-16T07:51:01.000Z
notebooks/circular_hole.ipynb
TEAR-ERC/tandem
2db44f324c3205e81d28633ede261ad179334aa3
BSD-3-Clause
2020-05-18T14:51:13.000Z
2021-12-16T12:56:31.000Z
notebooks/circular_hole.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