A novel semi-implicit scheme for elastodynamics and wave propagation in nearly and truly incompressible solids

This paper presents a novel semi-implicit scheme for elastodynamics and wave propagation problems in nearly and truly incompressible material models. The proposed methodology is based on the efficient computation of the Schur complement for the mixed displacement-pressure formulation using a lumped mass matrix for the displacement field. By treating the deviatoric stress explicitly and the pressure field implicitly, the critical time step is made to be limited by shear wave speed rather than the bulk wave speed. The convergence of the proposed scheme is demonstrated by computing error norms for the recently proposed LBB-stable BT2/BT1 element. Using the numerical examples modelled with nearly and truly incompressible Neo-Hookean and Ogden material models, it is demonstrated that the proposed semi-implicit scheme yields significant computational benefits over the fully explicit and the fully implicit schemes for finite strain elastodynamics simulations involving incompressible materials. Finally, the applicability of the proposed scheme for wave propagation problems in nearly and truly incompressible material models is illustrated.


Introduction
Incompressible material models are commonly encountered in science and engineering for modelling soils, clays, elastomers, rubbers, and biological soft tissues. In spite of the tremendous efforts towards the development of finite element schemes for computational solid mechanics, accurate and computationally efficient schemes for the simulation of elastodynamics under high-frequency loading conditions as well as wave propagation in incompressible materials are still scarce. With the increased need for the simulation of elastodynamics and wave propagation in biological soft tissues, for example, in elastography, see Sarvazyan et al. [48], Ophir et al. [46], Bercoff et al. [2], Ye el al. [55], Li and Cao [35], and references therein, the need for accurate and computationally efficient numerical schemes for incompressible material models has never been more important.
For problems governed by high strain rates or high-speed loading, for example, blast and impact loads, and wave propagation problems, explicit time integration schemes have proven to be computationally beneficial. However, the fundamental issue associated with the explicit schemes when used in combination with the nearly incompressible material models is the severe restriction on the critical time step posed by the bulk (or compressional) wave speed. As the Poisson's ratio, ν, approaches 0.5, the critical time step based on the bulk wave speed (c κ ) decreases by several orders of magnitude in comparison to that based on the shear wave speed (c μ ), as shown in Table 1.
One common approach for alleviating this problem for pure displacement formulations based on selective or reduced or selective-reduced integrations, see Malkus and Hughes [39] and Hughes [22], is selective mass scaling in which the density is selectively adjusted to reduce the high frequencies without detrimentally affecting C. Kadapa  Here, μ is the shear modulus and κ is the bulk modulus. Young's modulus, E = 1 and density, ρ 0 = 1 the low-frequency modes [45,53,56]. However, the main issue with selective mass scaling is the difficulty in computing the scaling parameter without compromising the accuracy and at the same time using lumpedmass matrices. Furthermore, such techniques are not applicable for truly incompressible material models for which hydrostatic pressure has to be computed as an additional solution variable in combination with the displacements.
The absence of the pressure term in the constraint equation for the truly incompressible case makes it impossible to develop a fully explicit numerical scheme for such models. For truly incompressible material models, fully implicit finite element schemes, for example, Rossi et al. [47], Scovazzi et al. [49,50], Liu and Marsden [36,37], and semi-implicit schemes based on fractional-step or projection methods which were originally developed for fluid problems, see Chorin [10,11], Temam [52], Kim and Moin [33], Guermond and Quartapelle [17], Guermond et al. [16], Lovrić et al. [38] are the only techniques available for modelling truly incompressible solids. In literature, Zienkiewicz et al. [59], Lahiri et al. [34], and Gil et al. [14] extended the fractional-step schemes to computational elastodynamics, and Caforio and Imperiale [7] extended a projection method to incompressible elastodynamics. While the implicit schemes become computationally expensive for large-scale models due to the need for inversion of large-scale sparse matrices during the iterative solution procedure at every time step, semi-implicit schemes offer several computational advantages over fully implicit schemes due to the substantial reduction in the size of matrix systems to be solved for.
Nevertheless, although fractional-step and projection schemes have been successfully adapted for elastodynamics problems modelled with incompressible materials, the fundamental problem associated with these schemes, especially in the context of solid mechanics problems, is the uncertainty over boundary conditions to be used during different steps of the projection schemes. Moreover, while the cost of computing correction terms for fractional-step schemes based on finite-difference based spatial discretisation is comparatively low, such schemes require numerical quadrature, which is often computationally expensive, for finite element schemes, especially for linear quadrilateral/hexahedral elements and higher-order triangular/tetrahedral finite elements. Additionally, contact modelling using fractional-step schemes is still an open problem. To the knowledge of the author, there is no published literature on modelling solid-solid contacts using fractional-step or projection schemes.
Another motivation behind the proposed work is the computational cost associated with the fully implicit scheme for elastodynamics using mixed formulations. In the literature on mixed formulations for computational solid mechanics, fully implicit schemes are well established, refer to Zienkiewicz and Taylor [61], Kadapa [28], Scovazzi et al. [47,49], and Janz et al. [24], and references therein. While these schemes do not pose severe restrictions on the time step size, the fact that large matrix systems need to be solved for at every iteration of each time step makes fully implicit schemes computationally very expensive, especially when the time step still has to be chosen small enough to obtain accurate results for problems which undergo extremely large strains which is often the case with hyperelastic material models. The proposed semi-implicit scheme also serves as a computationally efficient alternative to the fully implicit schemes, particularly for problems discretised with uniform meshes.
The proposed scheme is based on the mixed displacement-pressure formulation and exploits the lumped mass matrices for the displacement field in efficiently evaluating the Schur complement. This work is an extension of the unified finite element formulations for computational solid mechanics based on Bézier elements proposed in Kadapa [27,28]. In this scheme, deviatoric stress is treated explicitly and the pressure term implicitly, thus lifting the severe restrictions on the critical time step posed by the bulk wave speed. Moreover, the proposed scheme not only utilises the same finite element framework used for the implicit mixed formulation but also poses no ambiguity over the boundary conditions to be used. Additionally, the extension of the proposed scheme for modelling contacts using either Lagrange multipliers or penalty method or Nitsche method is straightforward, see [54,61] for details. The proposed formulation is equally applicable for mixed formulation with LBB-stable elements [28,29] as well as for equal-order elements with pressure stabilisation [8,9,40,42,43,47,49,57], as long as the mass matrix for the displacement field can be approximated as a lumped mass matrix. For the sake of simplicity, only the LBB-stable BT2/BT1 element proposed in Kadapa [28] is considered in the present work.
The outline of the paper is as follows. The mixed displacement-pressure formulation and the proposed semi-implicit scheme for problems in the small-strain regime are discussed in Sect. 2. The extension of the proposed semi-implicit scheme to finite strain regime is presented in Sect. 3. Efficient solution of the coupled matrix system using Schur complement is discussed in Sect. 4, followed by the discussion on the finite element discretisation in Sect. 5. The stability of the semi-implicit scheme is presented in Sect. 6. The performance of the proposed scheme and its applicability to elastodynamics and wave propagation in nearly and truly incompressible materials are assessed in Sect. 7 using several numerical examples. The paper is concluded with Sect. 8 by summarising the key features of the proposed work.

