Mathematics behind the Hybrid Continuum-Discontinuum Method#
In this document, the following conventions are assumed for the notations: (a) Gibbs notation is employed for tensor algebra and the order is denoted by the number of underdots (e.g. \(\underset{˙}{\mathbf{u}} = u_{i}\mathbf{e}_{\mathbf{i}}\),\(\underset{¨}{\mathbf{\sigma}} = \sigma_{ij}\mathbf{e}_{\mathbf{i}} \otimes \mathbf{e}_{\mathbf{j}}\)); (b) the dot symbol means the contraction between two tensors (e.g. \(\underset{˙}{\mathbf{T}} = \underset{¨}{\mathbf{\sigma}} \cdot \underset{˙}{\mathbf{n}}\)); (3) the colon symbol denotes for the double contraction of two second-order tensors (e.g. \(\underset{¨}{\mathbf{\sigma}} = \underset{-}{\mathbf{D}} : \underset{¨}{\mathbf{\varepsilon}}\)); (4) the wedge symbol \(\land\) is for the cross product (e.g. \(\mathbf{e}_{3} = \mathbf{e}_{1} \land \mathbf{e}_{2}\)) and (5) the tensorial product \(\underset{¨}{\mathbf{\sigma}} = \underset{˙}{\mathbf{x}} \otimes \underset{˙}{\mathbf{F}}\) states for the linear production two tensors.
2.1 Weak formulation of updated Lagrange#
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\), \(\underset{˙}{\mathbf{x}}\) is the current position vector of the point and \(\underset{˙}{\mathbf{u}}\left( \underset{˙}{\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 \(\underset{¨}{\mathbf{\varepsilon}}\left( \underset{˙}{\mathbf{x}}\ ,\ t \right) \in \mathbb{R}^{D \times D}\).
Balance of mass:
in \(\forall\underset{˙}{\mathbf{x}} \in \Omega\)
Balance of momentum:
\(\frac{d}{dt}\int_{\Omega}^{}\rho\underset{˙}{\dot{\mathbf{u}}}\left( \underset{˙}{\mathbf{x}} \right)d\Omega = \int_{\Omega}^{}\rho\underset{˙}{\ddot{\mathbf{u}}}\left( \underset{˙}{\mathbf{x}} \right)d\Omega = \int_{\partial\Omega}^{}\underset{˙}{\mathbf{t}}d\partial\Omega + \int_{\Omega}^{}\underset{˙}{\mathbf{b}}d\Omega\) in \(\forall\underset{˙}{\mathbf{x}} \in \Omega\)
\(\rho\) is the density of the solid body. The force is split into tractions \(\underset{˙}{\mathbf{t}}\) acted on boundaries \(\partial\Omega\)and into body force \(\underset{˙}{\mathbf{b}}\) in volume \(\Omega \mapsto \mathbb{R}^{D}\). Similarly, the conservation of the moment of momentum can be written as
\(\int_{\Omega}^{}\rho\underset{˙}{\ddot{\mathbf{u}}}\left( \underset{˙}{\mathbf{x}} \right) \land \underset{˙}{\mathbf{x}}d\Omega = \int_{\partial\Omega}^{}\underset{˙}{\mathbf{t}} \land \underset{˙}{\mathbf{x}}d\partial\Omega + \int_{\Omega}^{}{\underset{˙}{\mathbf{b}} \land \underset{˙}{\mathbf{x}}}d\Omega\) in \(\forall\underset{˙}{\mathbf{x}} \in \Omega\)
In terms of the Cauchy stress tensor, the traction force satisfies \(\underset{˙}{\mathbf{t}} = \underset{¨}{\mathbf{\sigma}} \cdot \underset{˙}{\mathbf{n}}\), Eq is rewritten according to Gausss theorem \(\int_{\partial\Omega}^{}\underset{˙}{\mathbf{t}}d\partial\Omega = \int_{\partial\Omega}^{}{\underset{¨}{\mathbf{\sigma}} \cdot \underset{˙}{\mathbf{n}}}d\partial\Omega = \int_{\Omega}^{}{\underset{¨}{\mathbf{\sigma}} \cdot \nabla}d\Omega\). Since the stress tensor is symmetric based on Eq as \(\underset{¨}{\mathbf{\sigma}} = {\underset{¨}{\mathbf{\sigma}}}^{T}\), the conservation of momentum is
\(\rho\underset{˙}{\ddot{\mathbf{u}}}\left( \underset{˙}{\mathbf{x}} \right) - \nabla \cdot \underset{¨}{\mathbf{\sigma}} - \underset{˙}{\mathbf{b}} = 0\) in \(\forall\underset{˙}{\mathbf{x}} \in \Omega\)
Boundary conditions:
\(\underset{¨}{\mathbf{\sigma}} \cdot \underset{˙}{\mathbf{n}} = \underset{˙}{\mathbf{t}}\) on \(\forall\underset{˙}{\mathbf{x}} \in \partial\Omega_{t}\)
\(\underset{˙}{\mathbf{u}} = \overline{\underset{˙}{\mathbf{u}}}\) on \(\forall\underset{˙}{\mathbf{x}} \in \partial\Omega_{u}\)
Initial conditions:
\(\underset{˙}{\mathbf{u}}\left( \underset{˙}{\mathbf{x}}\ ,\ t\ = \ 0 \right) = \overline{\underset{˙}{\mathbf{u}}}\left( \underset{˙}{\mathbf{x}} \right)\) on \(\forall\underset{˙}{\mathbf{x}} \in \partial\Omega_{u}\)
\(\dot{\underset{˙}{\mathbf{u}}}\left( \underset{˙}{\mathbf{x}}\ ,\ t\ = \ 0 \right) = \overline{\dot{\underset{˙}{\mathbf{u}}}}\left( \underset{˙}{\mathbf{x}} \right)\) on \(\forall\underset{˙}{\mathbf{x}} \in \partial\Omega_{u}\)
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
\(\int_{\Omega}^{}{\left( \left( \nabla\ \cdot \ \underset{¨}{\mathbf{\sigma}} \right)\ \cdot \ \delta\ \underset{˙}{\mathbf{v}}\ - \ \left( \underset{˙}{\mathbf{b}}\ - \ \rho\ \underset{˙}{\ddot{\mathbf{u}}} \right)\ \cdot \ \delta\ \underset{˙}{\mathbf{v}} \right)d\Omega} = 0\) in \(\forall\underset{˙}{\mathbf{v}} \in \Omega\)
Integration by parts of the first term in Eq gives
\(\int_{\Omega}^{}\left( \nabla\ \cdot \ \underset{¨}{\mathbf{\sigma}} \right) \cdot \delta\underset{˙}{\mathbf{v}}d\Omega = \int_{\Omega}^{}\left( \underset{¨}{\mathbf{\sigma}}\ : \ \nabla\ \delta\ \underset{˙}{\mathbf{v}} \right)d\Omega - \int_{\partial\Omega}^{}\left( {\underset{¨}{\mathbf{\sigma}}}^{T}\ \delta\ \underset{˙}{\mathbf{v}}\ \cdot \ \underset{˙}{\mathbf{n}} \right)d\partial\Omega = \int_{\Omega}^{}\left( \underset{¨}{\mathbf{\sigma}}\ : \ \nabla\ \delta\ \underset{˙}{\mathbf{v}} \right)d\Omega - \int_{\partial\Omega_{t}}^{}\left( \underset{˙}{\mathbf{t}}\ \cdot \ \delta\ \underset{˙}{\mathbf{v}} \right)d\partial\Omega\) in \(\forall\underset{˙}{\mathbf{v}} \in \Omega\)
Substitution of Eq. back into Eq gives
\(\int_{\Omega}^{}{\left( \underset{¨}{\mathbf{\sigma}}\ : \ \nabla\ \delta\ \underset{˙}{\mathbf{v}}\ - \ \left( \underset{˙}{\mathbf{b}}\ - \ \rho\ \underset{˙}{\ddot{\mathbf{u}}} \right)\ \cdot \ \delta\ \underset{˙}{\mathbf{v}} \right)d\Omega = \int_{\partial\Omega_{t}}^{}{\underset{˙}{\mathbf{t}} \cdot}}\delta\underset{˙}{\mathbf{v}}d\partial\Omega\) in \(\forall\underset{˙}{\mathbf{v}} \in \Omega\)
2.2 Spatial discretization in finite domain#
In this section, Eq 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( \underset{˙}{\mathbf{x}} \right)\) associated with node \(i\) having coordinate \(\underset{˙}{\mathbf{x}}\). The global interpolation matrix is defined as [57]
\(\mathcal{N}^{g}\left( \underset{˙}{\mathbf{x}} \right) = \left\lbrack diag\ \left\lbrack \mathcal{N}_{1}^{g}\ \left( \underset{˙}{\mathbf{x}} \right) \right\rbrack\ \ diag\ \left\lbrack \mathcal{N}_{2}^{g}\ \left( \underset{˙}{\mathbf{x}} \right) \right\rbrack\ \cdots\ \ diag\ \left\lbrack \mathcal{N}_{n_{poin}}^{g}\ \left( \underset{˙}{\mathbf{x}} \right) \right\rbrack \right\rbrack\)
where \(diag\left\lbrack \mathcal{N}_{1}^{g}\ \left( \underset{˙}{\mathbf{x}} \right) \right\rbrack\) is the \(n_{\dim}\)×\(n_{\dim}\) diagonal matrix. The element displacement in global coordinate is
for \(\forall\underset{˙}{\mathbf{x}} \in \Omega^{e}\) \({}^{h}\underset{˙}{\mathbf{u}}\left( \underset{˙}{\mathbf{x}} \right) = \sum_{i = 1}^{n_{poi}}{\mathcal{N}_{i}\left( \underset{˙}{\mathbf{x}} \right)}\underset{˙}{\mathbf{u}} = \mathcal{N}^{g}\underset{˙}{\mathbf{u}}\)
Therefore, the discretised virtual work expression considering the kinetic is
\(\int_{{}^{h}\Omega}^{}\left\lbrack {\underset{¨}{\mathbf{\sigma}}}^{T}\ \mathcal{B}^{g}\ \delta\ \underset{˙}{\mathbf{v}}\ - \ \left( \underset{˙}{\mathbf{b}}\ - \ \rho\ \underset{˙}{\ddot{\mathbf{u}}} \right)\ \cdot \ \mathcal{N}^{g}\ \delta\ \underset{˙}{\mathbf{v}} \right\rbrack d\Omega - \int_{\partial^{h}\Omega_{t}}^{}\underset{˙}{\mathbf{t}} \cdot \mathcal{N}^{g}\delta\underset{˙}{\mathbf{v}}d\partial\Omega = 0\) in \(\forall\delta\underset{˙}{\mathbf{v}} \in \Omega\)
Since the above equation is satisfied for all vectors \(\delta\underset{˙}{\mathbf{v}}\), equation can be reformulated as
\(\mathbf{M}^{FE}\underset{˙}{\ddot{\mathbf{u}}} + {\underset{˙}{\mathbf{f}}}^{int}\left( \underset{˙}{\mathbf{u}} \right) - {\underset{˙}{\mathbf{f}}}^{ext} = 0\)
while the internal force \({\underset{˙}{\mathbf{f}}}^{int}\left( \underset{˙}{\mathbf{u}} \right)\), external force \({\underset{˙}{\mathbf{f}}}^{ext}\)and inertia mass \(\mathbf{M}^{FE}\) are
\({\underset{˙}{\mathbf{f}}}^{int}\left( \underset{˙}{\mathbf{u}} \right) = \int_{{}^{h}\Omega}^{}{\left( \mathcal{B}^{g} \right)^{T} \cdot \underset{¨}{\mathbf{\sigma}}\left( \underset{˙}{\mathbf{u}} \right)d\Omega}\)
\({\underset{˙}{\mathbf{f}}}^{ext} = \int_{{}^{h}\Omega}^{}{\left( \mathcal{N}^{g} \right)^{T}\rho\underset{˙}{\mathbf{b}}d\Omega} + \int_{\partial^{h}\Omega_{t}}^{}{\left( \mathcal{N}^{g} \right)^{T}\underset{˙}{\mathbf{t}}}d\partial\Omega = 0\)
\(\mathbf{M}^{FE} = \int_{{}^{h}\Omega}^{}{\left( \mathcal{N}^{g} \right)^{T}\rho\mathcal{N}^{g}d\Omega}\)
\(\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
\(\begin{array}{r} {\underset{˙}{\mathbf{f}}}^{int} = \overset{e = 1}{\underset{n_{elem}}{\mathbf{A}}}\left( {\underset{˙}{\mathbf{f}}}_{(e)}^{int} \right) \\ {\underset{˙}{\mathbf{f}}}^{ext} = \overset{e = 1}{\underset{n_{elem}}{\mathbf{A}}}\left( {\underset{˙}{\mathbf{f}}}_{(e)}^{ext} \right) \end{array}\)
where \(\overset{e = 1}{\underset{n_{elem}}{\mathbf{A}}}\) is the finite element assembly operator and the element force vector can be obtained as
\({\underset{˙}{\mathbf{f}}}_{(e)}^{int}\left( \underset{˙}{\mathbf{u}} \right) = \int_{{}^{e}\Omega}^{}{\left( \mathcal{B}^{e} \right)^{T} \cdot \underset{¨}{\mathbf{\sigma}}\left( \underset{˙}{\mathbf{u}} \right)d\Omega}\)
\({\underset{˙}{\mathbf{f}}}_{(e)}^{ext} = \int_{{}^{e}\Omega}^{}{\left( \mathcal{N}^{e} \right)^{T}\rho\underset{˙}{\mathbf{b}}d\Omega} + \int_{\partial^{e}\Omega_{t}}^{}{\left( \mathcal{N}^{e} \right)^{T}\underset{˙}{\mathbf{t}}}d\partial\Omega = 0\)
The superscripts \(e\) denotes the variables in local element.
2.3 Dynamic relaxation and time integration scheme#
The continuum-discontinuum method (CDM) employs the explicit time integration based on velocity verlet algorithm to solve Eq . \({}^{n + 1}\dot{\underset{˙}{\mathbf{u}}}\) and \({}^{n}\underset{˙}{\mathbf{u}}\) are the velocity vector and displacement vector at \({}^{n + 1}t\) and \({}^{n}t\) for a time step \(\mathrm{\Delta}t\)
\({}^{n}\ddot{\underset{˙}{\mathbf{u}}} = \frac{1}{2\mathrm{\Delta}t}\left( {}^{n + 1}\ \dot{\underset{˙}{\mathbf{u}}}\ -^{n}\ \dot{\underset{˙}{\mathbf{u}}} \right)\)
\({}^{n}\dot{\underset{˙}{\mathbf{u}}} =^{n + 1}\underset{˙}{\mathbf{u}} -^{n}\underset{˙}{\mathbf{u}}\)
If the damping matrix \(\mathbf{C}^{FE} = \alpha\mathbf{M}^{FE}\) is considered, the governing equation is written as
\(\mathbf{M}^{FE}\underset{˙}{\ddot{\mathbf{u}}} + \mathbf{C}^{FE}\dot{\underset{˙}{\mathbf{u}}} + {\underset{˙}{\mathbf{f}}}^{int}\left( \underset{˙}{\mathbf{u}} \right) - {\underset{˙}{\mathbf{f}}}^{ext} = 0\)
Eq. can be rewritten in the form of nodal velocities
\({}^{n + 1}\underset{˙}{\dot{\mathbf{u}}} = \left\lbrack \left( 1\ - \ 2\ \alpha\ \mathrm{\Delta}t \right)^{n}\ \underset{˙}{\dot{\mathbf{u}}}\ + \ \frac{2\mathrm{\Delta}t}{\mathbf{M}^{FE}}\ \left( {}^{n}\ {\underset{˙}{\mathbf{f}}}^{ext}\ -^{n}\ {\underset{˙}{\mathbf{f}}}^{int}\ \left( {}^{n}\ \underset{˙}{\mathbf{u}} \right) \right) \right\rbrack\)
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
\(\mathrm{\Delta}t_{cri} = \gamma\frac{2}{\omega}\left( \sqrt{\xi^{2} + 1}\ - \ \xi \right)\)
In CDM, each triangular element is considered as elastic with constant strain. During each timestep, the element strain rate is computed as
\(\dot{\underset{¨}{\mathbf{\varepsilon}}} = \frac{1}{2}\left( \nabla\ \dot{\underset{˙}{\mathbf{u}}}\ + \ \nabla^{T}\ \dot{\underset{˙}{\mathbf{u}}} \right),\dot{\underset{¨}{\mathbf{\theta}}} = \frac{1}{2}\left( \nabla\ \dot{\underset{˙}{\mathbf{u}}}\ - \ \nabla^{T}\ \dot{\underset{˙}{\mathbf{u}}} \right),\nabla\dot{\underset{˙}{\mathbf{u}}} \cong \frac{1}{2A}\sum_{s}^{}{\left( {\dot{\underset{˙}{\mathbf{u}}}}^{(a)}\ + \ {\dot{\underset{˙}{\mathbf{u}}}}^{(b)} \right) \cdot \mathbf{n}\mathrm{\Delta}s}\)
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
\(\mathrm{\Delta}\underset{¨}{\mathbf{\sigma}} = \lambda{\dot{\underset{¨}{\mathbf{\varepsilon}}}}_{v}\mathbf{\delta}\mathrm{\Delta}t + 2\mu\dot{\underset{¨}{\mathbf{e}}}\mathrm{\Delta}t\)
where \(\mathbf{\delta}\) is Kronecker symbol, \(\lambda\) and \(G\) are Lame constants. The updated stress and strain at \(n\) step are
\({}^{n + 1}\underset{¨}{\mathbf{\sigma}} =^{n}\underset{¨}{\mathbf{\sigma}} + \mathrm{\Delta}\underset{¨}{\mathbf{\sigma}},{}^{n + 1}\dot{\underset{¨}{\mathbf{\varepsilon}}} =^{n}\dot{\underset{¨}{\mathbf{\varepsilon}}} + \mathrm{\Delta}\underset{¨}{\mathbf{\varepsilon}}\)
2.4 Contact algorithm in CDM#
The contact algorithms in CDM is divided into: (a) neighbour search and (b) geometric resolution. Neighbour search is a rough search aims to provide a list of possible blocks in contact including global search e.g. an altering digital tree (ADT) [58], no binary search (NBS) contact detection algorithm [59], CGRID [60], DESS algorithm [61], master-slave algorithm[62] and local search, e.g. node-to-surface (NTS) algorithm[63], RIGID algorithm[64], inside-outside algorithm [65]. Geometric resolution algorithms strongly depend on complexity of the geometric representation of cells. The discrete function representation (DFR) algorithm has a computational complexity of order O(N) [66] and the common-plane (CP), which bisects the space between two contacting blocks, is advantageous for vertex-to-vertex or edge-to-vertex contacts where the definition of the contact normal is a non-trivial problem and the method has a complexity of order O(N) [67]. The number of iterations of CP is depended on the accuracy of the initial guess of the CP position, and then a fast common plane (FCP) method [68] and a shortest link (SL) method [69] are proposed to increase efficiency of the CP method.
In this study, the contact detection algorithm is developed for blocks with evident size using square bounding box method in cooperation with the NBS algorithm (shown in Figure 2), which aims to divide the NBS model to several groups (list in Table 1), the flow chart of the enhanced NBS algorithm is
\({\underset{˙}{\mathbf{x}}}^{k} = 1 + int\left( \frac{{\underset{˙}{\mathbf{x}}}_{i} - {\underset{˙}{\mathbf{x}}}_{\min}}{d(0)}\ + \ \frac{1}{2} \right)\)
Figure 2. The schematic algorithm of NBS contact detection in CDM.
The enhanced NBS contact detection method used in mGbCDM
Step 1: Loop the blocks and find the maximum size buffer box for the initial group box \(d(0) = d_{\max}\); |
|---|
Step 2: Divide the blocks into \(n\) groups with size of buffer box for the nth group box as \(d(n) = d_{\max} \cdot \alpha^{n - 1}\), \(\alpha \in \left( 0\ ,\ \left. \ 1 \right\rbrack \right.\ \); |
Step 3: All the blocks are mapped in to the grid space with edge length of \(d(0)\) as depicted in Figure 2 , the central point of the block \({\underset{˙}{\mathbf{x}}}^{k}\) is computed in Eq ; |
Step 4: Loop all the blocks and detect contacts for the first group, the contact couple groups is identified when \(\left| {\underset{˙}{\mathbf{x}}}^{(t)}\ - \ {\underset{˙}{\m athbf{x}}}^{(c)} \right| < \max\left( d^{(t)}\ + \ d^{(c)} \right)\), the contact state can be recognized as neighboring contacts or center contacts; |
Step 5: Repeat step 3 and step 4 for all groups of the remaining blocks and identify the states of the contacts; |