Skip to content
Snippets Groups Projects
Unverified Commit 8954fffe authored by Lars Bilke's avatar Lars Bilke
Browse files

[nb] Fix image.

parent 61c0b52f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Tests/Data/Mechanics/CooksMembrane/CooksMembraneBbar.ipynb
+ 2
− 2
View file @ 8954fffe
%% Cell type:raw id:73c13b4b-fee8-44b4-8a30-427afac95d32 tags:
 
+++
title = "Cook's membrane example"
date = "2024-06-11"
author = "Wenqing Wang"
image = "cooks_membrane.png"
image = "figures/cooks_membrane.png"
web_subsection = "small-deformations"
weight = 3
+++
 
%% Cell type:markdown id:286c4e5e-eb58-409e-ab9b-f2c42fd388b4 tags:
 
$$
\newcommand{\B}{\text{B}}
\newcommand{\F}{\text{F}}
\newcommand{\I}{\mathbf I}
\newcommand{\intD}[1]{\int_{\Omega_e}#1\mathrm{d}\Omega}
$$
 
# Cook's membrane example for nearly icompressible solid
 
## B bar method
Considering a strain decomposition: $\mathbf\epsilon = \underbrace{\mathbf\epsilon- \frac{1}{3}(\epsilon:\mathbf I)}_{\text{deviatoric}}\I + \underbrace{\frac{1}{3}(\epsilon:\mathbf I)}_{\text{dilatational}} \I$.
The idea of the B bar method is to use another quadrature rule to interpolate the dilatational part, which leads to a modified B matrix [1]:
$$
\bar\B = \underbrace{\B - \B^{\text{dil}}}_{\text{original B elements}}+ \underbrace{{\bar\B}^{\text{dil}}}_{\text{by another quadrature rule} }
$$
There are several methods to form ${\bar\B}^{\text{dil}}$ such as selective integration, generalization of the mean-dilatation formulation. In the current OGS, we use the latter, which reads
$$
{\bar\B}^{\text{dil}} = \frac{\intD{\B^{\text{dil}}(\xi)}}{\intD{}}
$$
 
## Example
To verify the implementation of the B bar method, the so called Cook's membrane is used as a benchmark.
Illustrated in the following figure, this example simulates a tapered
and swept panel of unit thickness. The left edge is clamped and the right edge is applied with a distributed shearing load $F$ = 100 N/mm. The plane strain condition is considered. This numerical model is exactly the same as that is presented in the paper by T. Elguedj et al [1,2].
 
<img src="./figures/cooks_membrane.png" alt="Cook's membrane" width="320" height="320" />
<img src="figures/cooks_membrane.png" alt="Cook's membrane" width="320" height="320" />
 
## Reference
 
[1] T.J.R. Hughes (1980). Generalization of selective integration procedures to anisotropic and nonlinear media. International Journal for Numerical Methods in Engineering, 15(9), 1413-1418.
 
[2] T. Elguedj, Y. Bazilevs, V.M. Calo, T.J.R. Hughes (2008),
$\bar\B$ and $\bar\F$ projection methods for nearly incompressible linear and non-linear elasticity and plasticity using higher-order NURBS elements, Computer Methods in Applied Mechanics and Engineering, 197(33--40), 2732-2762.
 
%% Cell type:code id:f0ae03a5-4cb3-43aa-91f1-5f25e2de61bc tags:
 
``` python
import os
from pathlib import Path
 
from ogs6py.ogs import OGS
 
out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
if not out_dir.exists():
out_dir.mkdir(parents=True)
```
 
%% Cell type:code id:35a460be-3080-4f1f-9285-9e682fe20d38 tags:
 
``` python
import xml.etree.ElementTree as ET
 
import pyvista as pv
```
 
%% Cell type:code id:2bbde8c5-9907-43c2-a26c-61a5ab512a6d tags:
 
``` python
def get_last_vtu_file_name(pvd_file_name):
tree = ET.parse(Path(out_dir) / pvd_file_name)
root = tree.getroot()
# Get the last DataSet tag
last_dataset = root.findall(".//DataSet")[-1]
 
# Get the 'file' attribute of the last DataSet tag
file_attribute = last_dataset.attrib["file"]
return f"{out_dir}/" + file_attribute
 
 
def get_top_uy(pvd_file_name):
top_point = (48.0e-3, 60.0e-3, 0)
file_name = get_last_vtu_file_name(pvd_file_name)
mesh = pv.read(file_name)
p_id = mesh.find_closest_point(top_point)
u = mesh.point_data["displacement"][p_id]
return u[1]
```
 
%% Cell type:code id:a6d91f6e-b0b7-4ed0-a673-4cc3581a19bb tags:
 
``` python
def run_single_test(mesh_name, output_prefix, use_bbar="false"):
model = OGS(INPUT_FILE="CooksMembrane.prj", PROJECT_FILE=f"{out_dir}/modified.prj")
model.replace_text(mesh_name, xpath="./mesh")
model.replace_text(use_bbar, xpath="./processes/process/use_b_bar")
model.replace_text(output_prefix, xpath="./time_loop/output/prefix")
model.replace_text(
"BiCGSTAB", xpath="./linear_solvers/linear_solver/eigen/solver_type"
)
model.replace_text("ILUT", xpath="./linear_solvers/linear_solver/eigen/precon_type")
vtu_file_name = output_prefix + "_ts_1_t_1.000000.vtu"
model.replace_text(vtu_file_name, xpath="./test_definition/vtkdiff[1]/file")
model.replace_text(vtu_file_name, xpath="./test_definition/vtkdiff[2]/file")
model.replace_text(vtu_file_name, xpath="./test_definition/vtkdiff[3]/file")
 
model.write_input()
 
# Run OGS
model.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir} -m .")
 
# Get uy at the top
return get_top_uy(output_prefix + ".pvd")
```
 
%% Cell type:code id:37127e0d-10dd-4c7a-b263-0fe3e42b6b64 tags:
 
``` python
mesh_names = [
"mesh.vtu",
"mesh_n10.vtu",
"mesh_n15.vtu",
"mesh_n20.vtu",
"mesh_n25.vtu",
"mesh_n30.vtu",
]
output_prefices_non_bbar = [
"cooks_membrane_sd_edge_div_4_non_bbar",
"cooks_membrane_sd_refined_mesh_10_non_bbar",
"cooks_membrane_sd_refined_mesh_15_non_bbar",
"cooks_membrane_sd_refined_mesh_20_non_bbar",
"cooks_membrane_sd_refined_mesh_25_non_bbar",
"cooks_membrane_sd_refined_mesh_30_non_bbar",
]
 
uys_at_top_non_bbar = []
for mesh_name, output_prefix in zip(mesh_names, output_prefices_non_bbar):
uy_at_top = run_single_test(mesh_name, output_prefix)
uys_at_top_non_bbar.append(uy_at_top)
 
print(uys_at_top_non_bbar)
```
 
%% Output
 
OGS finished with project file output/modified.prj.
Execution took 0.21584200859069824 s
OGS finished with project file output/modified.prj.
Execution took 0.08944129943847656 s
OGS finished with project file output/modified.prj.
Execution took 0.09354853630065918 s
OGS finished with project file output/modified.prj.
Execution took 0.10343599319458008 s
OGS finished with project file output/modified.prj.
Execution took 0.1191704273223877 s
OGS finished with project file output/modified.prj.
Execution took 0.1564652919769287 s
[0.0021645867841231024, 0.0022603329644579387, 0.0023752958560671676, 0.002519725590136147, 0.00265152941337909, 0.0028682896170252165]
 
%% Cell type:code id:99c34461-1d0a-42f7-a716-86e9bac7c046 tags:
 
``` python
output_prefices = [
"cooks_membrane_sd_edge_div_4",
"cooks_membrane_sd_refined_mesh_10",
"cooks_membrane_sd_refined_mesh_15",
"cooks_membrane_sd_refined_mesh_20",
"cooks_membrane_sd_refined_mesh_25",
"cooks_membrane_sd_refined_mesh_30",
]
 
uys_at_top_bbar = []
for mesh_name, output_prefix in zip(mesh_names, output_prefices):
uy_at_top = run_single_test(mesh_name, output_prefix, "true")
uys_at_top_bbar.append(uy_at_top)
 
print(uys_at_top_bbar)
```
 
%% Output
 
OGS finished with project file output/modified.prj.
Execution took 0.07093119621276855 s
OGS finished with project file output/modified.prj.
Execution took 0.07613062858581543 s
OGS finished with project file output/modified.prj.
Execution took 0.08341574668884277 s
OGS finished with project file output/modified.prj.
Execution took 0.10603880882263184 s
OGS finished with project file output/modified.prj.
Execution took 0.12837624549865723 s
OGS finished with project file output/modified.prj.
Execution took 0.16935420036315918 s
[0.006798855415340229, 0.007728027781081195, 0.00787252293068606, 0.007934707855031697, 0.007963259983774562, 0.007989988696891803]
 
%% Cell type:code id:7c55e390-e329-45ac-b3d0-b4f785059672 tags:
 
``` python
import matplotlib.pyplot as plt
import numpy as np
 
ne = [4, 10, 15, 20, 25, 30]
 
 
def plot_data(ne, u_y_bbar, uy_non_bbar, file_name=""):
# Plotting
plt.rcParams["figure.figsize"] = [5, 5]
 
if len(u_y_bbar) != 0:
plt.plot(
ne, np.array(u_y_bbar) * 1e3, marker="o", linestyle="dashed", label="B bar"
)
if len(uy_non_bbar) != 0:
plt.plot(
ne,
np.array(uy_non_bbar) * 1e3,
marker="x",
linestyle="dashed",
label="non B bar",
)
 
plt.xlabel("Number of elements per side")
plt.ylabel("Top right corner displacement /mm")
plt.legend()
 
plt.tight_layout()
if file_name != "":
plt.savefig(file_name)
plt.show()
```
 
%% Cell type:markdown id:c74f1383-596b-4e2c-93f9-61b41248ca2f tags:
 
## Result
 
### Vertical diplacement at the top point
 
The following figure shows that the convergence of the solutions obtained by using the B bar method follows the one presented in the paper by T. Elguedj et al [1]. However, the results obtained without the B bar method are quit far from the converged solution with the finest mesh.
 
%% Cell type:code id:96ac5068-feae-4c82-90d9-ee4535b08291 tags:
 
``` python
plot_data(ne, uys_at_top_bbar, uys_at_top_non_bbar, "b_bar_linear.png")
```
 
%% Output
 
 
%% Cell type:markdown id:399d3f5a-e5bd-4c10-91f5-16a083b100b9 tags:
 
### Contour plot
 
%% Cell type:code id:2199724b-4bb3-4184-8e77-c40fc5c87832 tags:
 
``` python
import matplotlib.tri as tri
import vtuIO
 
nedges = ["4", "10", "15", "20", "25", "30"]
 
 
def contour_plot(pvd_file_name, title):
file_name = get_last_vtu_file_name(pvd_file_name)
m_plot = vtuIO.VTUIO(file_name, dim=2)
triang = tri.Triangulation(m_plot.points[:, 0], m_plot.points[:, 1])
triang = tri.Triangulation(m_plot.points[:, 0], m_plot.points[:, 1])
s_plot = m_plot.get_point_field("sigma")
s_trace = s_plot[:, 0] + s_plot[:, 1] + s_plot[:, 2]
u_plot = m_plot.get_point_field("displacement")
 
fig, ax = plt.subplots(ncols=2, figsize=(8, 3))
ax[0].set_title(title, loc="left", y=1.12)
plt.subplots_adjust(wspace=0.5)
 
contour_stress = ax[0].tricontourf(triang, s_trace, cmap="viridis")
contour_displacement = ax[1].tricontourf(triang, u_plot[:, 1], cmap="gist_rainbow")
fig.colorbar(contour_stress, ax=ax[0], label="Stress trace / MPa")
fig.colorbar(contour_displacement, ax=ax[1], label="Dispplacement / m")
fig.tight_layout()
plt.savefig(pvd_file_name + ".png")
plt.show()
```
 
%% Cell type:markdown id:a7cabdce-6174-4ab8-96ad-d70cc190929b tags:
 
#### Results obtained without the B bar method:
 
%% Cell type:code id:79199f70-9d17-42f8-83b6-6223247e7080 tags:
 
``` python
for nedge, output_prefix in zip(nedges, output_prefices_non_bbar):
contour_plot(output_prefix + ".pvd", "Number of elements per side: " + nedge)
```
 
%% Output
 
 
 
 
 
 
 
%% Cell type:markdown id:9ffed34d-66ed-4549-92e5-14903c6dd5a5 tags:
 
#### Results obtained with the B bar method:
 
%% Cell type:code id:c65992aa-5533-4eea-9568-eb9d0226439e tags:
 
``` python
for nedge, output_prefix in zip(nedges, output_prefices):
contour_plot(output_prefix + ".pvd", "Number of elements per side: " + nedge)
```
 
%% Output
 
 
 
 
 
 
 
%% Cell type:markdown id:312f9540-d590-4215-ac39-a6d652bf82f7 tags:
 
The contour plots show that even with the coarsest mesh, the B bar method still gives reasonable stress result.
 
%% Cell type:code id:327e543b-b055-4c3f-b604-db4760c2be25 tags:
 
``` python
```
......
web/content/docs/benchmarks/small-deformations/cooksmembranebbar/figures/cooks_membrane.png deleted 120000 → 0
+ 0
− 1
View file @ 61c0b52f
../../../../../../../Tests/Data/Mechanics/CooksMembrane/figures/cooks_membrane.png
\ No newline at end of file
web/content/docs/benchmarks/small-deformations/figures/cooks_membrane.png deleted 120000 → 0
+ 0
− 1
View file @ 61c0b52f
../../../../../../Tests/Data/Mechanics/CooksMembrane/figures/cooks_membrane.png
\ No newline at end of file
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