An improved updated Lagrangian SPH method for structural modelling

Abstract

This paper presents a set of novel refined schemes to enhance the accuracy and stability of the updated Lagrangian SPH (ULSPH) for structural modelling. The original ULSPH structure model was first proposed by Gray et al. (Comput Methods Appl Mech Eng 190:6641–6662, 2001) and has been utilised for a wide range of structural analyses including metal, soil, rubber, ice, etc., although the model often faces several drawbacks including unphysical numerical damping, high-frequency noise in reproduced stress fields, presence of several artificial terms requiring ad hoc tunings and numerical instability in the presence of tensile stresses. In these regards, this study presents a set of enhanced schemes corresponding to (1) consistency correction on discretisation schemes for differential operators, (2) a numerical diffusive term incorporated in the continuity or the density rate equation, (3) tuning-free stabilising term based on Riemann solution and (4) careful control/switch of stress divergence differential operator model under tensile stresses. Qualitative/quantitative validations are conducted through several well-known benchmark tests.

Graphical abstract

Appendices

Appendix A. Theoretical solution of elastic wave propagation (Sect. 4.5)

In this appendix, the analytical solution of 1D elastic wave propagation corresponding to the benchmark test in Sect. 4.5 is derived based on the elementary rod theory [52]. In this theory, the rod is assumed to be long and slender so that the lateral contraction could be negligibly small.

The equation of motion is obtained as follows. Let us consider an element inside a rod. Let q(x,t) be the external force per unit volume, u(x,t) be the displacement in x-direction, A be the area of cross section, Δx be the length of the element in x-direction as an infinitesimal value, ρ be the density of the rod, and F be the force acting on the face of cross section, as shown in Fig. 

Fig. 34

An element of a rod with subjected loads

34.

Considering the force balance for an element of the rod as:

$$ \left( {\rho A\Delta x} \right)\frac{{\partial^{2} u}}{{\partial t^{2} }} = - F + \left( {F + \frac{\partial F}{{\partial x}}\Delta x} \right) + qA\Delta x $$

(A.1)

where ρA denotes the mass per unit length of the rod. The above momentum equation can be expressed as:

$$ \frac{\partial F}{{\partial x}} = \rho A\frac{{\partial^{2} u}}{{\partial t^{2} }} - qA $$

(A.2)

The relationship between the strain ε and displacement u can be written as:

$$ \varepsilon = \frac{\partial u}{{\partial x}} $$

(A.3)

The Hooke’s law in one-dimensional form can be described as:

$$ \sigma = \frac{F}{A} = E\varepsilon $$

(A.4)

where $\sigma$ stands for stress. Substituting Eqs. (A.3) and (A.4) into Eq. (A.2), we obtain Eq. (A.5) as:

$$ \frac{\partial }{\partial x}\left( {EA\frac{\partial u}{{\partial x}}} \right) = \rho A\frac{{\partial^{2} u}}{{\partial t^{2} }} - qA $$

(A.5)

If the external force q(x,t) does not exist and the rod is a homogeneous material, we can obtain the below wave equation:

$$ \frac{{\partial^{2} u}}{{\partial t^{2} }} - C^{2} \frac{{\partial^{2} u}}{{\partial x^{2} }} = 0;\quad C = \sqrt {\frac{E}{\rho }} $$

(A.6)

where C represents the speed of sound corresponding to the velocity of the propagating wave. We can solve this wave equation under the provided initial conditions and boundary conditions.

For a composite rod, i.e. a rod consisting of physically different materials with clear phase interface, we consider general d’Alembert solution of Eq. (A.6) for each material phase. In case of the test case in Sect. 4.5, left component (phase 1: aluminium) and right component (phase 2: copper) are connected at the centre of the rod with a clear material interface. The wave reflects at the interface due to the discontinuity of materials. Considering that displacement and stress are continuous across the material interface, we can find the reflection coefficient $A_{1}^{{{\text{ref}}}}$ and transmission coefficient $A_{2}^{{{\text{inc}}}}$.

For phase 1, the displacement and the stress induced by the incident wave (coefficient of $A_{1}^{{{\text{inc}}}}$) can be written with using angular frequency ω and wave number k as:

$$ u_{1}^{{{\text{inc}}}} \left( {x,t} \right) = A_{1}^{{{\text{inc}}}} \exp \left( {i\left( {\omega t - k_{1} x} \right)} \right) $$