Governing equations
For linear isotropic elastic material, the Cauchy stress, σ , is related to the infinitesimal strain tensor, ε, as where μ is the shear modulus, κ is bulk modulus, ε dev is the deviatoric component of ε, I is the identity tensor of order two, and u is the displacement vector. The above relation (1) between the Cauchy stress tensor (σ ) and the small-strain tensor (ε) produces unphysical values for the incompressible case, i.e., for ν = 0.5 for which κ = ∞. Moreover, for the truly incompressible case, the finite element formulation based on displacements alone is no longer a viable choice. To overcome these issues, a mixed displacement-pressure formulation in which hydrostatic pressure (p) is computed as an additional solution variable, is employed c.f. [26,28,29,31,32].
The main idea behind the mixed displacement-pressure formulation is to split the stress into deviatoric and volumetric components and replace the volumetric part by an improved value. Accordingly, the modified Cauchy stress tensor is defined as where and p is the independent approximation for hydrostatic pressure. For nearly incompressible materials in small strains, p is related to the volumetric strain (ε v ) by the relation which becomes for the truly incompressible case, i.e., when κ = ∞. Therefore, for the truly incompressible case, pressure must be determined as an additional solution variable, together with displacement, for enforcing the volumetric constraint (5) on the displacement field. Using the above relations, the equations governing the elastodynamics of a solid body in the small strain regime can be written as where B 0 is the domain of the solid in the original configuration; X is an arbitrary point in the domain B 0 ; t is the time variable; T is the total time span; ∇ X is the original configuration gradient operator; ρ 0 is the density of the solid in the original configuration; v(= du dt ) is the velocity vector; a(= d 2 u dt 2 ) is the acceleration vector; u 0 is the initial displacement vector; v 0 is the initial velocity vector; n 0 is the unit outward normal on the boundary, ∂B 0 , of B 0 ; f 0 is the body force in the original configuration; u D is the prescribed displacement field on the Dirichlet boundary ∂B D 0 ; and t 0 is the prescribed traction force on the Neumann boundary ∂B N 0 . The Dirichlet and Neumann boundaries are such that ∂B 0 = ∂B D 0 ∪ ∂B N 0 and ∂B D 0 ∩ ∂B N 0 = ∅.

