![]() |
2D FemElement based on a triangle with a constant thickness More...
#include <SurgSim/Physics/Fem2DElementTriangle.h>
Public Member Functions | |
Fem2DElementTriangle () | |
Constructor. More... | |
Fem2DElementTriangle (std::array< size_t, 3 > nodeIds) | |
Constructor. More... | |
Fem2DElementTriangle (std::shared_ptr< FemElementStructs::FemElementParameter > elementData) | |
Constructor for FemElement object factory. More... | |
void | setThickness (double thickness) |
Sets the triangle's thickness. More... | |
double | getThickness () const |
Gets the triangle's thickness. More... | |
void | initialize (const SurgSim::Math::OdeState &state) override |
Initialize the FemElement once everything has been set. More... | |
double | getVolume (const SurgSim::Math::OdeState &state) const override |
Gets the element volume based on the input state (in m-3) More... | |
SurgSim::Math::Vector | computeCartesianCoordinate (const SurgSim::Math::OdeState &state, const SurgSim::Math::Vector &naturalCoordinate) const override |
Computes a given natural coordinate in cartesian coordinates. More... | |
SurgSim::Math::Vector | computeNaturalCoordinate (const SurgSim::Math::OdeState &state, const SurgSim::Math::Vector &cartesianCoordinate) const override |
Computes a natural coordinate given a global coordinate. More... | |
![]() | |
FemElement () | |
Constructor. More... | |
virtual | ~FemElement () |
Virtual destructor. More... | |
size_t | getNumDofPerNode () const |
Gets the number of degree of freedom per node. More... | |
size_t | getNumNodes () const |
Gets the number of nodes connected by this element. More... | |
size_t | getNodeId (size_t elementNodeId) const |
Gets the elementNodeId-th node id. More... | |
const std::vector< size_t > & | getNodeIds () const |
Gets the node ids for this element. More... | |
void | setYoungModulus (double E) |
Sets the Young modulus (in N.m-2) More... | |
double | getYoungModulus () const |
Gets the Young modulus (in N.m-2) More... | |
void | setPoissonRatio (double nu) |
Sets the Poisson ratio (unitless) More... | |
double | getPoissonRatio () const |
Gets the Poisson ratio (unitless) More... | |
void | setMassDensity (double rho) |
Sets the mass density (in Kg.m-3) More... | |
double | getMassDensity () const |
Gets the mass density (in Kg.m-3) More... | |
double | getMass (const SurgSim::Math::OdeState &state) const |
Gets the element mass based on the input state (in Kg) More... | |
virtual void | addForce (SurgSim::Math::Vector *F, double scale=1.0) const |
Adds the element force (computed for a given state) to a complete system force vector F (assembly) More... | |
virtual void | addMass (SurgSim::Math::SparseMatrix *M, double scale=1.0) const |
Adds the element mass matrix M (computed for a given state) to a complete system mass matrix M (assembly) More... | |
virtual void | addDamping (SurgSim::Math::SparseMatrix *D, double scale=1.0) const |
Adds the element damping matrix D (= -df/dv) (comuted for a given state) to a complete system damping matrix D (assembly) More... | |
virtual void | addStiffness (SurgSim::Math::SparseMatrix *K, double scale=1.0) const |
Adds the element stiffness matrix K (= -df/dx) (computed for a given state) to a complete system stiffness matrix K (assembly) More... | |
virtual void | addFMDK (SurgSim::Math::Vector *F, SurgSim::Math::SparseMatrix *M, SurgSim::Math::SparseMatrix *D, SurgSim::Math::SparseMatrix *K) const |
Adds the element force vector, mass, stiffness and damping matrices (computed for a given state) into a complete system data structure F, M, D, K (assembly) More... | |
virtual void | addMatVec (double alphaM, double alphaD, double alphaK, const SurgSim::Math::Vector &x, SurgSim::Math::Vector *F) const |
Adds the element matrix-vector contribution F += (alphaM.M + alphaD.D + alphaK.K).x (computed for a given state) into a complete system data structure F (assembly) More... | |
bool | isValidCoordinate (const SurgSim::Math::Vector &naturalCoordinate) const |
Determines whether a given natural coordinate is valid. More... | |
template<typename DerivedSub , typename T , int Opt, typename Index > | |
void | assembleMatrixBlocks (const DerivedSub &subMatrix, const std::vector< size_t > blockIds, size_t blockSize, Eigen::SparseMatrix< T, Opt, Index > *matrix, bool initialize=true) const |
Helper method to add a sub-matrix made of squared-blocks into a matrix, for the sake of clarity. More... | |
void | updateFMDK (const Math::OdeState &state, int options) |
Update the FemElement based on the given state. More... | |
Protected Member Functions | |
void | initializeMembers () |
Initializes variables needed before Initialize() is called. More... | |
SurgSim::Math::Matrix33d | computeRotation (const SurgSim::Math::OdeState &state) |
Computes the triangle element's rotation given a state. More... | |
virtual void | computeLocalStiffness (const SurgSim::Math::OdeState &state, Eigen::Matrix< double, 18, 18 > *localStiffnessMatrix) |
Computes the triangle's local stiffness matrix. More... | |
void | computeStiffness (const SurgSim::Math::OdeState &state, SurgSim::Math::Matrix *stiffnessMatrix) |
Computes the triangle's stiffness matrix. More... | |
virtual void | computeLocalMass (const SurgSim::Math::OdeState &state, Eigen::Matrix< double, 18, 18 > *localMassMatrix) |
Computes the triangle's local mass matrix. More... | |
void | computeMass (const SurgSim::Math::OdeState &state, SurgSim::Math::Matrix *massMatrix) |
Computes the triangle's mass matrix. More... | |
void | doUpdateFMDK (const Math::OdeState &state, int options) override |
Update the FemElement based on the given state. More... | |
void | computeShapeFunctionsParameters (const SurgSim::Math::OdeState &restState) |
Compute the various shape functions (membrane and plate deformations) parameters. More... | |
std::array< double, 9 > | batozDhxDxi (double xi, double eta) const |
Batoz derivative dHx/dxi. More... | |
std::array< double, 9 > | batozDhxDeta (double xi, double eta) const |
Batoz derivative dHx/deta. More... | |
std::array< double, 9 > | batozDhyDxi (double xi, double eta) const |
Batoz derivative dHy/dxi. More... | |
std::array< double, 9 > | batozDhyDeta (double xi, double eta) const |
Batoz derivative dHy/deta. More... | |
Matrix39Type | batozStrainDisplacement (double xi, double eta) const |
Batoz strain displacement matrix evaluated at a given point. More... | |
![]() | |
void | setNumDofPerNode (size_t numDofPerNode) |
Sets the number of degrees of freedom per node. More... | |
void | initializeFMDK () |
Initialize f, M, D, K variables. More... | |
virtual void | doInitializeFMDK () |
Function to be overridden by the derived classes to initialize the f, M, D, K variables. More... | |
Protected Attributes | |
Eigen::Matrix< double, 18, 1 > | m_x0 |
The element's rest state. More... | |
SurgSim::Math::Matrix33d | m_initialRotation |
Initial rotation matrix for the element. More... | |
Eigen::Matrix< double, 18, 18 > | m_MLocal |
Stiffness matrix (in local coordinate frame) More... | |
Eigen::Matrix< double, 18, 18 > | m_KLocal |
Stiffness matrix (in local coordinate frame) More... | |
double | m_restArea |
The triangle rest area. More... | |
double | m_thickness |
Thickness of the element. More... | |
SurgSim::Math::Matrix33d | m_membraneShapeFunctionsParameters |
Membrane (in-plane) deformation. More... | |
SurgSim::Math::Vector3d | m_xij |
Batoz variable \(x_{ij} = x_i - x_j\). More... | |
SurgSim::Math::Vector3d | m_yij |
Batoz variable \(y_{ij} = y_i - y_j\). More... | |
SurgSim::Math::Vector3d | m_lij_sqr |
Batoz variable \(l_{ij}^2 = x_{ij}^2 + y_{ij}^2\). More... | |
SurgSim::Math::Vector3d | m_ak |
Batoz variable \(a_{k} = -x_{ij}/l_i^2\). More... | |
SurgSim::Math::Vector3d | m_bk |
Batoz variable \(b_{k} = 3/4 x_{ij} y_{ij}/l_{ij}2\). More... | |
SurgSim::Math::Vector3d | m_ck |
Batoz variable \(c_{k} = (1/4 x_{ij}^2 - 1/2 y_{ij}^2)/l_{ij}^2\). More... | |
SurgSim::Math::Vector3d | m_dk |
Batoz variable \(d_{k} = -y_{ij}/l_{ij}^2\). More... | |
SurgSim::Math::Vector3d | m_ek |
Batoz variable \(e_{k} = (1/4 y_{ij}^2 - 1/2 x_{ij}^2)/l_{ij}^2\). More... | |
SurgSim::Math::Vector3d | m_Pk |
Batoz variable \(P_{k} = -6x_{ij}/l_{ij}^2 = 6.\textrm{m_a}_k\). More... | |
SurgSim::Math::Vector3d | m_qk |
Batoz variable \(q_{k} = 3x_{ij}y_{ij}/l_{ij}^2 = 4.\textrm{m_b}_k\). More... | |
SurgSim::Math::Vector3d | m_tk |
Batoz variable \(t_{k} = -6y_{ij}/l_{ij}^2 = 6.\textrm{m_d}_k\). More... | |
SurgSim::Math::Vector3d | m_rk |
Batoz variable \(r_{k} = 3y_{ij}^2/l_{ij}^2\). More... | |
SurgSim::Math::Matrix | m_integral_dT_d |
Plate mass matrix: integral terms related to the dof \((z)\). More... | |
SurgSim::Math::Matrix | m_integralHyiHyj |
Plate mass matrix: integral terms related to the dof \((\theta_x)\). More... | |
SurgSim::Math::Matrix | m_integralHxiHxj |
Plate mass matrix: integral terms related to the dof \((\theta_y)\). More... | |
![]() | |
size_t | m_numDofPerNode |
Number of degree of freedom per node for this element. More... | |
std::vector< size_t > | m_nodeIds |
Node ids connected by this element. More... | |
double | m_rho |
Mass density (in Kg.m-3) More... | |
double | m_E |
Young modulus (in N.m-2) More... | |
double | m_nu |
Poisson ratio (unitless) More... | |
SurgSim::Math::Vector | m_f |
The force vector. More... | |
SurgSim::Math::Matrix | m_M |
The mass matrix. More... | |
SurgSim::Math::Matrix | m_D |
The damping matrix. More... | |
bool | m_useDamping |
Flag to specify of the damping is used. More... | |
SurgSim::Math::Matrix | m_K |
The stiffness matrix. More... | |
Private Types | |
typedef Eigen::Matrix< double, 3, 3 > | Matrix33Type |
typedef Eigen::Matrix< double, 3, 6 > | Matrix36Type |
typedef Eigen::Matrix< double, 6, 6 > | Matrix66Type |
typedef Eigen::Matrix< double, 3, 9 > | Matrix39Type |
typedef Eigen::Matrix< double, 9, 9 > | Matrix99Type |
Private Member Functions | |
void | computeLocalMembraneMass (const SurgSim::Math::OdeState &state, Eigen::Matrix< double, 18, 18 > *localMassMatrix) |
Computes the triangle's local membrane part mass matrix. More... | |
void | computeLocalPlateMass (const SurgSim::Math::OdeState &state, Eigen::Matrix< double, 18, 18 > *localMassMatrix) |
Computes the triangle's local plate part mass matrix. More... | |
void | computeIntegral_dTd () |
Computes the integral terms d^T.d over the parametrized triangle area. More... | |
void | computeIntegral_HxHxT () |
Computes the integral terms Hy.Hy^T over the parametrized triangle area. More... | |
void | computeIntegral_HyHyT () |
Computes the integral terms Hx.Hx^T over the parametrized triangle area. More... | |
Additional Inherited Members | |
![]() | |
typedef SurgSim::Framework::ObjectFactory1< FemElement, std::shared_ptr< FemElementStructs::FemElementParameter > > | FactoryType |
![]() | |
static FactoryType & | getFactory () |
2D FemElement based on a triangle with a constant thickness
The triangle is modeled as a shell (6DOF) which is decomposed into a membrane (in-plane 2DOF \((x, y)\)) and a plate (bending/twisting 3DOF \((z, \theta_x, \theta_y)\)). The thin-plate assumption does not consider the drilling dof \(\theta_z\). The system includes the drilling DOF for completeness but does not assign any mass or stiffness to it.
The membrane (in-plane) equations (mass and stiffness) are following "Theory of Matrix Structural Analysis" from J.S. Przemieniecki.
The thin-plate (bending) equations (mass and stiffness) are following "A Study Of Three-Node Triangular Plate Bending Elements", Jean-Louis Batoz Numerical Methods in Engineering, vol 15, 1771-1812 (1980).
The plate mass matrix is not detailed in the above paper, but the analytical equations have been derived from it. Moreover, to account for contribution of the displacement along z to the plate mass matrix, we used a cubic expression of this displacement given in: "Shell elements: modelizations DKT, DST, DKTG and Q4g", Code_Aster, 2013, Thomas De Soza.
|
private |
|
private |
|
private |
|
private |
|
private |
SurgSim::Physics::Fem2DElementTriangle::Fem2DElementTriangle | ( | ) |
Constructor.
|
explicit |
Constructor.
nodeIds | An array of 3 node ids defining this triangle element with respect to a DeformableRepresentationState which is passed to the initialize method. |
|
explicit |
Constructor for FemElement object factory.
elementData | A FemElement struct defining this triangle element with respect to a DeformableRepresentationState which is passed to the initialize method. |
SurgSim::Framework::AssertionFailure | if nodeIds has a size different than 3 |
|
protected |
Batoz derivative dHx/deta.
xi,eta | The parametric coordinate (in [0 1] and xi+eta<1.0) |
|
protected |
Batoz derivative dHx/dxi.
xi,eta | The parametric coordinate (in [0 1] and xi+eta<1.0) |
|
protected |
Batoz derivative dHy/deta.
xi,eta | The parametric coordinate (in [0 1] and xi+eta<1.0) |
|
protected |
Batoz derivative dHy/dxi.
xi,eta | The parametric coordinate (in [0 1] and xi+eta<1.0) |
|
protected |
Batoz strain displacement matrix evaluated at a given point.
xi,eta | The parametric coordinate (in [0 1] and xi+eta<1.0) |
|
overridevirtual |
Computes a given natural coordinate in cartesian coordinates.
state | The state at which to transform coordinates |
naturalCoordinate | The coordinates to transform |
Implements SurgSim::Physics::FemElement.
|
private |
Computes the integral terms d^T.d over the parametrized triangle area.
This integral is required in the plate mass matrix computation. The displacement along z is w(x, y) = [d1 d2 d3 d4 d5 d6 d7 d8 d9].U = d.U with di cubic shape functions and U nodal plate displacements.
|
private |
Computes the integral terms Hy.Hy^T over the parametrized triangle area.
This integral is required in the plate mass matrix computation. The displacement along thetay is Thetay(x, y) = -dw/dx = betax = Hx^T.U with Hxi quadratic shape functions and U nodal plate displacements.
|
private |
Computes the integral terms Hx.Hx^T over the parametrized triangle area.
This integral is required in the plate mass matrix computation. The displacement along thetax is Thetax(x, y) = dw/dy = -betay = -Hy^T.U with Hyi quadratic shape functions and U nodal plate displacements.
|
protectedvirtual |
Computes the triangle's local mass matrix.
state | The state to compute the local mass matrix from | |
[out] | localMassMatrix | The local mass matrix to store the result into |
|
private |
Computes the triangle's local membrane part mass matrix.
state | The state to compute the local membrane part mass matrix from | |
[out] | localMassMatrix | The local mass matrix to store the result into |
The mass matrix is derived from the kinetic energy, which depends on the displacements \(\mathbf{u}\).
In the case of triangle membrane, the displacements are a subset of the full shell displacements:
\[ \displaystyle \mathbf{u}(x, y) = \left( \begin{array}{c} u_x(x, y) \\ u_y(x, y) \end{array} \right) = \left( \begin{array}{c} \sum_{i=0}^2 N_i(x, y).u_x^i \\ \sum_{i=0}^2 N_i(x, y).u_y^i \end{array} \right) = \left( \begin{array}{cccccc} N_0(x, y) & 0 & N_1(x, y) & 0 & N_2(x, y) & 0 \\ 0 & N_0(x, y) & 0 & N_1(x, y) & 0 & N_2(x, y) \end{array} \right) \left( \begin{array}{c} u_x^0 \\ u_y^0 \\ u_x^1 \\ u_y^1 \\ u_x^2 \\ u_y^2 \end{array} \right) = a.\mathbf{U} \]
with \(\mathbf{U}\) the triangle membrane nodal displacements vector (subset of the full shell nodal displacements) and \(N_i(x, y) = a_i + b_i.x + c_i.y\) the triangle linear shape functions.
Considering a constant mass density \(\rho\) and a constant thickness \(t\), the kinetic energy of the triangle membrane is \(T = \int_V \frac{1}{2} \rho(x, y, z) \mathbf{\dot{u}}^T.\mathbf{\dot{u}} \, dV = \frac{1}{2}.\rho.t.\int_A \left(\dot{\mathbf{U}}^T.a^T\right).\left(a.\dot{\mathbf{U}}\right) \, dA\).
The triangle membrane mass matrix can be derived from the kinetic energy using the Lagrange equation of motion:
\[ \displaystyle \frac{d \frac{\partial T}{\partial \dot{\mathbf{U}}}}{d t} = \frac{1}{2}.\rho.t.\frac{d \frac{\partial \displaystyle \int_A \left(\dot{\mathbf{U}}^T.a^T\right).\left(a.\dot{\mathbf{U}}\right) \, dA}{\partial \dot{\mathbf{U}}}}{d t} = \frac{1}{2}.\rho.t.\frac{d \displaystyle \int_A 2.a^T.a.\dot{\mathbf{U}} \, dA}{d t} = \left( \rho.t.\int_A a^T.a \, dA \right).\ddot{\mathbf{U}} \]
Therefore, the mass matrix is:
\[ M = \rho.t.\int_A a^T.a \, dA \]
To integrate over the triangle area, we operate the following variables change:
\[ \left\{ \begin{array}{ccc} x(\xi, \eta) & = & x_0 + \xi.(x_1 - x_0) + \eta.(x_2 - x_0) \\ y(\xi, \eta) & = & y_0 + \xi.(y_1 - y_0) + \eta.(y_2 - y_0) \end{array} \right. \]
Therefore, we have \( J = \left( \begin{array}{cc} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \end{array} \right) = \left( \begin{array}{cc} (x_1 - x_0) & (x_2 - x_0) \\ (y_1 - y_0) & (y_2 - y_0) \\ \end{array} \right) \) and \(|J| = 2A\). Therefore the mass matrix becomes
\[ M = \rho.t.\int_0^1 \int_0^{1-\eta} a^T.a |J| \, d\xi \, d\eta = 2.A.\rho.t.\int_0^1 \int_0^{1-\eta} \left( \begin{array}{cccccc} N0.N0 & 0 & N0.N1 & 0 & N0.N2 & 0 \\ 0 & N0.N0 & 0 & N0.N1 & 0 & N0.N2 \\ N1.N0 & 0 & N1.N1 & 0 & N1.N2 & 0 \\ 0 & N1.N0 & 0 & N1.N1 & 0 & N1.N2 \\ N2.N0 & 0 & N2.N1 & 0 & N2.N2 & 0 \\ 0 & N2.N0 & 0 & N2.N1 & 0 & N2.N2 \end{array} \right) \, d\xi \, d\eta \]
Considering the triangle in its local space (i.e. \(x_0=y_0=y_1=0\), \(x_1>0\), \(y_2>0\)), we get:
\[ M = 2.A.\rho.t. \left( \begin{array}{cccccc} \frac{1.0}{12.0} & 0 & \frac{1.0}{24.0} & 0 & \frac{1.0}{24.0} & 0 \\ 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{24.0} & 0 & \frac{1.0}{24.0} \\ \frac{1.0}{24.0} & 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{24.0} & 0 \\ 0 & \frac{1.0}{24.0} & 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{24.0} \\ \frac{1.0}{24.0} & 0 & \frac{1.0}{24.0} & 0 & \frac{1.0}{12.0} & 0 \\ 0 & \frac{1.0}{24.0} & 0 & \frac{1.0}{24.0} & 0 & \frac{1.0}{12.0} \end{array} \right) = \rho.A.t. \left( \begin{array}{cccccc} \frac{1.0}{6.0} & 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{12.0} & 0 \\ 0 & \frac{1.0}{6.0} & 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{12.0} \\ \frac{1.0}{12.0} & 0 & \frac{1.0}{6.0} & 0 & \frac{1.0}{12.0} & 0 \\ 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{6.0} & 0 & \frac{1.0}{12.0} \\ \frac{1.0}{12.0} & 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{6.0} & 0 \\ 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{6.0} \end{array} \right) \]
which needs to be assembled properly in the full shell element mass matrix.
|
private |
Computes the triangle's local plate part mass matrix.
state | The state to compute the local plate part mass matrix from | |
[out] | localMassMatrix | The local mass matrix to store the result into |
The mass matrix is derived from the kinetic energy, which depends on the displacements \(\mathbf{u}\).
In the case of triangle plate, the displacements are a subset of the full shell displacements:
\[ \displaystyle \mathbf{u}(x, y, z) = \left( \begin{array}{c} u(x, y, z) = z.\beta_x(x, y) \\ v(x, y, z) = z.\beta_y(x, y) \\ w(x, y) \end{array} \right) \]
Note that we use the plate formulation DKT (Discrete Kirchhoff Theory) developped by Batoz in 1980. In this formulation, the displacement \(w\) is never formaly expressed, we only know that it is a cubic polynomial function of \((x, y)\). To complete the plate mass matrix calculation, we tried to use the expression of \(w\) from Przemieniecki book "Theory of Matrix Structural Analysis" Chapter 11.10, equation 11.57 (i.e. \(w(x, y) = \mathbf{d}.C^{-1}.\mathbf{U}\)). Unfortunately, it can occur that the matrix \(C\) becomes singular (example of a right isoceles triangle, i.e. \(x_0=y_0=y_1=x_2=0\) and \(x_1=y_2=1\)). To overcome this issue, we prefered the Code_Aster approach developped in "Shell elements: modelizations DKT, DST, DKTG and Q4g", Code_Aster, 2013, Thomas De Soza, where a cubic well defined expression of \(w(x, y)\) is given as: \(w(x, y) = \mathbf{d}.\mathbf{U}\) where \(d_i\) are cubic shape functions.
Therefore, it leads to the following plate displacements:
\[ \mathbf{u}(x, y, z) = \left( \begin{array}{c} u(x, y, z) = z.\beta_x(x, y) = z.\mathbf{H_x}^T(\xi, \eta).\mathbf{U} \\ v(x, y, z) = z.\beta_y(x, y) = z.\mathbf{H_y}^T(\xi, \eta).\mathbf{U} \\ w(x, y) = \mathbf{d}.\mathbf{U} \end{array} \right) = \left( \begin{array}{c} z.\mathbf{H_x}^T(\xi, \eta) \\ z.\mathbf{H_y}^T(\xi, \eta) \\ \mathbf{d} \end{array} \right) . \mathbf{U} = a.\mathbf{U} \]
with \(\mathbf{U}\) the triangle plate nodal displacements vector (subset of the full shell nodal displacements) and \(\mathbf{H_x}\) and \(\mathbf{H_y}\) are the Batoz shape functions.
Considering a constant mass density \(\rho\) and a constant thickness \(t\), the kinetic energy of the triangle plate is \(T = \int_V \frac{1}{2} \rho(x, y, z) \mathbf{\dot{u}}^T.\mathbf{\dot{u}} \, dV = \frac{1}{2}.\rho.\int_{-t/2}^{+t/2} \int_A \left(\dot{\mathbf{U}}^T.a^T\right).\left(a.\dot{\mathbf{U}}\right) \, dA \, dz \).
The triangle plate mass matrix can be derived from the kinetic energy using the Lagrange equation of motion:
\[ \displaystyle \frac{d \frac{\partial T}{\partial \dot{\mathbf{U}}}}{d t} = \frac{1}{2}.\rho.\frac{d \frac{\partial \displaystyle \int_{-t/2}^{+t/2} \int_A \left(\dot{\mathbf{U}}^T.a^T\right).\left(a.\dot{\mathbf{U}}\right) \, dA \, dz}{\partial \dot{\mathbf{U}}}}{d t} = \frac{1}{2}.\rho.\frac{d \displaystyle \int_{-t/2}^{+t/2} \int_A 2.a^T.a.\dot{\mathbf{U}} \, dA \, dz}{d t} \\ = \left( \rho.\int_{-t/2}^{+t/2} \int_A a^T.a \, dA \, dz \right).\ddot{\mathbf{U}} = \left( \rho.\int_{-t/2}^{+t/2} \int_A (z.\mathbf{H_x}).(z.\mathbf{H_x}^T) + (z.\mathbf{H_y}).(z.\mathbf{H_y}^T) + \mathbf{d}^T).\mathbf{d} \, dA \, dz \right).\ddot{\mathbf{U}} \]
Therefore, the mass matrix is:
\[ M = \rho.\int_{-t/2}^{+t/2} \int_A (z.\mathbf{H_x}).(z.\mathbf{H_x}^T) \, dA \, dz + \rho.\int_{-t/2}^{+t/2} \int_A (z.\mathbf{H_y}).(z.\mathbf{H_y}^T) \, dA \, dz + \rho.\int_{-t/2}^{+t/2} \int_A \mathbf{d}^T.\mathbf{d} \, dA \, dz \]
It clearly appears that the triangle plate mass matrix is composed of 3 matrices, each providing the contribution of a displacement (i.e. \((u, v, w)\)) associated to a given degree of freedom (i.e. respectively \((\theta_y, \theta_x, z)\)). Let's note these 3 matrices \(M_u\), \(M_v\) and \(M_w\).
To integrate over the triangle area, we operate the following variables change:
\[ \left\{ \begin{array}{ccc} x(\xi, \eta) & = & x_0 + \xi.(x_1 - x_0) + \eta.(x_2 - x_0) \\ y(\xi, \eta) & = & y_0 + \xi.(y_1 - y_0) + \eta.(y_2 - y_0) \end{array} \right. \]
Therefore, we have \( J = \left( \begin{array}{cc} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \end{array} \right) = \left( \begin{array}{cc} (x_1 - x_0) & (x_2 - x_0) \\ (y_1 - y_0) & (y_2 - y_0) \\ \end{array} \right) \) and \(|J| = 2A\). Therefore the mass matrices become (noting that the shape functions \(d_i\) do not depend on \(z\)):
\[ \begin{array}{ccc} M_u & = & \rho.\int_{-t/2}^{+t/2} \int_A (z.\mathbf{H_x}).(z.\mathbf{H_x}^T) \, dA \, dz = \rho.\int_{-t/2}^{+t/2} z^2 \, dz.\int_A \mathbf{H_x}.\mathbf{H_x}^T \, dA \\ & = & \rho.\frac{t^3}{12}.\int_0^1 \int_0^{1-\eta} \mathbf{H_x}.\mathbf{H_x}^T |J| \, d\xi d\eta = 2.A.\rho.\frac{t^3}{12}.\int_0^1 \int_0^{1-\eta} \mathbf{H_x}.\mathbf{H_x}^T \, d\xi d\eta \\ M_v & = & \rho.\int_{-t/2}^{+t/2} \int_A (z.\mathbf{H_y}).(z.\mathbf{H_y}^T) \, dA \, dz = \rho.\int_{-t/2}^{+t/2} z^2 \, dz.\int_A \mathbf{H_y}.\mathbf{H_y}^T \, dA \\ & = & \rho.\frac{t^3}{12}.\int_0^1 \int_0^{1-\eta} \mathbf{H_y}.\mathbf{H_y}^T |J| \, d\xi d\eta = 2.A.\rho.\frac{t^3}{12}.\int_0^1 \int_0^{1-\eta} \mathbf{H_y}.\mathbf{H_y}^T \, d\xi d\eta \\ M_w & = & \rho.\int_{-t/2}^{+t/2} \int_A \mathbf{d}^T.\mathbf{d} \, dA \, dz = 2.A.\rho.t.\int_0^1 \int_0^{1-\eta} \mathbf{d}^T.\mathbf{d} \, d\xi d\eta \end{array} \]
To evaluate the 3 integral terms
\[ \begin{array}{ccccc} \int_0^1 \int_0^{1-\eta} \mathbf{H_x}.\mathbf{H_x}^T \, d\xi d\eta & , & \int_0^1 \int_0^{1-\eta} \mathbf{H_y}.\mathbf{H_y}^T \, d\xi d\eta & , & \int_0^1 \int_0^{1-\eta} \mathbf{d}^T.\mathbf{d} \, d\xi d\eta \end{array} \]
we used a formal mathematical package and tested the corresponding code against a numerically exact evaluation using the appropriate Gauss quadrature on the triangle.
The resulting triangle plate mass matrix needs to be assembled properly in the full shell element mass matrix.
|
protectedvirtual |
Computes the triangle's local stiffness matrix.
state | The state to compute the local stiffness matrix from | |
[out] | localStiffnessMatrix | The local stiffness matrix to store the result into |
|
protected |
Computes the triangle's mass matrix.
state | The state to compute the mass matrix from | |
[out] | massMatrix | The mass matrix to store the result into |
|
overridevirtual |
Computes a natural coordinate given a global coordinate.
state | The state at which to transform coordinates |
cartesianCoordinate | The coordinates to transform |
Implements SurgSim::Physics::FemElement.
|
protected |
Computes the triangle element's rotation given a state.
state | The state to compute the rotation from |
|
protected |
Compute the various shape functions (membrane and plate deformations) parameters.
restState | the rest state to compute the shape functions parameters from |
|
protected |
Computes the triangle's stiffness matrix.
state | The state to compute the stiffness matrix from | |
[out] | stiffnessMatrix | The stiffness matrix to store the result into |
|
overrideprotectedvirtual |
Update the FemElement based on the given state.
state | \((x, v)\) the current position and velocity to evaluate the various terms with |
options | Flag to specify which of the F, M, D, K needs to be updated |
Implements SurgSim::Physics::FemElement.
double SurgSim::Physics::Fem2DElementTriangle::getThickness | ( | ) | const |
Gets the triangle's thickness.
|
overridevirtual |
Gets the element volume based on the input state (in m-3)
state | The state to compute the volume with |
Implements SurgSim::Physics::FemElement.
|
overridevirtual |
Initialize the FemElement once everything has been set.
state | The state to initialize the FemElement with |
Reimplemented from SurgSim::Physics::FemElement.
|
protected |
Initializes variables needed before Initialize() is called.
void SurgSim::Physics::Fem2DElementTriangle::setThickness | ( | double | thickness | ) |
Sets the triangle's thickness.
thickness | The thickness of the triangle |
|
protected |
Batoz variable \(a_{k} = -x_{ij}/l_i^2\).
|
protected |
Batoz variable \(b_{k} = 3/4 x_{ij} y_{ij}/l_{ij}2\).
|
protected |
Batoz variable \(c_{k} = (1/4 x_{ij}^2 - 1/2 y_{ij}^2)/l_{ij}^2\).
|
protected |
Batoz variable \(d_{k} = -y_{ij}/l_{ij}^2\).
|
protected |
Batoz variable \(e_{k} = (1/4 y_{ij}^2 - 1/2 x_{ij}^2)/l_{ij}^2\).
|
protected |
Initial rotation matrix for the element.
|
protected |
Plate mass matrix: integral terms related to the dof \((z)\).
|
protected |
Plate mass matrix: integral terms related to the dof \((\theta_y)\).
|
protected |
Plate mass matrix: integral terms related to the dof \((\theta_x)\).
|
protected |
Stiffness matrix (in local coordinate frame)
|
protected |
Batoz variable \(l_{ij}^2 = x_{ij}^2 + y_{ij}^2\).
|
protected |
Membrane (in-plane) deformation.
DOF simulated: (x, y). "Theory of Matrix Structural Analysis" from J.S. Przemieniecki. Shape functions are \(f_i(x, y) = a_i + b_i.x + c_i.y\), here we store \((a_i, b_i, c_i)\) on each row.
|
protected |
Stiffness matrix (in local coordinate frame)
|
protected |
Batoz variable \(P_{k} = -6x_{ij}/l_{ij}^2 = 6.\textrm{m_a}_k\).
|
protected |
Batoz variable \(q_{k} = 3x_{ij}y_{ij}/l_{ij}^2 = 4.\textrm{m_b}_k\).
|
protected |
The triangle rest area.
|
protected |
Batoz variable \(r_{k} = 3y_{ij}^2/l_{ij}^2\).
|
protected |
Thickness of the element.
|
protected |
Batoz variable \(t_{k} = -6y_{ij}/l_{ij}^2 = 6.\textrm{m_d}_k\).
|
protected |
The element's rest state.
|
protected |
Batoz variable \(x_{ij} = x_i - x_j\).
|
protected |
Batoz variable \(y_{ij} = y_i - y_j\).