(A.7)

$$ \sigma_{1}^{{{\text{inc}}}} \left( {x,t} \right) = E_{1}^{{}} \frac{{\partial u_{1}^{{{\text{inc}}}} }}{\partial x} = - A_{1}^{{{\text{inc}}}} E_{1}^{{}} k_{1}^{{}} \exp \left( {i\left( {\omega t - k_{1} x} \right)} \right) $$

(A.8)

The displacement and the stress occurred by reflection wave are obtained as:

$$ u_{1}^{{{\text{ref}}}} \left( {x,t} \right) = A_{1}^{{{\text{ref}}}} \exp \left( {i\left( {\omega t + k_{1} x} \right)} \right) $$

(A.9)

$$ \sigma_{1}^{{{\text{ref}}}} \left( {x,t} \right) = E_{1}^{{}} \frac{{\partial u_{1}^{{{\text{ref}}}} }}{\partial x} = A_{1}^{{{\text{ref}}}} E_{1}^{{}} k_{1}^{{}} \exp \left( {i\left( {\omega t + k_{1} x} \right)} \right) $$

(A.10)

For phase 2, the displacement and the stress can be calculated as:

$$ u_{2}^{{{\text{inc}}}} \left( {x,t} \right) = A_{2}^{{{\text{inc}}}} \exp \left( {i\left( {\omega t - k_{2} x} \right)} \right) $$

(A.11)

$$ \sigma_{2}^{{{\text{inc}}}} \left( {x,t} \right) = E_{2}^{{}} \frac{{\partial u_{2}^{{{\text{inc}}}} }}{\partial x} = - A_{2}^{{{\text{inc}}}} E_{2}^{{}} k_{2}^{{}} \exp \left( {i\left( {\omega t - k_{2} x} \right)} \right) $$

(A.12)

Since the displacement and stress are continuous across the material interface, the following relationship would hold:

$$ A_{1}^{{{\text{inc}}}} + A_{1}^{{{\text{ref}}}} = A_{2}^{{{\text{inc}}}} $$

(A.13)

$$ - A_{1}^{{{\text{inc}}}} E_{1}^{{}} k_{1}^{{}} + A_{1}^{{{\text{ref}}}} E_{1}^{{}} k_{1}^{{}} = - A_{2}^{{{\text{inc}}}} E_{2}^{{}} k_{2}^{{}} $$

(A.14)

Thus, the reflection coefficient and transmission coefficient can be obtained as:

$$ A_{1}^{{{\text{ref}}}} = A_{1}^{{{\text{inc}}}} \frac{{E_{1}^{{}} k_{1}^{{}} - E_{2}^{{}} k_{2}^{{}} }}{{E_{1}^{{}} k_{1}^{{}} + E_{2}^{{}} k_{2}^{{}} }} $$

(A.15)

$$ A_{2}^{{{\text{inc}}}} = 2A_{1}^{{{\text{inc}}}} \frac{{E_{1}^{{}} k_{1}^{{}} }}{{E_{1}^{{}} k_{1}^{{}} + E_{2}^{{}} k_{2}^{{}} }} $$

(A.16)

Since the components would have difference in both Young’s modulus and density, Eqs. (A.15) and (A.16) will be further formulated to the form including density. We now consider the continuity of angular frequency across material interface as:

$$ \omega = k_{1}^{{}} C_{1}^{{}} = k_{2}^{{}} C_{2}^{{}} $$

(A.17)

We obtain the relationship of wave numbers between two phases as:

$$ \frac{{k_{1}^{{}} }}{{k_{2}^{{}} }} = \frac{{C_{2}^{{}} }}{{C_{1}^{{}} }} = \sqrt {\frac{{E_{2}^{{}} \rho_{1}^{{}} }}{{E_{1}^{{}} \rho_{2}^{{}} }}} $$

(A.18)

Equation (A.15) can be rewritten by substituting Eq. (A.18) as:

