In this document, the following conventions are assumed for the notations:
Gibbs notation is employed for tensor algebra and the order is denoted by the number of underdots (e.g. \({\mathbf{u}} = u_{i}\mathbf{e}_{\mathbf{i}}\),\({\mathbf{\sigma}} = \sigma_{ij}\mathbf{e}_{\mathbf{i}} \otimes \mathbf{e}_{\mathbf{j}}\));
the dot symbol means the contraction between two tensors (e.g. \({\mathbf{T}} = {\mathbf{\sigma}} \cdot {\mathbf{n}}\));
the colon symbol denotes for the double contraction of two second-order tensors (e.g. \({\mathbf{\sigma}} = {\mathbf{D}} : {\mathbf{\varepsilon}}\));
the wedge symbol \(\land\) is for the cross product (e.g. \(\mathbf{e}_{3} = \mathbf{e}_{1} \land \mathbf{e}_{2}\)) and,
the tensorial product \({\mathbf{\sigma}} = {\mathbf{x}} \otimes {\mathbf{F}}\) states for the linear production two tensors.
Consider an arbitrary body \(\Omega \subset \mathbb{R}^{D}\left( D\ \in \ [\ 1\ ,\ 2\ ,\ 3\ ] \right)\)
with external boundary \(\partial\Omega\) and the internal
discontinuity boundary \(\Gamma\), \({\mathbf{x}}\)
is the current position vector of the point and
\({\mathbf{u}}\left( {\mathbf{x}}\ ,t \right) \subset \mathbb{R}^{D}\)
is the displacement at time
\(t \in \left\lbrack 0\ ,\ t_{all} \right\rbrack\). The
infinitesimal strain tensor based on small deformation and deformation
gradients assumptions is characterized as
\({\mathbf{\varepsilon}}\left( {\mathbf{x}}\ ,\ t \right) \in \mathbb{R}^{D \times D}\).
where \(\rho\) is the density of the solid body. The force is split into
tractions \({\mathbf{t}}\) acted on boundaries
\(\partial\Omega\) and into body force
\({\mathbf{b}}\) in volume
\(\Omega \mapsto \mathbb{R}^{D}\). Similarly, the conservation of
the moment of momentum can be written as
In terms of the Cauchy stress tensor, the traction force satisfies
\({\mathbf{t}} = {\mathbf{\sigma}} \cdot {\mathbf{n}}\),
The equation above is rewritten according to Gausss theorem
\(\int_{\partial\Omega}^{}{\mathbf{t}}d\partial\Omega = \int_{\partial\Omega}^{}{{\mathbf{\sigma}} \cdot {\mathbf{n}}}d\partial\Omega = \int_{\Omega}^{}{{\mathbf{\sigma}} \cdot \nabla}d\Omega\).
Since the stress tensor is symmetric as
\({\mathbf{\sigma}} = {{\mathbf{\sigma}}}^{T}\),
the conservation of momentum is
The spatial version of the principle of virtual work states that the body \(\Omega\) is in equilibrium when the Cauchy stress satisfies the initial condition
In this section, Eq.8 is discretised to derive the mass matrix, internal
force and external force which is convenient to the matrix notation. Let
the generic finite element \(e \in \Omega\) is defined by
\(n_{node}\) nodes with shape function
\(\mathcal{N}^{(e)}\left( {\mathbf{x}} \right)\)
associated with node \(i\) having coordinate
\({\mathbf{x}}\). The global interpolation matrix is
defined as
where
\(diag\left\lbrack \mathcal{N}_{1}^{g}\ \left( {\mathbf{x}} \right) \right\rbrack\)
is the \(n_{\dim}\)×\(n_{\dim}\) diagonal matrix. The
element displacement in global coordinate is
while the internal force
\({{\mathbf{f}}}^{int}\left( {\mathbf{u}} \right)\),
external force \({{\mathbf{f}}}^{ext}\)and inertia
mass \(\mathbf{M}^{FE}\) are
\(\mathcal{B}^{g}\) and \(\mathcal{N}^{g}\) are matrices
incorporating the interpolation functions and their spatial derivatives
in global, respectively. The actual force vectors are assembled as
The continuum-discontinuum method (CDM) employs the explicit time
integration based on velocity verlet algorithm to solve Eq .12,
\({}^{n + 1}\dot{{\mathbf{u}}}\) and
\({}^{n}{\mathbf{u}}\) are the velocity vector and
displacement vector at \({}^{n + 1}t\) and \({}^{n}t\) for a
time step \(\mathrm{\Delta}t\)
while the damping coefficient is \(\alpha = 2\xi\omega\),
\(\xi\) is the damping ratio and \(\omega\) is the frequency
parameter. The stability of the explicit scheme is depended on the
critical timestep, which is determined by
The incremental strain tensor is denoted as
\(\mathrm{\Delta}\varepsilon_{ij} = {\dot{\varepsilon}}_{ij}\mathrm{\Delta}t\),
\(\mathbf{n}\) is unit vector in normal direction, \(A\) is the
area and \(s\) is the length of edge, the incremental stress is
updated according to the constitutive relations of the isotropic elastic
blocks as
The classical FDEM framework is limited to linear elements, as both its
cohesive element formulation and contact algorithm rely on the linear-edge
assumption. OpenFDEM generalizes it to an arbitrary higher-order element
framework. For simplicity, the quadratic triangular element is used as a
primary example here; the shape functions of other higher-order elements can
be constructed similarly.
The shape functions serve as a fundamental mathematical basis for
interpolating field variables within a bulk element from nodal values. Shape
functions of arbitrary order of triangular and rectangular elements can be
constructed through approaches based on the barycentric coordinates method and
the Lagrange polynomials method.
For a triangular element of arbitrary order \(p\), the shape functions
have a universal form as
where the subscript denotes the corresponding nodes shown in Figure 2. It is
simple to employ quadratic order shape functions to compute the internal nodal
force of the bulk element (Eq. 18) and the nodal external force (Eq. 19), in
the same way used in the finite element method. The strain matrix
\(\mathcal{B}\) of quadratic order in natural coordinates can also be
derived accordingly.