Finite element formulation
With w and q as the test functions for the displacement and pressure fields, respectively, the weak form for the mixed displacement-pressure formulation can be written as where dV and dA are the elemental volume and area, respectively. For the finite element analysis, the approximations for the displacement and pressure fields, u and p, and their test functions w and q are taken as w = N u w; u = N u u; q = N p q; and p = N p p, (8) where u and w are the vectors of displacement degrees of freedom (DOFs), and p and q are the vectors of the pressure degrees of freedom; N u and N p are the matrices of basis functions, respectively, for the displacement and pressure fields. By using the discretisation in Eq. (8), the semi-discrete equations for the mixed formulation (7) may be written as where M uu is the mass matrix and a is the vector of accelerations at the displacement DOFs. The internal force vector for the displacement DOFs, F int u , the external force vector for the displacement DOFs, F ext u , and the internal force vector for the pressure DOFs, F int p , are given by Here, D X and G X , respectively, are the divergence-displacement and gradient-displacement matrices with respect to the reference configuration, which, for the ath basis function of an element are given by Accordingly, the Cauchy stress tensor is represented as a column vector as The mass matrix, M uu , can be either consistent-mass matrix, M CM uu , or lumped-mass matrix, M LM uu . The consistent mass matrix is given as and the lumped-mass matrix, in this work, is computed using the row-sum lumping technique [60], according to which with n as the rank of the matrix M CM uu . The semi-discrete equations (9) can be solved in time by using either a fully implicit or fully explicit time integration scheme as presented in Kadapa [28]. While the fully explicit schemes are computationally appealing because of their matrix-free implementations, their performance is hindered by severe restrictions on the stable time step size due to bulk wave speed, especially when ν → 0.5, in addition to the fact that fully explicit schemes are not valid for the truly incompressible case. On the other hand, although the fully implicit scheme poses no such restrictions on the time step size, they become computationally expensive for finite strain problems due to the need for matrix solvers for the solution of large-scale sparse matrix systems during the iterative solution procedure at every time step, as will be demonstrated in the latter part of this paper. Therefore, we propose a novel semi-implicit scheme for which the time step depends only on the shear wave speed, and the size of the matrix system to be solved for is significantly lower than that of the fully implicit scheme.
For the sake of completeness and also to highlight the similarities and subtle differences between the semiimplicit scheme and the other two schemes, the fully implicit and fully explicit schemes are discussed briefly in Appendices A, B, and C. Note that, although the material in small strain regime is assumed to be linear elastic, the schemes are presented in the incremental form for the ease of extension to finite strain problems in Sect. 3.

Proposed semi-implicit scheme
For the semi-implicit scheme, the contribution to the internal force F int u in Eq. (9.1) due to pressure is treated implicitly while the contribution from the deviatoric stress is treated explicitly. Using the improved version of the explicit scheme of Chung and Lee [12] as presented in [27], the system of equations for the semi-implicit scheme for the mixed formulation (9) can be written as where t is the time step size, and α m , β, and γ are the time integration parameters. In Eq. (18) and in the subsequent parts of the paper, subscripts n and n + 1 refer to the previous and current time instants, t n and t n+1 , respectively. Using the spectral analysis [19], it can be shown that the explicit scheme is third-order accurate for The free parameter α m provides the user with control over the amount of numerical damping. Unless stated otherwise explicitly, α m = 1 in the numerical examples presented in this work. Now, with and by using the corresponding approximations from Eq. (8), the following matrix system results for the incremental displacement DOFs, u, and incremental pressure DOFs, p, at time step n + 1: where Equation (22) is in the exact block-matrix format of the mixed formulation using the fully implicit scheme (A.1). Therefore, the proposed semi-implicit scheme does not pose any additional requirements on the imposition of boundary conditions similar to the ones faced during various steps of projection-based and fractional-step schemes. Furthermore, with the proposed semi-implicit scheme, modelling of solid-solid contact using either penalty or Nitsche or Lagrange multiplier methods is straightforward, unlike fractional-step methods for which modelling of contacts is unclear.