$$ \begin{aligned} A_{1}^{{{\text{ref}}}} & = A_{1}^{{{\text{inc}}}} \frac{{E_{1}^{{}} k_{1}^{{}} - E_{2}^{{}} k_{2}^{{}} }}{{E_{1}^{{}} k_{1}^{{}} + E_{2}^{{}} k_{2}^{{}} }} \\ & = A_{1}^{{{\text{inc}}}} \frac{{1 - \frac{{E_{2}^{{}} k_{2}^{{}} }}{{E_{1}^{{}} k_{1}^{{}} }}}}{{1 + \frac{{E_{2}^{{}} k_{2}^{{}} }}{{E_{1}^{{}} k_{1}^{{}} }}}} \\ & = A_{1}^{{{\text{inc}}}} \frac{{1 - \sqrt {\frac{{E_{2}^{{}} \rho_{2}^{{}} }}{{E_{1}^{{}} \rho_{1}^{{}} }}} }}{{1 + \sqrt {\frac{{E_{2}^{{}} \rho_{2}^{{}} }}{{E_{1}^{{}} \rho_{1}^{{}} }}} }} \\ & = A_{1}^{{{\text{inc}}}} \frac{{1 - \alpha_{r} }}{{1 + \alpha_{r} }} \\ \end{aligned} $$

(A.19)

$$ \alpha_{r} = \sqrt {\frac{{E_{2}^{{}} \rho_{2}^{{}} }}{{E_{1}^{{}} \rho_{1}^{{}} }}} $$

(A.20)

Substituting Eq. (A.19) into Eq. (A.13), we can obtain the amplitude of the transmission wave as:

$$ \begin{aligned} A_{2}^{{{\text{inc}}}} & = A_{1}^{{{\text{inc}}}} + A_{1}^{{{\text{ref}}}} \\ & = A_{1}^{{{\text{inc}}}} + A_{1}^{{{\text{inc}}}} \frac{{1 - \alpha_{r} }}{{1 + \alpha_{r} }} \\ & = A_{1}^{{{\text{inc}}}} \frac{2}{{1 + \alpha_{r} }} \\ \end{aligned} $$

(A.21)

Appendix B. Discussion on the similarity among Riemann diffusive term, artificial viscosity and δ-SPH

The continuity equation including the δ or the density diffusive term can be written as (also discussed in Sect. 3.2):

$$ \frac{{{\text{D}}\rho_{i}^{{}} }}{{{\text{D}}t}} = - \rho_{i}^{{}} \sum\limits_{j} {{\varvec{u}}_{ij}^{{}} \cdot \nabla_{i}^{{}} w_{ij} V_{j} } + D_{i} $$

(B.1)

$$ D_{i} = \delta hc_{0}^{{}} \sum\limits_{j} {\psi_{ij} \frac{{{\varvec{r}}_{ij} \cdot \nabla_{i}^{{}} w_{ij} }}{{\left| {{\varvec{r}}_{ij} } \right|^{2} }}V_{j} } $$

(B.2)

$$ \psi_{ij}^{{{\text{0th}}}} = 2\rho_{ij} $$

(B.3)

$$ \psi_{ij}^{{{\text{1st}}}} = 2\rho_{ij} - \left( {\left\langle {\nabla \rho } \right\rangle_{i} + \left\langle {\nabla \rho } \right\rangle_{j} } \right) \cdot {\varvec{r}}_{ij} $$

(B.4)

Equations (B.3) and (B.4) correspond to zeroth- and first-order corrected functions (Refs. [12, 58], respectively).

In Riemann SPH, for the case of linear reconstruction of variables, the continuity equation is written as:

$$ \frac{{{\text{D}}\rho_{i}^{{}} }}{{{\text{D}}t}} = - 2\rho_{i}^{{}} \sum\limits_{j} {\left( {{\varvec{u}}_{{}}^{*} - {\varvec{u}}_{i}^{{}} } \right) \cdot \nabla_{i}^{{}} w_{ij} V_{j} } $$

(B.5)

$$ {\varvec{u}}_{{}}^{*} = u_{{}}^{*} \frac{{{\varvec{r}}_{ij} }}{{\left| {{\varvec{r}}_{ij} } \right|}} + \left\{ {\frac{{{\varvec{u}}_{i}^{{}} + {\varvec{u}}_{j}^{{}} }}{2} - \frac{{u_{i}^{{}} + u_{j}^{{}} }}{2}\frac{{{\varvec{r}}_{ij} }}{{\left| {{\varvec{r}}_{ij} } \right|}}} \right\} $$

(B.6)

$$ u_{{}}^{*} = \frac{1}{2}\left[ {u_{i} + u_{j} + \frac{1}{{C_{ij} }}\left( {p_{i} - p_{j} } \right)} \right] $$

