This page was generated from docs\source\user_guide/user_guide.ipynb.

User Guide#

This part of the documentation teaches you how to use the tools provided by the library through simple examples, with some background information.

The typical workflow in any simulation in the field of computational solid mechanics consists of the following steps:

  1. Setting up the solid domains and defining material properties.

    The tools in the library are limited to simply supported beams and plates.

  2. Defining boundary conditions (loads and supports).

    You only have to set up loads here, since the essential boundary conditions (displacement supports) are automatically satisfied by the approximation model.

  3. Performing some calculations.

    Currently, the library supports linear elastic calculations.

  4. Visualization of the results.

    Check out the examples gallery.

Setting up the domains#

Euler-Bernoulli and Timoshenko beams#

The library provides solutions for both Euler-Bernoully and Timoshenko beams, and setting up any of the two is quite similar. If you provide shear stiffness, the beam is a Timoshenko one, otherwise it is Euler-Bernoulli. The following code snippet illustrates the steps using the NavierBeam class.

[1]:

from sigmaepsilon.solid.fourier import NavierBeam

# geometry
beam_length = 1000.0
width, height = 20.0, 80.0  # rectangular cross-section

# material properties
young_modulus, poisson_ratio = 210000.0, 0.25  # material

# stiffness properties
inertia = width * height**3 / 12
bending_stiffness = young_modulus * inertia

# solution parameters
number_of_modes = 100

bernoulli_beam = NavierBeam(beam_length, number_of_modes, EI=bending_stiffness)

If your beam is of the Timoshenko type, you also have to provide the shear stiffness of the cross section.

[2]:

# calculate shear stiffness
area = width * height
shear_modulus = young_modulus / (2 * (1 + poisson_ratio))
shear_correction_factor = 5 / 6
shear_stiffness = shear_modulus * area * shear_correction_factor

timoshenko_beam = NavierBeam(
    beam_length, number_of_modes, EI=bending_stiffness, GA=shear_stiffness
)

In both cases, you also had to specify the number of harmonic terms involved in the solution at creating the instances for your beams.

Kirchoff-Love and Uflyand-Mindlin plates#

You can define Kirchhoff-Love or Uflyand-Mindlin plates using the NavierPlate class. Things are a bit more complicated, as there are more stiffness parameters to calculate and both the bending and shear stiffness terms have to be provided with matrices, depending on the type of plate you want. If you know how to, you can calculate these stiffness matrices on your own, or you can use tools from other solutions from the SigmaEpsilon ecosystem. We definitively suggest you to utilize all that SigmaEpsilon has to offer.

First we create a Kirchhoff-Love plate. For this, you have to specify the material stiffness matrix for the plate, which, for a linearly isotropic material takes the form of

\[\begin{split}D = \frac{E h^3}{12 (1 - \nu^2)} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1 - \nu}{2} \end{bmatrix}\end{split}\]

where \(E\) is the Young’s modulus, \(\nu\) is the Poisson’s ratio and \(h\) is the thickness. Of course, this easy to do manually.

[3]:

import numpy as np
from sigmaepsilon.solid.fourier import NavierPlate

# geometry
length_X, length_Y = (600.0, 800.0)
thickness = 25.0

# material properties
young_modulus = 2890.0
poisson_ratio = 0.2

# solution parameters
number_of_modes_X = 20
number_of_modes_Y = 20

# stiffness
D_factor = young_modulus * thickness**3 / (12 * (1 - poisson_ratio**2))
bending_stiffness = D_factor * np.array(
    [[1, poisson_ratio, 0], [poisson_ratio, 1, 0], [0, 0, (1 - poisson_ratio) / 2]]
)

kirchhoff_plate = NavierPlate(
    (length_X, length_Y),
    (number_of_modes_X, number_of_modes_Y),
    D=bending_stiffness,
)

If you want more power, you can use the tools from other libraries in the SigmaEpsilon ecosystem.

[4]:

from sigmaepsilon.math.linalg import ReferenceFrame