Kinematics, stress measures, and governing equations
With B 0 and B t , respectively, as the original/reference and current configurations of the solid body under consideration, we can define the nonlinear deformation map that takes a point X ∈ B 0 to a point x ∈ B t as X : B 0 → B t . Now, by making use of the definition of the displacement, the deformation gradient (F) is defined as For modelling hyperelastic materials in the incompressible finite strain regime, the deformation gradient, F, is decomposed into deviatoric and volumetric components [3] as where Using the above definitions, modified versions of some important strain and stress measures are defined as Right Cauchy-Green deformation tensor, Green-Lagrange strain tensor, E := First Piola-Kirchhoff stress tensor, P : Cauchy stress tensor, σ : where J = det(F), and Ψ dev is the deviatoric part of the energy function of the material under consideration. The strain energy density functions for the material models considered in this work are assumed to be decomposed into the deviatoric part, Ψ dev , and the volumetric part, Ψ vol (J ), as For the truly incompressible material models, the total volume change is zero, which can be represented mathematically as The important consequence of the above relation is that the volumetric energy function (Ψ vol (J )) vanishes for the materials that are truly incompressible. For such material models, the pressure field must be determined as an additional solution variable that enforces the incompressibility constraint (38). For the comprehensive details on finite strain continuum mechanics, the reader is referred to Bonet and Wood [3], Zienkiewicz and Taylor [61], and Holzapfel [20]. For the mixed formulation in the finite strain regime, the effective first Piola-Kirchhoff stress is defined as where p is the independent approximation for hydrostatic pressure. From the above equation, the effective Cauchy stress tensor becomes Now, using the above definitions for strain and stress measures, the governing equations for the elastodynamics of incompressible material models in the finite strain regime can be written in the original configuration as where ρ 0 is the density of the solid in the original configuration B 0 ; ρ is the density of the solid in the current configuration B t ; and v(= du dt X ) and a(= d 2 u ) are the material time derivatives of displacement (u).

Finite element formulation
The linearised consistent mixed displacement-pressure formulation recently proposed by the author for elastostatic problems in Kadapa and Mokarram [31] is extended to elastodynamics problems in this work. Following [31], the semi-discrete equations for the mixed formulation in the finite strain regime can be written analogous to Eq. (9) as where M uu , a and F ext u are the same as those in the small-strain formulation discussed in Sect. 2.2, and F int u and F int p are given by The matrices G x and D x in Eq. (43) are the gradient-displacement and divergence-displacement matrices with respect to the current configuration. Similar to (14), G x for the ath basis function of an element is given as The quantities J and ϑ in Eq. (44) are evaluated as where J n is J evaluated from the displacement from the previously converged time step. Note here that J and ϑ are constants at the current load step. This approach results in a symmetric global matrix as demonstrated in Kadapa and Mokarram [31]. Since the fully implicit scheme for the solution of nonlinear equations (42) is not the primary focus of the proposed work, it is presented briefly in Appendix B for the sake of completeness. The semi-implicit scheme proposed in this paper is discussed in the following Section.

Proposed semi-implicit scheme
Similar to the semi-implicit scheme for the small strain problem as discussed in Sect. 2.3, the semi-implicit scheme for the mixed formulation in finite strain regime can be written as which, after linearising and then using the approximations in Eq. (8), can be written in the form of a coupled matrix system for the incremental displacement and pressure DOFs as where Remark 1 Note that the displacement and pressure DOFs at t n+1 for the semi-implicit scheme are computed in a non-iterative manner and that the contribution from the linearisation of configuration-dependent terms, for example, D x in Eq. (47.1) to the matrix K uu is ignored. The main reason behind this decision is to achieve the demonstrated computational benefits by exploiting the lumped approximation for the mass matrix in the efficient evaluation of the Schur complement, as discussed in Sect. 4. Moreover, due to the stability condition, the time step sizes for the semi-implicit scheme are often small. Since the fully implicit scheme requires no more than two or three iterations for such small time step sizes, ignoring the linearisation of the configuration terms has no significant effect on the accuracy of the results, as demonstrated with the numerical examples.