(B.7)

$$ u_{i}^{{}} = {\varvec{u}}_{i}^{{}} \cdot \frac{{{\varvec{r}}_{ij} }}{{\left| {{\varvec{r}}_{ij} } \right|}}_{{}} ;\quad u_{j}^{{}} = {\varvec{u}}_{j}^{{}} \cdot \frac{{{\varvec{r}}_{ij} }}{{\left| {{\varvec{r}}_{ij} } \right|}} $$

(B.8)

$$ C_{ij} = \frac{{c_{i} \rho_{i} \sqrt {\rho_{i} } + c_{j} \rho_{j} \sqrt {\rho_{j} } }}{{\sqrt {\rho_{i} } + \sqrt {\rho_{j} } }} $$

(B.9)

From the equation of state, we have:

$$ p_{i}^{{}} = c_{0}^{2} \left( {\rho_{i}^{{}} - \rho_{0}^{{}} } \right) $$

(B.10)

Reformulating equations by substituting Eqs. (B.6) to (B.10) into Eq. (B.5), then:

$$ \frac{{{\text{D}}\rho_{i}^{{}} }}{{{\text{D}}t}} = - \rho_{i}^{{}} \sum\limits_{j} {{\varvec{u}}_{ij}^{{}} \cdot \nabla_{i}^{{}} w_{ij} V_{j} } + \rho_{i}^{{}} \sum\limits_{j} {\frac{{c_{0}^{2} \rho_{ij} }}{{C_{ij} }}\frac{{{\varvec{r}}_{ij} }}{{\left| {{\varvec{r}}_{ij} } \right|}} \cdot \nabla_{i}^{{}} w_{ij} V_{j} } $$

(B.11)

Comparing Eqs. (B.1) and (B.11), the density diffusion by Riemann solution D^R can be written as:

$$ D_{i}^{{\text{R}}} = hc_{0}^{{}} \sum\limits_{j} {\left( {\frac{{\rho_{i}^{{}} c_{0}^{{}} }}{{C_{ij} }}\frac{{\left| {{\varvec{r}}_{ij} } \right|}}{h}} \right)\psi_{ij}^{{\text{R}}} \frac{{{\varvec{r}}_{ij} \cdot \nabla_{i}^{{}} w_{ij} }}{{\left| {{\varvec{r}}_{ij} } \right|^{2} }}V_{j} } $$

(B.12)

$$ \psi_{ij}^{{\text{R}}} = 2\rho_{ij} $$

(B.13)

Therefore, the density diffusion by Riemann solution has close similarity with the δ-SPH diffusion term without a first-order correction [58].

Linear momentum equation with the artificial viscosity term can be written as:

$$ \frac{{{\text{D}}{\varvec{u}}_{i}^{{}} }}{{{\text{D}}t}} = \sum\limits_{j} {m_{j}^{{}} \left( {\frac{{{{\varvec{\sigma}}}_{j}^{{}} + {{\varvec{\sigma}}}_{i}^{{}} }}{{\rho_{i}^{{}} \rho_{j}^{{}} }}} \right) \cdot \nabla_{i}^{{}} w_{ij} } + {{\varvec{\varPi}}}_{i}^{{{\text{AV}}}} $$

(B.14)

$$ {{\varvec{\varPi}}}_{i}^{{{\text{AV}}}} = \sum\limits_{j} {m_{j}^{{}} \Pi_{ij}^{{{\text{AV}}}} \nabla_{i}^{{}} w_{ij} } $$

(B.15)