from sigmaepsilon.solid.material import (
    KirchhoffPlateSection,
    ElasticityTensor,
    LinearElasticMaterial,
    HuberMisesHenckyFailureCriterion_SP,
)
from sigmaepsilon.solid.material.utils import elastic_stiffness_matrix

from sigmaepsilon.solid.fourier import NavierPlate

# GEOMETRY
length_X, length_Y = (600.0, 800.0)
thickness = 25.0

# MATERIAL PROPERTIES
young_modulus = 2890.0
poisson_ratio = 0.2
yield_strength = 2.0

# SOLUTION PARAMETERS
number_of_modes_X = 20
number_of_modes_Y = 20

# SETTING UP HOOKE'S LAW
hooke = elastic_stiffness_matrix(E=young_modulus, NU=poisson_ratio)
frame = ReferenceFrame(dim=3)
stiffness = ElasticityTensor(hooke, frame=frame, tensorial=False)
failure_model = HuberMisesHenckyFailureCriterion_SP(yield_strength=yield_strength)
material = LinearElasticMaterial(stiffness=stiffness, failure_model=failure_model)

# CALCULATING SECTION STIFFNESS
kirchhoff_section = KirchhoffPlateSection(
    layers=[
        KirchhoffPlateSection.Layer(material=material, thickness=thickness),
    ]
)
bending_stiffness = kirchhoff_section.elastic_stiffness_matrix()

kirchhoff_plate = NavierPlate(
    (length_X, length_Y),
    (number_of_modes_X, number_of_modes_Y),
    D=bending_stiffness,
)

Similarly to beams, if you want to set up an Uflyand-Mindlin plate, the only difference is that you also have to provide the shear stiffness terms at instantiation. These are described by the matrix

\[\begin{split}G = \kappa \frac{E h}{2 (1 + \nu)} \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}\end{split}\]

where \(\kappa\) is the shear correction factor.

[5]:

D_factor = young_modulus * thickness**3 / (12 * (1 - poisson_ratio**2))
bending_stiffness = D_factor * np.array(
    [[1, poisson_ratio, 0], [poisson_ratio, 1, 0], [0, 0, (1 - poisson_ratio) / 2]]
)
G_factor = (5 / 6) * young_modulus * thickness / (2 * (1 + poisson_ratio))
shear_stiffness = G_factor * np.array([[1, 0], [0, 1]])

mindlin_plate = NavierPlate(
    (length_X, length_Y),
    (number_of_modes_X, number_of_modes_Y),
    D=bending_stiffness,
    S=shear_stiffness,
)

Again, you can insted use the tools from the SigmaEpsilon ecosystem.

[6]:

from sigmaepsilon.solid.material import MindlinPlateSection
from numpy import ascontiguousarray as ascont

mindlin_section = MindlinPlateSection(
    layers=[
        MindlinPlateSection.Layer(material=material, thickness=thickness),
    ]
)
ABDS_matrix = mindlin_section.elastic_stiffness_matrix()
bending_stiffness, shear_stiffness = (
    ascont(ABDS_matrix[:3, :3]),
    ascont(ABDS_matrix[3:, 3:]),
)

mindlin_plate = NavierPlate(
    (length_X, length_Y),
    (number_of_modes_X, number_of_modes_Y),
    D=bending_stiffness,
    S=shear_stiffness,
)

Specifying boundary conditions#

As mentioned before, you only have to care about defining the natural boundary conditions (aka. loads) of the problem. Regardless of the type of domain you have, defining the loads is quite similar, but there are some subtle differences. In both cases you need to organize your loads into groups.

Loads for beams#

For beams, you can define concentrated, or linear loads using the classes PointLoad and LineLoad classes. In both cases, you have to provide values for the vertical load $F_y$ and the bending moment $M_z$.

[7]:

from sigmaepsilon.solid.fourier import LoadGroup, PointLoad, LineLoad

beam_loads = LoadGroup(
    concentrated=LoadGroup(
        LC1=PointLoad(beam_length / 2, [1.0, 0.0]),
        LC2=PointLoad(beam_length / 2, [0.0, 1.0]),
    ),
    distributed=LoadGroup(
        LC3=LineLoad([0, beam_length], [1.0, 0.0]),
        LC4=LineLoad([beam_length / 2, beam_length], [0.0, 1.0]),
    ),
)

With the LineLoad class, you can also define arbitrary functions by passing strings for the load values. For instance, to define a sinusoidal vertical load on the right side of the beam:

[8]:

load_case = LineLoad([beam_length / 2, beam_length], ["sin(x)", 0.0])

In these cases, the load coefficients are calculated using a Monte-Carlo simulation.

Loads for plates#

For plates, you can define concentrated, or distributed loads over a rectangle using the classes PointLoad and RectangleLoad classes. In these cases, you must provide values for the vertical load $F_z$ and the two moments $M_y$ and $M_z$.

[9]:

from sigmaepsilon.solid.fourier import RectangleLoad

plate_loads = LoadGroup(
    LC1=PointLoad([length_X / 2, length_Y / 2], [-1.0, 0.0, 0.0]),
    LC2=RectangleLoad([[0, 0], [length_X, length_Y]], [-0.1, 0, 0]),
)

Linear Static Analysis#

Besides the loads, you also have to provide the points of evaluation where you want the results to be calculated. To perform a linear static analysis, call the linear_static_analysis instance method.

[10]:

import numpy as np

# the locations of the points where the solution is calculated
beam_points = np.linspace(0, beam_length, 500)

bernoulli_solution = bernoulli_beam.linear_static_analysis(beam_points, beam_loads)
timoshenko_solution = timoshenko_beam.linear_static_analysis(beam_points, beam_loads)

For plates, we can calculate the locations of the points of evaluation using numpy.meshgrid.

[11]:

x = np.linspace(0, length_X, 30)
y = np.linspace(0, length_Y, 40)
xv, yv = np.meshgrid(x, y)
plate_points = np.stack((xv.flatten(), yv.flatten()), axis=1)

kirchhoff_results = kirchhoff_plate.linear_static_analysis(plate_points, plate_loads)

Understanding the results#

The linear_static_analysis method in both cases returns a dictionary with an identical layout as the loads. So for instance, if a load case was accessible as beam_loads["concentrated", "LC1"], then the corresponding results are accessible as in bernoulli_solution["concentrated", "LC1"]. Each of these entries is an instance of BeamLoadCaseResultLinStat.

[12]:

type(bernoulli_solution)

[12]:

sigmaepsilon.deepdict.deepdict.DeepDict

[13]:

type(bernoulli_solution["concentrated", "LC1"])

[13]:

sigmaepsilon.solid.fourier.result.BeamLoadCaseResultLinStat

The results as a numpy.ndarray#

The out-of-the-box way of representing data is via NumPy arrays. This means a 2d NumPy array for each load case, where the first axis goes along the points of evaluation and the second along the components (strains, internal forces, etc.). This array is accessible through the data and values properties.

[14]:

# 6 results components for 500 points
bernoulli_solution["concentrated", "LC1"].data.shape

[14]:

(500, 6)

It is not very convenient to memorize the indices of the components. If you want a reminder, you can use the components property.

[15]:

bernoulli_solution["concentrated", "LC1"].components

[15]:

['UY', 'ROTZ', 'CZ', 'EXY', 'MZ', 'SY']

For more details, see the documentation of the classes BeamLoadCaseResultLinStat and PlateLoadCaseResultLinStat.

The results as an xarray.DataArray#

You can turn the results into an instance of xarray.DataArray by calling the to_xarray method of the instance (this requires xarray to be installed):

[16]:

xarr = bernoulli_solution["concentrated", "LC1"].to_xarray()

If you are familiar with xarray, you know that it is basically the same as a NumPy array, except that the values are also accessible using labels, besides the usual index-based nature of NumPy. The indexes attribute tells you about the available options.

[17]:

xarr.indexes

[17]:

Indexes:
    index      Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
       ...
       490, 491, 492, 493, 494, 495, 496, 497, 498, 499],
      dtype='int32', name='index', length=500)
    component  Index(['UY', 'ROTZ', 'CZ', 'EXY', 'MZ', 'SY'], dtype='object', name='component')

You can see that the array is two dimensional and the two axes are called ‘index’ and ‘component’. Here ‘index’ refers to the integer index of the point of evaluation, while ‘component’ refers to result component associated with a beam.

Of course, for Euler-Bernoulli beams, the shear strains and shear forces are zero. To access the vertical displacement of the 10th evaluation point, use the loc method of the DataArray like this:

[18]:

xarr.loc[9, "UY"]

[18]:

<xarray.DataArray 'values' ()> Size: 8B
array(6.28775469e-06)
Coordinates:
    index      int32 4B 9
    component  <U4 16B 'UY'

The same query with the sel method:

[19]:

xarr.sel(index=9, component="UY")

[19]:

<xarray.DataArray 'values' ()> Size: 8B
array(6.28775469e-06)
Coordinates:
    index      int32 4B 9
    component  <U4 16B 'UY'

Important

Like pandas, label based indexing in xarray is inclusive of both the start and stop bounds. For instance, to access the same result component for the first 10 points, you have to slice the data with [:9, 'UY'] instead of [:10, 'UY']:

[20]:

xarr.loc[:9, "UY"]

[20]:

<xarray.DataArray 'values' (index: 10)> Size: 80B
array([0.00000000e+00, 6.98938448e-07, 1.39785457e-06, 2.09672599e-06,
       2.79553025e-06, 3.49424481e-06, 4.19284710e-06, 4.89131451e-06,
       5.58962453e-06, 6.28775469e-06])
Coordinates:
  * index      (index) int32 40B 0 1 2 3 4 5 6 7 8 9
    component  <U4 16B 'UY'

You can also use the sel method to to the same thing:

[21]:

xarr.sel(index=slice(9), component="UY")

[21]:

<xarray.DataArray 'values' (index: 10)> Size: 80B
array([0.00000000e+00, 6.98938448e-07, 1.39785457e-06, 2.09672599e-06,
       2.79553025e-06, 3.49424481e-06, 4.19284710e-06, 4.89131451e-06,
       5.58962453e-06, 6.28775469e-06])
Coordinates:
  * index      (index) int32 40B 0 1 2 3 4 5 6 7 8 9
    component  <U4 16B 'UY'

The results as a pandas.DataFrame#

You can also turn the results into a pandas.DataFrame:

[22]:

df = bernoulli_solution["concentrated", "LC1"].to_pandas()
df

[22]:
UY ROTZ CZ EXY MZ SY
index
0 0.000000e+00 3.487721e-07 0.000000e+00 0.0 0.000000e+00 0.496817
1 6.989384e-07 3.487666e-07 -5.558252e-12 0.0 -9.960388e-01 0.497428
2 1.397855e-06 3.487499e-07 -1.112927e-11 0.0 -1.994365e+00 0.499024
3 2.096726e-06 3.487220e-07 -1.672093e-11 0.0 -2.996390e+00 0.500996
4 2.795530e-06 3.486828e-07 -2.233318e-11 0.0 -4.002106e+00 0.502586
... ... ... ... ... ... ...
495 2.795530e-06 -3.486828e-07 -2.233318e-11 0.0 -4.002106e+00 -0.502586
496 2.096726e-06 -3.487220e-07 -1.672093e-11 0.0 -2.996390e+00 -0.500996
497 1.397855e-06 -3.487499e-07 -1.112927e-11 0.0 -1.994365e+00 -0.499024
498 6.989384e-07 -3.487666e-07 -5.558252e-12 0.0 -9.960388e-01 -0.497428
499 1.320568e-20 -3.487721e-07 8.853895e-26 0.0 1.586618e-14 -0.496817

500 rows × 6 columns

Then you can query the object quite similarly to a DataArray.

[23]:

df.loc[9, 'UY']

[23]:

6.2877546906128156e-06

[24]:

df.loc[:9, 'UY']

[24]:

index
0    0.000000e+00
1    6.989384e-07
2    1.397855e-06
3    2.096726e-06
4    2.795530e-06
5    3.494245e-06
6    4.192847e-06
7    4.891315e-06
8    5.589625e-06
9    6.287755e-06
Name: UY, dtype: float64
On this page
Show Source