Efficient solution of the coupled system using the Schur complement
The coupled matrix systems in Eqs. (22) and (48) can be solved using either a direct solver or an iterative solver using block preconditioners as discussed in Benzi et al. [1]. However, such approaches for the solution of the coupled systems become computationally expensive for large-scale problems in three dimensions. To derive a computationally efficient scheme for the solution of the coupled matrix systems (22) and (48), in this work the approach of Schur complement reduction [1] is adapted. The novelty of the proposed work lies in the adaptation of the Schur complement, along with a lumped matrix, to solve the coupled matrix in a computationally efficient manner. Using the Schur complement reduction, the equations for the solution of incremental displacement and pressure in (22) and (48) can be written as where the Schur complement, S, is given as Although evaluation of the Schur complement is an expensive operation in general, even with a sparse K uu , it is worth highlighting at this point that it is not the case with the proposed semi-implicit scheme. The computational complexity of evaluating the Schur complement is drastically simplified in this work, thanks to the lumped-mass matrix M LM uu which renders K uu in Eqs. (22) and (48) as a diagonal matrix. Therefore, with a diagonal K uu , only sparse matrix multiplication and addition, which are significantly less expensive when compared to the inversion of sparse matrices, are the matrix arithmetic operations required for evaluating the Schur complement. Moreover, (i) since K −1 uu is a diagonal matrix, the Schur complement is still a sparse matrix, and (ii) since the size of the Schur complement is considerably lower than that of the coupled matrix systems (22) and (48), the computational cost of solving Eq. (54) is substantially lower than that of coupled systems, as demonstrated with numerical examples in Sect. 7.2.

Remark 2
The condition number of the coupled matrix systems (22) and (48), (A.1) and (B.1) deteriorates as the time step size is reduced. This ill-conditioning poses difficulties in adapting block preconditioners for the saddle point problem, see [1]. However, since the incremental DOFs in the semi-implicit scheme are solved using (54.1) and (54.2) in which the diagonal nature of the K uu is exploited to efficiently invert the Schur complement (S), and since the Schur complement is still symmetric, sparse and significantly reduced in size when compared with the original matrix system, the reduced system can be solved efficiently using either a sparse direct solver or Krylov iterative solvers with preconditioners, without worrying about the difficulties associated with solving a saddle-point system.

Finite element spaces for displacement and pressure
The proposed semi-implicit scheme is applicable to any combination of finite element spaces for the displacement and pressure fields for which the mass matrix for the displacement field can be approximated as a lumped mass matrix. However, such compatible displacement-pressure combinations are limited considerably by the inf-sup stability criterion for the saddle-point problems resulting from the mixed-formulations as well as the suitability of the displacement basis functions for lumped mass matrices. For the discussion on finite element spaces for the mixed formulations, the reader is referred to Brezzi and Falk [5], Brezzi and Fortin [6], and Brezzi, and Bathe [4].
The widely used inf-sup stable elements for incompressible solid mechanics and incompressible fluid flow problems are the Taylor-Hood family of elements [6,51]. However, because of the difficulties posed by the quadratic Lagrange elements for explicit schemes, the P2/P1 element is not a viable mixed element for the proposed semi-implicit scheme. Some suitable displacement-pressure combinations for the proposed semiimplicit scheme are: (i) the subdivision-stabilised NURBS proposed in Kadapa et al. [29] and (ii) the BT2/BT1 element proposed in Kadapa [28]. Because of the ease of mesh generation for complex geometries, see Kadapa [27], the BT2/BT1 element is considered in the present work. The application of the proposed scheme to any other suitable combination of displacement-pressure discretisation is akin to the BT2/BT1 element.
Sparse matrix patterns for the K pp matrix and the Schur complement matrix, S, for the BT2/BT1 element are presented in Fig. 1b, c, respectively, for the unstructured tetrahedral mesh shown in Fig. 1a. As illustrated, the Schur complement (55) resulting from the use of lumped-mass matrix, M LM uu , is not a full matrix; the Schur complement in the present work is a sparse matrix.

Stability of the semi-implicit scheme
For the fully explicit schemes, i.e., when both the deviatoric and the volumetric components in the internal force vector in Eq. (10) are treated explicitly, the stable time step is determined by the time taken by the bulk wave to travel through the smallest element in the finite element mesh [25,28,55]. However, since only the deviatoric term is treated explicitly in the proposed scheme, and since deviatoric stress is a function of shear modulus but not bulk modulus, the stable time step is limited by the propagation of the shear wave rather than the bulk wave. To prove this, Eqs.
Following Joly [25], we have where c μ (= √ μ/ρ 0 ) is the speed of the shear wave and h e is the element characteristic length. Now, using the condition (60), the sufficient condition for the stability of the proposed semi-implicit scheme is Therefore, the critical time step for the proposed semi-implicit scheme depends only on the shear modulus but not the bulk modulus.