$$ \Pi_{i}^{{{\text{AV}}}} = \left\{ {\begin{array}{*{20}l} {\frac{{\alpha^{{{\text{AV}}}} hc_{ij}^{{}} }}{{\rho_{ij}^{{}} }}\frac{{{\varvec{u}}_{ij} \cdot {\varvec{r}}_{ij} }}{{\left| {{\varvec{r}}_{ij} } \right|^{2} + \left( {0.1h} \right)^{2} }}} \hfill & {\quad \left( {{\varvec{u}}_{ij} \cdot {\varvec{r}}_{ij} < 0} \right)} \hfill \\ 0 \hfill & {\quad \left( {{\varvec{u}}_{ij} \cdot {\varvec{r}}_{ij} \ge 0} \right)} \hfill \\ \end{array} } \right.;\quad \rho_{ij}^{{}} = \frac{{\rho_{i}^{{}} + \rho_{j}^{{}} }}{2},\quad c_{ij}^{{}} = \frac{{c_{i}^{{}} + c_{j}^{{}} }}{2} $$

(B.16)

In the Riemann SPH, linear momentum continuity equation is written as (see Sect. 3.3):

$$ \frac{{{\text{D}}{\varvec{u}}_{i}^{{}} }}{{{\text{D}}t}} = \sum\limits_{j} {m_{j}^{{}} \left( {\frac{{{{\varvec{\sigma}}}_{j}^{{}} + {{\varvec{\sigma}}}_{i}^{{}} }}{{\rho_{i}^{{}} \rho_{j}^{{}} }}} \right) \cdot \nabla_{i}^{{}} w_{ij} } + {{\varvec{\varPi}}}_{i}^{{\text{R}}} $$

(B.17)

$$ {{\varvec{\Pi}}}_{i}^{{\text{R}}} = \sum\limits_{j} {m_{j}^{{}} \left( {\frac{{2C_{ij} \left( {u_{j} - u_{i} } \right)}}{{\rho_{i}^{{}} \rho_{j}^{{}} }}} \right)\nabla_{i}^{{}} w_{ij} } $$

(B.18)

$$ C_{ij} = \frac{{c_{i} \rho_{i} \sqrt {\rho_{i} } + c_{j} \rho_{j} \sqrt {\rho_{j} } }}{{\sqrt {\rho_{i} } + \sqrt {\rho_{j} } }} $$

(B.19)

$$ u_{i} = {\varvec{u}}_{i} \cdot \frac{{{\varvec{r}}_{ij} }}{{\left| {{\varvec{r}}_{ij} } \right|}};\quad u_{j} = {\varvec{u}}_{j} \cdot \frac{{{\varvec{r}}_{ij} }}{{\left| {{\varvec{r}}_{ij} } \right|}} $$

(B.20)

$$ \begin{aligned}u_{i} &= \left( {{\varvec{u}}_{i} + \frac{1}{2}\left\langle {\nabla {\varvec{u}}} \right\rangle_{i}^{c} \cdot {\varvec{r}}_{ij} } \right) \cdot \frac{{{\varvec{r}}_{ij} }}{{\left| {{\varvec{r}}_{ij} } \right|}};\quad\\ u_{j} &= \left( {{\varvec{u}}_{j} - \frac{1}{2}\left\langle {\nabla {\varvec{u}}} \right\rangle_{j}^{c} \cdot {\varvec{r}}_{ij} } \right) \cdot \frac{{{\varvec{r}}_{ij} }}{{\left| {{\varvec{r}}_{ij} } \right|}} \end{aligned}$$

(B.21)

where Eqs. (B.20) and (B.21) correspond to linear and second-order constructions of variables. Here, considering linear construction of variables (Eq. B.20), since we have

$$ u_{j} - u_{i} = {\varvec{u}}_{ij} \cdot \frac{{{\varvec{r}}_{ij} }}{{\left| {{\varvec{r}}_{ij} } \right|}} $$

thus, we can reformulate Eq. (B.18) as:

$$ {{\varvec{\Pi}}}_{i}^{{\text{R}}} = \sum\limits_{j} {m_{j}^{{}} \Pi_{ij}^{{\text{R}}} \nabla_{i}^{{}} w_{ij} } $$

(B.22)

$$ \Pi_{ij}^{{\text{R}}} = \left( {\frac{{2\left| {{\varvec{r}}_{ij} } \right|C_{ij} }}{{\rho_{i}^{{}} \rho_{j}^{{}} }}} \right)\left( {\frac{{{\varvec{u}}_{ij} \cdot {\varvec{r}}_{ij} }}{{\left| {{\varvec{r}}_{ij} } \right|^{2} }}} \right) $$

(B.23)

Therefore, one can see that the AV term and the Riemann term have similarity. Indeed, Meng et al. [22] discussed this similarity and proposed a limiter function for Riemann stabilisation term in order to avoid excessive dissipation. Note that a concise summary of δ versus Riemann density diffusive terms as well as artificial viscosity versus Riemann momentum diffusive terms is presented in Fig. 

Fig. 35

A concise summary of δ versus Riemann density diffusive terms as well as artificial viscosity versus Riemann momentum diffusive terms

35.