Boundary Condition Mode

This example will explain the basic usage of the boundary condition mode with saenopy.

This mode solves a boundary condition situation, where the displacement (Dirichlet conditon, fixed nodes) or the force (von Neumann condition, free nodes) on each node is defined.

In other words, the user can provide target displacements for some nodes (e.g. the borders of the mesh with 0 displacement) and target forces for all other nodes (e.g. nodes with 0 force in the bulk of the material).

The input is:

• the material parameters (see Material)

• the mesh (nodes and connectivity) (see Mesh)

• the displacements \(U_i\) for fixed nodes

• the forces \(f_i\) for free nodes

The output is:

• the forces \(f_i\) for fixed nodes

• the displacements \(U_i\) for free nodes

First the basic theoretical background is explained, followed by the code to use it.

Theory

Due to the high non-linearity of the material, the boundary condition problem has to be solved iteratively.

First the stiffness matrix \(K_{ij}(U)\), the nodal forces \(f_i(U)\) and the total energy \(E(U)\) are calculated for all nodes (indices \(i\), \(j\)) for the current displacements \(U\).

\[\begin{split}K_{ij}(U) &= \frac{\partial E}{\partial u_i\partial u_j}\\ f_{i}(U) &= \frac{\partial E}{\partial u_i}\end{split}\]

Then the following linear equation (\(A_{ij}\cdot x_j = b_i\)) is solved with a constant displacement \(U\), for an infinitessimal displacement shift \(\Delta u\).

\[\underbrace{K_{ij}(U)}_{A_{ij}} \cdot \underbrace{\Delta u_j}_{x_j} = \underbrace{f_i(U) - f^\mathrm{ext}_i}_{b_i}\]

Index \(i\) runs over all free nodes, \(j\) over all nodes. Where \(K_{ij}\) is the stiffness matrix, \(\Delta u_i\) an infinitessimal displacement, \(f_i\) the force on the node \(i\) and \(f^\mathrm{ext}_i\) the external force on the node \(i\) (boundary condition).

This equation is solved using the conjugate gradient method to obtain a value for \(\Delta u_i\). A tiny fraction ( stepper) of this is applied to the displacement \(U_i\):

\[U^\prime_i = U_i + \mathrm{stepper} \cdot \Delta u_i.\]

With the new displacement \(U^\prime\), the stiffness matrix \(K_{ij}(U^\prime)\), the nodal forces \(f_i(U^\prime)\), and the total energy \(E(U^\prime)\) are updated.

This procedure is iterated until the total energy \(E\) of the system converges. The convergence criterion is:

\[\frac{\mathrm{std}(\{E^{t-6}, E^{t-5}, ..., E^{t-1}, E^t\})}{\mathrm{mean}(\{E^{t-6}, E^{t-5}, ..., E^{t-1}, E^t\})} \leq \mathrm{rel\_conv\_crit}\]

So when the standard deviation of \(E\) divided by the mean of \(E\) for the last 6 iterations is lower than the threshold rel_conv_crit.

Example

The following code example (the code of all boxed joined in a single file) can be downloaded at boundarycondition.py.

First, import the Solver class. This is the central class of saenopy.

from saenopy import Solver
import saenopy
# initialize the object
M = Solver()

Set the material model (see Material):

from saenopy.materials import SemiAffineFiberMaterial

# provide a material model
material = SemiAffineFiberMaterial(1645, 0.0008, 1.0075, 0.033)
M.setMaterialModel(material)

Define the mesh (see Mesh).

import numpy as np

# define the coordinates of the nodes of the mesh
# the array has to have the shape N_n x 3
R = np.array([[0., 0., 0.],  # 0
              [0., 1., 0.],  # 1
              [1., 1., 0.],  # 2
              [1., 0., 0.],  # 3
              [0., 0., 1.],  # 4
              [0., 1., 1.],  # 5
              [1., 1., 1.],  # 6
              [1., 0., 1.]]) # 7

# define the concetivity of the mesh (only tetrahedra are allowed)
# the array has to have the shape N_t x 4
# every entry is an index referencing a node in R (indices start with 0)
T = np.array([[0, 1, 3, 5],
              [1, 2, 3, 5],
              [0, 5, 3, 4],
              [4, 5, 3, 7],
              [5, 2, 3, 6],
              [3, 5, 6, 7]])

# provide the node data
M.setNodes(R)
# the tetrahedron data
M.setTetrahedra(T)

And define the boundary conditions (see Boundary Conditions).

# the displacement boundary conditions of the nodes
# if a displacement boundary condition is given, the node will be fixed
U = np.array([[  0.  ,   0.  ,   0.  ],  # 0
              [  0.  ,   0.  ,   0.  ],  # 1
              [np.nan, np.nan, np.nan],  # 2
              [np.nan, np.nan, np.nan],  # 3
              [  0.  ,   0.  ,   0.  ],  # 4
              [  0.  ,   0.  ,   0.  ],  # 5
              [np.nan, np.nan, np.nan],  # 6
              [np.nan, np.nan, np.nan]]) # 7

# the force boundary conditions of the nodes
# if a target force boundary condition is given, the node will be free
# this is the force that the material applies after solving onto the nodes
# therefore for a pull to the right (positive x-direction) we have to provide
# a target force to the left (negative x-direction)
F_ext = np.array([[np.nan, np.nan, np.nan],  # 0
                  [np.nan, np.nan, np.nan],  # 1
                  [-2.5  ,  0.   ,  0.   ],  # 2
                  [-2.5  ,  0.   ,  0.   ],  # 3
                  [np.nan, np.nan, np.nan],  # 4
                  [np.nan, np.nan, np.nan],  # 5
                  [-2.5  ,  0.   ,  0.   ],  # 6
                  [-2.5  ,  0.   ,  0.   ]]) # 7

# and the boundary condition
M.setBoundaryCondition(U, F_ext)

After defining all input parameters, we can start with the solving process (for parameters see solve_boundarycondition()).

# relax the mesh and move the "varible" nodes
M.solve_boundarycondition();

We can now view the result with the integrated mesh viewer. The nodes (yellow points) are connected by tetrahedra (yellow lines). The displacements (green vectors) cause the reaction forces (red vectors).

# visualize the meshes
M.viewMesh(50, 0.1)