Numerical examples
The time step sizes for the semi-implicit and the full-explicit schemes are computed as where CFL is the Courant-Friedrichs-Lewy number which, in this work, also takes into account the parameters All the simulations presented in this work are performed using the in-house C++ code on a single Intel i7-8750H CPU on a personal computer. For all the examples, the spectral radius for the fully implicit scheme is taken as ρ ∞ = 0. The free parameter for the semi-implicit scheme is α m = 1, unless stated otherwise explicitly. The mass matrix used for the fully explicit and the semi-implicit schemes is the lumped-mass matrix while the mass matrix for the fully implicit scheme is the consistent-mass matrix, if not specified otherwise.

Material models for finite strain problems
The proposed semi-implicit scheme for finite strain problems is applicable to any nearly or truly hyperelastic material models, for example, Neo-Hookean, Mooney-Rivlin, Veronda-Westmann, Ishihara and Ogden, etc., that are used for modelling polymeric materials and biological soft tissues [18,21,41,44]. However, for the purpose of demonstration, the present work is limited to Neo-Hookean and Ogden material models. The deviatoric part of the energy function, and the corresponding material properties used in the present work for the Neo-Hookean and Ogden models are given as: • Neo-Hookean model: where I C = tr(C). Young's modulus for this material model is assumed to be E = 1.2 × 10 6 Pa.
• Ogden model: Here, λ 1 , λ 2 , and λ 3 are the principal stretches. For the Ogden material, the shear modulus results from The material parameters for this material model are taken as The volumetric energy function for the nearly incompressible case is assumed as

Convergence studies
The ability of the proposed scheme in computing accurate numerical results is demonstrated by evaluating the error norms in displacement, pressure, and stress using the example of a manufactured solution. The problem domain consists of unit square/cube for 2D/3D problems with Dirichlet boundary conditions specified on all the boundary conditions. The Young's modulus and density of the material are assumed to be E = 100 and ρ 0 = 1. The analytical solution used for the 2D problem is and for the 3D problem is The corresponding body force, f 0 , is computed by evaluating the modified Cauchy stress (2) first and then substituting it in the governing equation (6.1). Sample finite element meshes used for this example are presented in Fig. 2. Results obtained with the proposed semi-implicit scheme are compared against the results obtained with the fully implicit scheme. L 2 norms of error in displacement, pressure, and stress at t = 1.5 for the 2D and 3D problems with two different values of Poisson's ratio, ν = 0.499 and ν = 0.5, are presented in Fig. 3, along with the corresponding values obtained with the fully implicit scheme. As illustrated, the solution obtained with the proposed semiimplicit scheme converges with optimal convergence rates in displacement, pressure, and stress fields.

Twisting of a column: nearly incompressible case
In this example, the computational benefits of the proposed semi-implicit scheme over a fully explicit scheme are assessed by studying the problem of twisting of a column modelled with the nearly incompressible Neo-Hookean material model. The column is of size 2 cm × 12 cm × 2 cm and is clamped at its bottom end, and it is excited with an initial velocity, v = 1500 sin(π y/12) (z, 0, −x) cm/s. This example is studied with two meshes M1 and M2 shown in Fig. 4. The density is taken as ρ 0 = 1.1 g/cm 3 and simulations are performed for different values of Poisson's ratio using the fully explicit and semi-implicit schemes. The performance is assessed by studying the evolution of Y-displacement of point A and the amount of time required for each simulation to reach 10 ms. The Y-displacement of point A obtained with the proposed semiimplicit scheme is in good agreement with the corresponding values obtained with the fully explicit scheme, as presented in Fig. 5. The computational time required for the simulation to complete 10 ms for each mesh using both the schemes and different values of Poisson's ratio is presented in Table 2. As shown, the computational cost of the semi-implicit scheme remains almost constant irrespective of the value of Poisson's ratio while the computational cost of the fully explicit scheme increases with increasing value of Poisson's ratio; the cost of the fully explicit scheme exceeds that of the semi-implicit scheme for ν ≥ 0.45, and for ν = 0.49999, the fully explicit scheme is about 100 times computationally more expensive than the semi-implicit scheme. Even though the proposed semi-implicit scheme requires matrix-vector, matrix-matrix multiplications, and sparse matrix solvers, significant computational benefits can be gained by using the proposed semi-implicit scheme over the matrix-free fully explicit scheme.

