Commit 5c947896 authored by Linda Günther's avatar Linda Günther Committed by Lars Bilke
Browse files

explanation stress concentration and running model

parent 9a3c1c76
......@@ -89,7 +89,7 @@
},
{
"cell_type": "code",
"execution_count": 1,
"execution_count": 2,
"id": "d7b9c91c",
"metadata": {
"tags": []
......@@ -119,7 +119,7 @@
},
{
"cell_type": "code",
"execution_count": 2,
"execution_count": 3,
"id": "c30056fc",
"metadata": {
"tags": []
......@@ -147,7 +147,7 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 4,
"id": "6166f65f",
"metadata": {
"tags": []
......@@ -191,7 +191,7 @@
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 5,
"id": "9488928e",
"metadata": {
"tags": []
......@@ -235,7 +235,7 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 6,
"id": "59145ce9",
"metadata": {
"tags": []
......@@ -278,7 +278,7 @@
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": 7,
"id": "90a1c879-bda2-43c8-abbb-6f05f6afbfc5",
"metadata": {},
"outputs": [],
......@@ -288,7 +288,7 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": 8,
"id": "328f47df-1d12-4bad-acc4-7437528cb734",
"metadata": {},
"outputs": [],
......@@ -298,7 +298,7 @@
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": 9,
"id": "48ecaf5c-e7de-474e-9626-4d32e7c804c1",
"metadata": {},
"outputs": [
......@@ -352,18 +352,17 @@
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": 10,
"id": "3de96148-11b0-4777-8f56-5193bca3555b",
"metadata": {},
"outputs": [],
"source": [
"from ogs6py import ogs\n",
"import vtuIO"
"from ogs6py import ogs"
]
},
{
"cell_type": "code",
"execution_count": 10,
"execution_count": 11,
"id": "390fcbf7-8c0e-4396-888b-4589e815300f",
"metadata": {},
"outputs": [],
......@@ -373,7 +372,7 @@
},
{
"cell_type": "code",
"execution_count": 11,
"execution_count": 12,
"id": "753661c8-105f-4b94-a65f-a6ba6841833a",
"metadata": {},
"outputs": [
......@@ -382,7 +381,7 @@
"output_type": "stream",
"text": [
"OGS finished with project file ../disc_with_hole.prj.\n",
"Execution took 0.878124475479126 s\n"
"Execution took 0.8345866203308105 s\n"
]
}
],
......
%% Cell type:markdown id:4796061d-62e3-489e-ac20-aea05dae1a5d tags:
|<div style="width:330px"><img src="https://www.ufz.de/static/custom/weblayout/DefaultInternetLayout/img/logos/ufz_transparent_de_blue.png" width="300"/></div>|<div style="width:330px"><img src="https://discourse.opengeosys.org/uploads/default/original/1X/a288c27cc8f73e6830ad98b8729637a260ce3490.png" width="300"/></div>|<div style="width:330px"><img src="https://tu-freiberg.de/sites/default/files/media/freiberger-alumni-netzwerk-6127/wbm_orig_rgb_0.jpg" width="300"/></div>|
|---|---|--:|
%% Cell type:markdown id:ce763078-44cd-4729-b895-0f86fd997ee3 tags:
# Linear elasticity: disc with a hole
%% Cell type:markdown id:74db7348 tags:
In this benchmark example we consider a linear elastic small deformation problem. More specifically, a plate with a central hole that is put under tension on its top boundary is simulated. By exploiting symmetries, below we evaluate this problem just for the top right quarter of this disc.
Doing this it is important to define boundary conditions for the regarded section of the plate. For the bottom boundary ($x$-axis, $\theta=-90°$) and the left boundary ($y$-axis, $\theta=0°$) we prescribe Dirichlet boundary conditions that constrain the normal displacement along the edge to be zero. On the top, where a tensile traction is applied to the plate, Neumann boundary conditions are prescribed.
For the description of the plate and its load, some dimensions need to be defined. The quarter disc under consideration can be seen as a square with an edge length of $30\, \text{mm}$. The radius of the circular hole is $a = 4\, \text{mm}$ and the applied tension on the top boundary has a value of $\sigma = 10\, \text{MPa}$.
To fully capture and understand the behaviour of stress and strain distributions around the hole, it is necessary to also define material properties. In case of isotropic linear elasticity, the relevant parameters are the Young's modulus and Poisson's ratio. The following parameters are chosen here:
$ E = 1\,\text{MPa} \qquad \nu = 0.3$
For verification of the numerical implementation, the numerical solution of the problem will be compared to the analytical solution.
%% Cell type:markdown id:8361723f tags:
## Analytical solution
The overall stress distributions in the plate around the hole can be represented by Kirsch's Solution, which is expressed here in cylindrical coordinates. The following equations are valid for an infinitely extended plate. Since the hole is very small compared to the dimension of the plate, we can consider this specification as rendered. The parameter $\sigma$ stands for the applied tension whereas $a$ identifies the radius of the inner hole. $r$ and $\theta$ are the cylindrical coordinates and describe the distance of a point to the centre of the hole ($r$) and the angle ($\theta$) under which it is deflected with respect to the axis along which the traction is applied.
%% Cell type:markdown id:515b8575 tags:
\begin{align}
\sigma_{rr} \left( r , \theta \right) &=
\frac{\sigma}{2}
\left[ \left( 1 - \frac{a^2}{r^2} \right) + \left( 1 + 3 \frac{a^4}{r^4} - 4 \frac{a^2}{r^2} \right) \cos \left( 2 \theta \right) \right]
\\
\sigma_{\theta \theta} \left( r , \theta \right) &=
\frac{\sigma}{2}
\left[ \left( 1 + \frac{a^2}{r^2} \right) - \left( 1 + 3 \frac{a^4}{r^4} \right) \cos \left( 2 \theta \right) \right]
\\
\sigma_{r \theta} \left( r , \theta \right) &=
- \frac{\sigma}{2}
\left[ \left( 1 - 3 \frac{a^4}{r^4} + 2 \frac{a^2}{r^2} \right) \sin \left( 2 \theta \right) \right]
\end{align}
%% Cell type:markdown id:0370aa5c tags:
In the chosen OGS model, the traction is applied as mentioned before on the top boundary of the plate along the $y$-axis, so the $x$-axis is at $\theta = -90°$.
To visualise the problem, we plot the stress distribution along the $x$- and $y$-axes as well as along the diagonal ($\theta = -45°$).
As we can see below, the cavity causes characteristic stress distributions. We now want to take a closer look at those along the $x$-axis. For larger distances from the hole, the stresses are approximately distributed as if the plate was continuous. As the hole is approached, the tangential stress is increasing until it reaches its maximum value directly at the contour. Interestingly, that value is three times higher than the applied traction. It can be concluded, that the hole leads to a threefold stress concentration in the plate.
%% Cell type:code id:d7b9c91c tags:
``` python
#HIDDEN
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
#Some plot settings
plt.style.use('seaborn-deep')
plt.rcParams['lines.linewidth']= 2.0
plt.rcParams['lines.color']= 'black'
plt.rcParams['legend.frameon']=True
plt.rcParams['font.family'] = 'serif'
plt.rcParams['legend.fontsize']=14
plt.rcParams['font.size'] = 14
plt.rcParams['axes.spines.right'] = False
plt.rcParams['axes.spines.top'] = False
plt.rcParams['axes.spines.left'] = True
plt.rcParams['axes.spines.bottom'] = True
plt.rcParams['axes.axisbelow'] = True
plt.rcParams['figure.figsize'] = (8, 6)
```
%% Cell type:code id:c30056fc tags:
``` python
#HIDDEN
def sig_rr(sig,r,theta,a):
return 0.5*sig*((1-a**2/r**2)+(1+3*np.power(a,4)/np.power(r,4)-4*a**2/r**2)*np.cos(2*np.pi*theta/180)) * np.heaviside(r-a,1)
def sig_tt(sig,r,theta,a):
return 0.5*sig*((1+a**2/r**2)-(1+3*np.power(a,4)/np.power(r,4))*np.cos(2*np.pi*theta/180)) * np.heaviside(r-a,1)
def sig_rt(sig,r,theta,a):
return -0.5*sig*((1-3*np.power(a,4)/np.power(r,4)+2*a**2/r**2)*np.sin(2*np.pi*theta/180)) * np.heaviside(r-a,1)
```
%% Cell type:markdown id:19db8675 tags:
### Stress distribution along the $x$-axis ($\theta= -90°$, orthogonal to the load)
%% Cell type:code id:6166f65f tags:
``` python
r = np.linspace(4,30,1000)
plt.plot(r, sig_rr(10,r,-90,4), label = "$\sigma_{rr}$")
plt.plot(r, sig_tt(10,r,-90,4), label = "$\sigma_{\\theta \\theta}$")
plt.plot(r, sig_rt(10,r,-90,4), label = "$\sigma_{r \\theta}$")
plt.plot([4,4], [0,30], color="0.6", linestyle = "--", label = "$Hole\ Radius$")
plt.legend(loc="upper right")
plt.tight_layout()
plt.grid(True)
plt.xlabel('$r$ / mm')
plt.xlim(0,30)
plt.ylabel('$\\sigma$ / MPa')
plt.ylim(0,30);
```
%%%% Output: display_data
%% Cell type:markdown id:1aa5734f tags:
### Stress distribution along the diagonal ($\theta= -45°$)
%% Cell type:code id:9488928e tags:
``` python
r = np.linspace(4,30,1000)
plt.plot(r, sig_rr(10,r,-45,4), label = "$\sigma_{rr}$")
plt.plot(r, sig_tt(10,r,-45,4), label = "$\sigma_{\\theta \\theta}$")
plt.plot(r, sig_rt(10,r,-45,4), label = "$\sigma_{r \\theta}$")
plt.plot([4,4], [0,10], color="0.6", linestyle = "--", label = "$Hole\ Radius$")
plt.legend(loc="upper right")
plt.tight_layout()
plt.grid(True)
plt.xlabel('$r$ / mm')
plt.xlim(0,30)
plt.ylabel('$\\sigma$ / MPa')
plt.ylim(0,10);
```
%%%% Output: display_data
%% Cell type:markdown id:b95e7cd8 tags:
### Stress distribution along the $y$-axis ($\theta=0°$, parallel to the applied tension)
%% Cell type:code id:59145ce9 tags:
``` python
r = np.linspace(4,30,1000)
plt.plot(r, sig_rr(10,r,0,4), label = "$\sigma_{rr}$")
plt.plot(r, sig_tt(10,r,0,4), label = "$\sigma_{\\theta \\theta}$")
plt.plot(r, sig_rt(10,r,0,4), label = "$\sigma_{r \\theta}$")
plt.plot([4,4], [-10,10], color="0.6", linestyle = "--", label = "$Hole\ Radius$")
plt.legend(loc="lower right")
plt.tight_layout()
plt.grid(True)
ax = plt.gca()
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.spines['bottom'].set_position(('data',0))
ax.yaxis.set_ticks_position('left')
ax.spines['left'].set_position(('data',0))
plt.xlabel('$r$ / mm')
plt.xlim(0,30)
plt.ylabel('$\\sigma$ / MPa')
plt.ylim(-10,10);
```
%%%% Output: display_data
%% Cell type:code id:90a1c879-bda2-43c8-abbb-6f05f6afbfc5 tags:
``` python
cart_to_cyl = lambda x,y: [np.sqrt(x**2+y**2), np.rad2deg(np.arctan(x/y))]
```
%% Cell type:code id:328f47df-1d12-4bad-acc4-7437528cb734 tags:
``` python
X, Y = np.meshgrid(np.linspace(.1,30,1000),np.linspace(0.1,30,1000))
```
%% Cell type:code id:48ecaf5c-e7de-474e-9626-4d32e7c804c1 tags:
``` python
#cylindical coordinates from Cartesian grid
r, theta = cart_to_cyl(X,Y)
eps = 1e-2
fig, ax = plt.subplots(ncols=3,figsize=(18,6))
l1=ax[0].contourf(X,Y,np.ma.masked_array(sig_rr(10,r,-theta,4),mask=r<4))
l2=ax[1].contourf(X,Y,np.ma.masked_array(sig_tt(10,r,-theta,4),mask=r<4))
l3=ax[2].contourf(X,Y,np.ma.masked_array(sig_rt(10,r,-theta,4),mask=r<4))
fig.colorbar(l1,ax=ax[0])
fig.colorbar(l2,ax=ax[1])
fig.colorbar(l3,ax=ax[2])
for i in range(3):
ax[i].set_aspect('equal')
ax[i].set_xlabel('$x$ / mm')
ax[i].set_ylabel('$y$ / mm')
ax[0].set_title('$\\sigma_{rr}$ / MPa')
ax[1].set_title('$\\sigma_{\\theta\\theta}$ / MPa')
ax[2].set_title('$\\sigma_{r\\theta}$ / MPa')
fig.tight_layout();
```
%%%% Output: display_data
%% Cell type:markdown id:03953cb7-452d-4e6a-9371-4a553c83473f tags:
**ToDo**
* add missing information (geometry, parameters)
* check consistency of units (e.g. mm, N, s, MPa=N/mm²)
* maybe we can use contour plots to show point-wise difference to numerical result
* confirm whether $x$ is at $+$ or $-90°$
%% Cell type:code id:3de96148-11b0-4777-8f56-5193bca3555b tags:
``` python
from ogs6py import ogs
import vtuIO
```
%% Cell type:code id:390fcbf7-8c0e-4396-888b-4589e815300f tags:
``` python
model=ogs.OGS(INPUT_FILE="../disc_with_hole.prj", PROJECT_FILE="../disc_with_hole.prj")
```
%% Cell type:code id:753661c8-105f-4b94-a65f-a6ba6841833a tags:
``` python
model.run_model(logfile="out.txt")
```
%%%% Output: stream
OGS finished with project file ../disc_with_hole.prj.
Execution took 0.878124475479126 s
Execution took 0.8345866203308105 s
%% Cell type:code id:831ac74d-e16f-49d6-82e1-1db5e6284237 tags:
``` python
```
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment