 |
NuTo
Numerics Tool
|
\(\newcommand\vec[1]{\boldsymbol{#1}}\)
From the principle of virtual displacements, it follows for a body in equilibrium:
\[ \delta W_{int}^{(t+1)} = \delta W_{ext}^{(t+1)}, \]
where the variation of the internal energy \(\delta W_{int}^{(t+1)}\) and the variation of the external energy \(\delta W_{ext}^{(t+1)}\) is given by
\begin{align*} \delta W_{int}^{(t+1)} &= \int_V\vec{\sigma}^T(\vec{d}^{(t+1)})\delta\vec{\varepsilon}^{(t+1)}\;dV\\ &= \int_V\left(\vec{\sigma}(\vec{d}^{(t)})+\Delta \vec\sigma\left(\vec{d}^{(t)}+\Delta \vec{d}^{(t+1)}\right)\right)^T\delta\vec\varepsilon^{(t+1)}\;dV \\ \delta W_{ext}^{(t+1)} &= \left(\vec F_{ext}^{(t)}+\Delta \vec{F}_{ext}^{(t+1)}\right)^T \delta \vec{d}^{(t+1)}. \end{align*}
Linearization of the system of equations gives
\[ \int_V \left( \dfrac{\partial \sigma(\vec{d}^{(t)})}{\partial \varepsilon}\Delta\vec\varepsilon^{(t+1)} \right)^T \delta\vec\varepsilon^{(t+1)}\;dV = \left( \Delta \vec F_{ext}^{(t+1)} \right)^T \delta \vec{d}^{(t+1)} - \int_V \left( \vec\sigma^{(t)} \right)^T \delta\vec\varepsilon^{(t+1)}\;dV + \left( \vec{F}_{ext}^{(t)} \right)^T \delta\vec d^{(t+1)}, \]
where the last two summands on the right hand side vanish, if the previous step is in equilibrium. In the general formulation, the vector \(\vec{d}\) includes all the nodal displacements. However, if kinematic boundary conditions have to be applied, the variations of the displacement state have to be kinematically compatible. In general, the linear constraint equations can be written as
\[ \vec{T}_1^{(t+1)} \vec{d}_1^{(t+1)}+\vec{T}_2^{(t+1)} \vec{d}_2^{(t+1)} = \bar{\vec{b}}^{(t+1)}, \]
where the nodal displacements are separated into a vector of independent, unknown DOFs \(\vec{d}_1\) and a vector of dependent DOFs \(\vec{d}_2\). For a simulation using the displacement control approach without additional coupling equations, the matrix \(\vec{T}_1^{(t+1)}\) vanishes and the matrix \(\vec{T}_2^{(t+1)}\) is equal to the identity matrix with \(\bar{\vec{b}}^{(t+1)}\) being the prescribed displacements. In the case of periodic boundary conditions (e.g. \(d_k+d_j=0\)) or in the case of coupling a "hanging node" via displacement constraints (e.g. \(d_l=0.5d_j+0.5d_k\)), both transformation matrices do not vanish. These matrices depend on the load step indicated by \(t\), since these conditions might change between successive load steps (e.g. \(\bar{\vec{b}}\) is modified in each increment of a displacement controlled analysis). The assignment of a degree of freedom to either the dependent or the independent DOFs is sometimes arbitrary. For example in the case of the periodic boundary conditions, either \(d_k\) or \(d_j\) are independent, whereas the other is then the dependent DOF. Using a Gauss elimination procedure, which is equivalent to an inversion of the matrix \(\vec T_2^{(t+1)}\), the vector of total displacements \(\vec{d}\) can be expressed as a function of the independent DOFs \(\vec{d}_1\) only:
\[ \vec{d}^{(t+1)} = \begin{bmatrix} \vec{d}_1^{(t+1)}\\ \vec{d}_2^{(t+1)} \end{bmatrix} = \begin{bmatrix} \vec{I}\\ -\vec{C}_{con} \end{bmatrix} \vec{d}_1^{(t+1)} + \begin{bmatrix} \vec{0}\\ \vec{b}^{(t+1)} \end{bmatrix}, \]
where it is assumed that the vector \(\vec d^{(t+1)}\) is reordered so that the independent parameters are at the top. Note the difference between the vector \(\bar{\vec{b}}^{(t+1)}\) and \(\vec{b}^{(t+1)}\). The matrix \(\vec{C}_{con}\) relates the dependent DOFs with the independent DOFs. In the further derivation it is assumed that the matrix \(C\) does not depend on the load step, since the increase of the loading is only obtained by modification of the right hand side \(\vec{b}\).
Consequently, the variation of the nodal displacement vector is given by
\[ \delta\vec{d}^{(t+1)} = \begin{bmatrix} \vec{I}\\ -\vec{C}_{con} \end{bmatrix} \delta\vec{d}_1^{(t+1)}. \]
Using the element shape functions, the relation between nodal displacements \(\vec{d}\) and strains \(\vec{\varepsilon}\) can be expressed as
\begin{align*} \Delta\vec{\varepsilon}^{(t+1)} &= \vec{\varepsilon}^{(t+1)}- \vec{\varepsilon}^{(t)}\\ &= \vec{B}\left( \begin{bmatrix} \vec{I}\\ -\vec{C}_{con} \end{bmatrix} \left( \vec{d}_1^{(t+1)}-\vec{d}_1^{(t)} \right) + \begin{bmatrix} \vec{0}\\ \vec{b}^{(t+1)} \end{bmatrix} - \begin{bmatrix} \vec{0}\\ \vec{b}^{(t)} \end{bmatrix} \right)\\ \delta\vec{\varepsilon}^{(t+1)} &= \vec{B}\left( \begin{bmatrix} \vec{I}\\ -\vec{C}_{con} \end{bmatrix} \delta\vec{d}_1^{(t+1)} \right). \end{align*}
Setting
\[ \Delta\vec{d}_1^{(t+1)}=\left(\vec{d}_1^{(t+1)}-\vec{d}_1^{(t)}\right) \]
and substituting of these equations into the linearized system of the virtual work equations gives
\[ \left( \begin{bmatrix} \vec{I}\\ -\vec{C}_{con} \end{bmatrix} \Delta\vec{d}_1^{(t+1)}+ \begin{bmatrix} \vec{0}\\ \vec{b}^{(t+1)} \end{bmatrix} - \begin{bmatrix} \vec{0}\\ \vec{b}^{(t)} \end{bmatrix} \right)^{T} \vec{K}^T \left( \begin{bmatrix} \vec{I}\\ -\vec{C}_{con} \end{bmatrix} \delta\vec{d}_1^{(t+1)} \right) = \left( \Delta \vec F_{ext}^{(t+1)} + \vec{F}_{ext}^{(t)} - \vec F_{int}^{(t)} \right)^T \begin{bmatrix} \vec{I}\\ -\vec{C}_{con} \end{bmatrix} \delta\vec{d}_1^{(t+1)}. \]
The integrated stiffness matrix and the internal and external force vectors are split into submatrices similar to \(\vec d_1\) and \(\vec d_2\):
\begin{align*} \int_V \vec{B}^T\dfrac{\partial \sigma(\vec{d}^{(t)})}{\partial \varepsilon}\vec{B}\;dV &= \vec K = \begin{bmatrix} \vec K_{11} & \vec K_{12}\\ \vec K_{21} & \vec K_{22} \end{bmatrix}\\ \int_V\vec{B}^T\vec\sigma^{(t)}\;dV&=\vec F_{int}^{(t)} = \begin{bmatrix} \vec F_{1,int}^{(t)}\\ \vec F_{2,int}^{(t)} \end{bmatrix}. \end{align*}
Canceling out \(\delta\vec{d}_1^{(t+1)}\), transposing and splitting the external and internal forces into two parts yields:
\[ \begin{bmatrix} \vec{I}& -\vec{C}_{con}^{T} \end{bmatrix} \begin{bmatrix} \vec K_{11} & \vec K_{12}\\ \vec K_{21} & \vec K_{22} \end{bmatrix} \left( \begin{bmatrix} \vec{I} -\vec{C}_{con} \end{bmatrix} \Delta\vec{d}_1^{(t+1)} + \begin{bmatrix} \vec{0}\\ \vec{b}^{(t+1)} \end{bmatrix} - \begin{bmatrix} \vec{0}\\ \vec{b}^{(t)} \end{bmatrix} \right) = \begin{bmatrix} \vec{I}& -\vec{C}^T_{con} \end{bmatrix} \begin{bmatrix} \Delta \vec F_{1,ext}^{(t+1)}+\vec F_{1,ext}^{(t)}-\vec F_{1,int}^{(t)}\\ \Delta \vec F_{2,ext}^{(t+1)}+\vec F_{2,ext}^{(t)}-\vec F_{2,int}^{(t)} \end{bmatrix}. \]
Expanding this equation finally gives:
\begin{align*} \vec{K}^{mod}\Delta\vec{d}_1^{(t+1)} = \left( \vec{K}_{12}-\vec{C}^T_{con}\vec{K}_{22} \right) \left( \vec{b}^{(t)}-\vec{b}^{(t+1)}\right)+\Delta \vec F_{1,ext}^{(t+1)}+\vec F_{1,ext}^{(t)}\\ -\vec F_{1,int}^{(t)} -\vec{C}^T_{con}\left(\Delta \vec F_2^{(t+1)}+\vec F_{2,ext}^{(t)}-\vec F_{2,int}^{(t)} \right) \end{align*}
with the modified stiffness matrix
\[ \vec K^{mod} = \vec{K}_{11}-\vec{C}^T_{con}\vec{K}_{21}-\vec{K}_{12}\vec{C}_{con}+\vec{C}_{con}^T\vec{K}_{22}\vec{C}_{con}. \]
NuTo::Constraint::Constraints has a container of NuTo::Constaint::Equation.NuTo::Constraint::Equation has a time dependent right-hand-side (std::function<double(double time)>) and a container of NuTo::Constraint::Term.NuTo::NodeBase, a component and a coefficient.That way, NuTo::Constraint::Constraints can build the matrix \(\vec T\) and the rhs \(\vec b\) of
\[ \vec T \vec u = \vec b(t) \]
The method NuTo::Assembler::BuildGlobalDof() performs a sparse (pivot?) Gauss elimination and identifies the dependent dofs \(\vec u_2\), performs renumbering, and so on... After that, the following members are valid:
BlockSparseMatrix mConstraintMatrix - corresponds to \(\vec T_2^{-1} \vec T_1\)BlockSparseMatrix mConstraintRhsMapping - corresponds to \(\vec T_2^{-1}\)BlockFullVector mConstraintRhs - corresponds to \(\vec T_2^{-1} \vec b(t)\)mConstraintRhs is updated to a new global time via NuTo::ConstraintUpdateRhs(time) which internally calculates \(\vec b(t)\) and applies \(\vec T_2^{-1}\) (which is equal, but way faster, than performing the Gauss elimination again ( \(O(n^3)\)?) - at every time step.)
1) Create as many instances of NuTo::Constraint::Equation as you need
mechanics/constraints/ConstraintCompanion.hNuTo::Constraint::Equation or a std::vector of those2) Store the equation
NuTo::Constraint::Constraints::Add(...)NuTo::StructureBase::Constraints() 
1.8.11