Twisting of a column: fully incompressible case
In this example, the performance of the proposed semi-implicit scheme in comparison to the fully implicit scheme is assessed by studying the column from the previous example using the truly incompressible Neo-Hookean and Ogden material models. The evolution of Y-displacement of point A using the fully implicit and the proposed semi-implicit scheme for the neo-Hookean material model shown in Fig. 6 proves the convergence of the proposed semi-implicit scheme as the mesh is refined. The difference in the numerical results between the fully implicit and semiimplicit schemes obtained with the coarse mesh M1 is due to the use of consistent mass matrices in the fully implicit scheme and the lumped mass matrix in the semi-implicit scheme, as demonstrated in the same Figure with the solution obtained with the fully implicit scheme using the lumped mass matrix. Contour plots of The critical time step size for the semi-implicit scheme is of the order of 10 −2 ms irrespective of the value of Poisson's ratio while that for the fully explicit scheme decreases from the order 10 −2 ms to 10 −4 ms as Poisson's ratio is increased from 0.3 to 0.49999. The matrix solver used for the SI scheme is the SuperLU sparse direct solver called from the opensource C++ library Eigen the pressure field presented in Fig. 7 for the M2 mesh using the neo-Hookean material model demonstrate excellent agreement between the pressure fields obtained with the proposed semi-implicit scheme and the fully implicit scheme. The results obtained with the Ogden material model, as shown in Figs. 8 and 9, illustrate the applicability of the proposed semi-implicit scheme to generic hyperelastic material models.

Computational efficiency
The amount of time taken for each simulation to reach 10 ms using the semi-implicit scheme and the fully implicit scheme with different time steps for both the meshes is presented in Table. 3, and the Y-displacement plots obtained with different time steps with M2 mesh are presented in Fig. 10. From these results, it is evident that the semi-implicit scheme outperforms the fully implicit scheme in terms of computational efficiency even though the time step for the semi-implicit scheme is limited by the shear wave speed and the smallest element size in the mesh. The computational cost of the fully implicit scheme for computing the numerical results that are of comparable accuracy to that of the semi-implicit scheme is an order of magnitude higher than that of the semi-implicit scheme. From the results shown in Table. 3, it is also apparent that the computational cost of the fully implicit scheme would be substantially higher for the same time step as that of the semi-implicit scheme. While this might not be the case for problems consisting of non-uniform meshes in which element sizes vary by several orders, the proposed semi-implicit scheme can still serve as an efficient alternative to the fully implicit scheme for problems with almost uniform meshes. The reason behind the computational efficiency of the semi-implicit

Stent model
This example is concerned with the application of the proposed scheme for modelling problems with complex geometries. For this purpose, a stent model with the geometry as shown in Fig. 11 is considered. Due to the symmetry, only 1/8th of the domain is considered for the analysis. The finite element mesh shown in Fig. 11 consists of 6815 nodes and 3113 elements. The material is assumed to be truly incompressible Ogden model. The symmetry boundary condition is applied on all the outer faces parallel to the coordinate axes, and uniform pressure of 5000 Pa is applied on the inner surface. The time step size is of order of 10 −3 ms. The values of the displacement and velocity of a corner point displayed in Fig. 12 show that the results obtained with the (a) (b) Fig. 9 Twisting column: contour plots of pressure field at a t = 2 ms, and b t = 6 ms obtained with the truly incompressible Ogden material model  The value in the brackets is the ratio of the time take by the FI scheme over the time taken by the SI scheme for the corresponding mesh. The stable uniform time step sizes used for the M1 and M2 meshes for the SI scheme are 0.04 ms and 0.02 ms, respectively proposed semi-implicit scheme match well with those obtained with the fully implicit scheme. Contour plots of pressure overlaid on the deformed configurations of the stent at two different time instants, as shown in Figs. 13 and 14, illustrate smooth pressure fields obtained using the proposed semi-implicit scheme. Thus, this example showcases the applicability of the proposed semi-implicit scheme for performing elastodynamics simulations of truly incompressible material models using unstructured tetrahedral meshes that can be readily generated using the existing mesh generators.

