Plate with a Hole#
Plane stress structural linear analysis.
create and mesh a plate with a hole
define a solid body with a linear-elastic plane stress material
create an external pressure load
export and plot stress results
A plate with length \(2L\), height \(2h\) and a hole with radius \(r\) is subjected to a uniaxial tension \(p=-100\) MPa. What is being looked for is the von Mises stress distribution and the concentration of normal stress \(\sigma_{11}\) over the hole.
Let’s create a meshed plate with a hole out of quad cells with the help of pygmsh
(install with pip install pygmsh
). Only a quarter model of the plate is considered. A boolean difference between a rectangle and a disk results in a plate with a hole. This plate is meshed with triangles. A so-called re-combination of the surface mesh gives a mesh with quad cells.
import felupe as fem
import pygmsh
h = 1.0
L = 2.0
r = 0.3
with pygmsh.occ.Geometry() as geom:
geom.characteristic_length_min = 0.02
geom.characteristic_length_max = 0.02
rectangle = geom.add_rectangle([0, 0, 0], L, h)
disk = geom.add_disk([0, 0, 0], r)
plate = geom.boolean_difference(rectangle, disk)
geom.set_recombined_surfaces(plate)
mesh = geom.generate_mesh()
The points and cells of the above mesh are used to initiate a FElupe mesh.
mesh = fem.Mesh(
points=mesh.points[:, :2], cells=mesh.cells[1].data, cell_type=mesh.cells[1].type
)
A numeric quad-region created on the mesh in combination with a vector-valued displacement field represents the plate. The Boundary conditions for the symmetry planes are generated on the displacement field.
region = fem.RegionQuad(mesh)
displacement = fem.Field(region, dim=2)
field = fem.FieldContainer([displacement])
boundaries = fem.dof.symmetry(displacement)
The material behaviour is defined through a built-in isotropic linear-elastic material formulation for plane stress problems. A solid body applies the linear-elastic material formulation on the displacement field.
umat = fem.LinearElasticPlaneStress(E=210000, nu=0.3)
solid = fem.SolidBody(umat, field)
The external uniaxial tension is applied by a pressure load on the right end at \(x=L\). Therefore, a boundary region in combination with a field has to be created at \(x=L\).
region_boundary = fem.RegionQuadBoundary(mesh, mask=mesh.points[:, 0] == L)
field_boundary = fem.FieldContainer([fem.Field(region_boundary, dim=2)])
load = fem.SolidBodyPressure(field_boundary, pressure=-100)
The equivalent stress von Mises, projected to mesh points, will be added to the result file. For the two-dimensional case it is calculated by
import numpy as np
def von_mises(field, **kwargs):
"Von Mises Stress projected to mesh points."
stress = solid.evaluate.gradient(field)[0]
vonmises = np.sqrt(
stress[0, 0] ** 2
+ stress[1, 1] ** 2
+ 3 * stress[0, 1] ** 2
+ stress[0, 0] * stress[1, 1]
)
return fem.project(vonmises, region)
The simulation model is now ready to be solved. The results are saved within a XDMF-file, where additional point-data is passed to the point_data
argument.
step = fem.Step(items=[solid, load], boundaries=boundaries)
job = fem.Job(steps=[step])
job.evaluate(filename="result.xdmf", point_data={"von Mises Stress": von_mises})
The stress distribution may be viewed either by opening the XDMF time-series file in ParaView or by creating a view on the field-container directly inside the notebook.
view = fem.View(field, point_data={"von Mises Stress": von_mises(field)})
plotter = view.plot("von Mises Stress", show_edges=False)
plotter.show(jupyter_backend="panel")
The normal stress distribution over the hole at \(x=0\) is plotted with matplotlib.
import matplotlib.pyplot as plt
left = mesh.points[:, 0] == 0
plt.plot(
fem.tools.project(solid.results.stress[0], region)[:, 0][left],
mesh.points[:, 1][left] / h,
"o-",
)
plt.xlabel(r"$\sigma_{11}(x=0, y)$ in MPa $\longrightarrow$")
plt.ylabel(r"$y/h$ $\longrightarrow$");