Circular shear wave propagation
This example has been previously studied in Ye et al. [56]. The geometry of the problem consists of 60 × 60 mm 2 square block, and it is studied using the plane-strain assumption and the small-strain formulation. Two meshes, M1 and M2, with 60 × 60 × 4 and 120 × 120 × 4 BT2/BT1 elements, respectively, are considered for this example. The M1 mesh and two points at which the response is evaluated are shown in Fig. 15a. Homogeneous Dirichlet boundary conditions are applied on all the boundaries, and traction forces defined by the following Ricker wavelet function are applied at the centre of the domain, as shown in Fig. 15b: where the frequency f = 300 Hz, t 0 = 0.003 s, and the amplitude is A = 0.0002 N for mesh M1 and A = 0.0004 N for mesh M2. The shear modulus is μ = 0.001 MPa and the density is ρ 0 = 1000 kg/m 3 . Two different values of Poisson's ratio, ν = 0.49995 and ν = 0.5, are considered in order to illustrate the applicability of the proposed semi-implicit scheme for simulating wave propagation in a nearly as well as a truly incompressible elastic continuum.
Nearly incompressible case (ν = 0.49995) For this case, results obtained with the semi-implicit and fully explicit schemes match well, as shown with the evolution of Y displacement at points A and B in Figs. 16 and 17, and the contour plots of shear stress in Fig. 19. The minor differences between the present results and those published in Ye et al. [56] are attributed to the differences in the formulations, spatial discretisations, and the time integration schemes used in the present work and Ye et al. [56].  Fig. 19a, b. For this case, the fully explicit scheme takes 153 and 3798 s for M1 and M2 meshes, respectively, while the semi-implicit scheme requires only 51 and 628 s. The computational gains of the proposed semi-implicit scheme would be even more significant for 3D problems.
Truly incompressible case (ν = 0.5) For the truly incompressible case, the fully explicit scheme is not valid; therefore, analysis is carried out only with the semi-implicit scheme. The evolution of Y displacement at points A and B obtained with M2 mesh is presented in Fig. 18 along with the corresponding values obtained with ν = 0.49995. As expected, the difference between the results obtained with nearly and truly incompressible models is negligible. This case proves the applicability of the proposed semi-implicit scheme for computing the numerical solutions for wave propagation in truly incompressible material models without the need for any ad-hoc penalisation terms in approximating or imposing the incompressibility constraint (Fig. 19). This example is concerned with wave propagation in a truly incompressible elastic continuum with voids of different sizes and shapes, as depicted in Fig. 20a. Plane-strain condition is assumed, and the material is considered as linear elastic with the material properties E = 10 GPa and ρ 0 = 2500 kg/m 3 . The domain is fixed on the left, bottom and right edges, and its top edge is traction free. The wave is initiated by applying a constant force (F) of magnitude 10 5 N in the negative Y-direction at the centre of the domain (point O in Fig. 20a) for a period of 0.2 s.
The finite element mesh used for the analysis is shown in Fig. 20b. The mesh consists of 31533 nodes and 15496 elements. To minimise oscillations due to high-frequency modes, the parameter α m is chosen as α m = 1.2. The displacement field and shear stress at four different time instants during the propagation of the wave are illustrated, respectively, in Figs. 21 and 22. These contour plots show that the solution obtained with the proposed scheme is free from any spurious oscillations. Thus, this example showcases the suitability of the proposed semi-implicit scheme for simulating wave propagation in problems with complex geometries using unstructured meshes that can be readily generated using existing mesh generators.  A novel semi-implicit scheme is proposed for the simulation of elastodynamics and wave propagation problems modelled with nearly and truly incompressible material models. The proposed scheme is based on the idea of efficient evaluation of the Schur complement by using a lumped mass matrix for the displacement field. Since the pressure field is treated implicitly, the critical time step is controlled by the shear wave speed rather than the bulk wave speed. The optimal convergence of the proposed scheme is evidenced by using the LBB-stable BT2/BT1 element. The results obtained for the example of twisting of a column prove that the proposed semi-implicit scheme yields significant computational gains when compared with both the fully explicit and the fully implicit schemes in computing numerical solutions of elastodynamics problems modelled with incompressible material models.   Fig. 22 Wave propagation in a complex elastic continuum: shear stress at four different time instants. The warp factor is 50 The applicability of the proposed scheme in performing elastodynamics simulations of problems with complex geometries is illustrated using the example of a stent.
With the examples of shear wave propagation in a simple and a complex domain with voids, the suitability of the proposed scheme in simulating wave propagation in nearly and truly incompressible elastic media is demonstrated.
In conclusion, the computational benefits of the proposed semi-implicit scheme make it a computationally efficient alternative to the fully implicit and fully explicit schemes in computing accurate numerical solutions of elastodynamics and wave propagation problems modelled with nearly and truly incompressible material models. The present work enhances the capabilities of the finite element framework based on Bézier elements [27,28] in computing computationally efficient solutions of elastodynamics and wave propagation in incompressible materials.
The proposed work offers numerous possibilities for its extension to applications in soft as well as smart materials because of their inherent incompressible nature. A straightforward extension of the proposed semiimplicit scheme would be for the simulation of wave propagation in soft periodic structures for applications in metamaterials [13,23,58].