Introduction

The test case of a manufactured torsional deformation is a relevant benchmark test in solid mechanics simulations including non-linear and finite strain dynamics for rigorous validation of accuracy and convergence. This test was originally proposed by Brannon et al. (2011) targeting Material Point Method (MPM), and was employed in the context of SPH by several researchers including Tang et al. (2024) and Gotoh et al. (2025). In the paper by Gotoh et al. (2025) the comprehensive derivations of theoretical solutions were presented including those for pressure and total energy. In this test, through consideration of the momentum equation and the considered non-linear material constitutive model, external body forces are derived to result in a prescribed or desired displacement field corresponding to a torsional deformation. Thus, the accuracy and convergence of the structure model can be verified in terms of displacement, stress and energy, in a scrupulous manner.

Physical phenomena

• solid mechanics

• non-linearity

• finite strain

Geometry

The geometry is a 2D circular annulus with an outer radius R₀ = 1.25 m and an inner radius R_i = 0.75 m, as shown below in Figure 1. The physical properties of the annulus are set as ϱ₀ = 1.0E+3 kg/m³, E = 1.0E+3 Pa and ν = 0.3 (Gotoh et al. 2025).

Boundary conditions

The traction free (surface stress free) boundary condition needs to be imposed on the surface of the annulus. This is achieved as a fixed boundary condition, i.e., by fixing the particles in the outermost and innermost layers of the annulus, or at least restricting their motion only in the circumferential direction.

Initial conditions

The initial set up is described in Fig. 1. The circular annulus is initially stationary. After the initial state, it experiences a torsional deformation of angle α (dependent on the radial distance, R, and time, t, caused by a time-varying external body force. In this test, the desired deformation field is specified by the following deformation angle α(R,t) = g(t) h(R) as follows:

g(t) = A/2 {1 - cos(2πt/T)} ; h(R) = (4R - 3)² (4R - 5)² (1)

The coefficients in the above formula are given as A = 0.3 rad and T = 1.0 s. In other words, the annulus is initially undeformed and reaches its maximum deformation at T/2 = 0.5 s, returning to its undeformed state at T = 1.0 s. This cycle is repeated with a period of T = 1.0 s (Gotoh et al. 2025).

The external body force b (with the dimension of acceleration) required to produce the specified deformation (by deformation angle α(R,t)) is determined as follows in polar coordinates:

b_r = - R{g'(t) h(R)}² + μ/ρ₀ R {g(t) h'(R)}² (2)

bθ = R g''(t) h(R) - μ/ρ₀ {3h'(R) + R h''(R)} g(t) (3)

where μ is one of Lame’s constants corresponding to the shear elastic modulus, g'(t) and g''(t) indicate the derivative and 2nd derivative of g with respect to the quantity t in parentheses, respectively.

Discretization

The initial distribution of the particles

The initial distribution of the particles is set in two ways. First, the particles are arranged in a regular and evenly spaced manner (referred to as the “unperturbed case”). Specifically, the particles are evenly placed on concentric circles with radius increment equal to the particle spacing. Second, the initial positions of the particles are randomly shifted from the unperturbed state by a maximum of a certain small value (0.1% in Gotoh et al. 2025) of the particle spacing d₀ (referred to as the “perturbed case”). Note that the particles in the outermost and innermost layers are shifted only in the circumferential direction since they define the edge of the annulus.

Specifically, the function for the perturbation of particle positions is given as follows:

(2RAN - 1) e_perturb d₀ (4)

Where RAN: a random number between 0 and 1, e_perturb: the coefficient setting the magnitude of perturbation, which is set 0.001 in Gotoh et al. (2025), d₀: particle spacing.

Error functions for quantitative validations

In order to quantitatively evaluate the solid model, the error functions (relative error) of displacement, pressure and total energy are defined as follows:

Error_Disp = (∑_iᴺ |(d_theo )_i - (d_num )_i| ) / (∑_iᴺ |(d_theo)_i|) (5)

Error_Pres = (∑_iᴺ |(p_theo)_i - (p_num )_i| ) / (∑_iᴺ |(p_theo)_i|) (6)

Error_T.E. = |T.E._theo - T.E._num| / |T.E._theo| (7)

where d is the displacement, p is the pressure, T.E. is the total energy of the entire system, N is the total number of particles, and i corresponds to each target particle. The subscripts “theo” and “num” indicate theoretical values and numerical results, respectively. The theoretical values and their derivation are described comprehensively in the next section “Theoretical solutions”. Note that the subscript “theo” is omitted in the next section “Theoretical solutions”. These quantitative error graphs are recommended to be presented along with some qualitative results such as pressure field.

Theoretical Solutions

Figure 2 illustrates the definition of position, angle, coordinates and bases (represented by vector e). The initial position vector X (magnitude R, azimuth Θ) is rotated by a torsional deformation of angle α to become vector r (magnitude R, azimuth Θ). The following relationship holds for angles: θ = Θ + α.

Theoretical solutions corresponding to the manufactured torsional deformation can be found here. These theoretical solutions comprise external body forces, pressure and total energy of the structural system corresponding to the prescribed deformation field.

Results specification

Presentations and comparisons of the simulation results (both “unperturbed case” and “perturbed case” described in the beginning of the section “Discretization”) with the theoretical values presented in the section “Theoretical solutions". Specific targets for comparison include:

• snapshot of pressure field

• convergence analysis with the error functions (relative error) of displacement, pressure and total energy defined by Eqs. (5)-(7) described in the section “Discretization”

Typical result of comparison is presented in the section “Benchmark results.”

Benchmark results

In this section, the examples of numerical results for this benchmark test are presented (Gotoh et al., 2025). First, the results of the unperturbed case are presented by a representative qualitative figure corresponding to reproduced and theoretical pressure field (Fig. 3) and a quantitative one related to error functions of displacement, pressure and energy (Fig. 4). Second, the corresponding information is presented in Figs. 5 and 6 for the case of initially perturbed particle distribution.

Fig. 3 presents the snapshots of particles along with pressure field for TLSPH-C1st (C1st is presented by Khayyer et al. 2021), TLSPH-C2nd, and TLSPH-C2nd-cR, together with the corresponding theoretical values at t = 3.74 s for the case of unperturbed particle distributions with an initial particle spacing of d₀ = 0.025 m. From the presented figure, the snapshot by TLSPH-C1st is characterized by a clearly disordered particle distribution at this instant due to presence of hourglass modes. On the other hand, through the incorporation of second-order consistent approximations of deformation gradient and corresponding acceleration, the TLSPH-C2nd presents a stable state of particle distribution together with a stable pressure field. For this initially unperturbed particle distribution case, the result by TLSPH-C2nd-cR is qualitatively identical to that by TLSPH-C2nd, both showing close similarity with the corresponding theoretical solution. Here, the pressure value is negative. As can be seen from Eq. (29), this indicates that the surrounding normal stresses are, on average, tensile stresses.

Fig. 4 illustrates the results of the three error functions defined by Eqs. (5)-(7) in order to confirm the accuracy and convergence of TLSPH-C1st, TLSPH-C2nd, and TLSPH-C2nd-cR at two time instants, t = 0.50 s (the instant of the first maximum displacement) and t = 3.74 s (the instant of the presented snapshot in Fig. 3), with the initial particle spacings of d₀ = 0.05 m, 0.025 m, 0.0125 m, 0.00625 m, and 0.003125 m for the case of unperturbed particle distributions. This figure portrays that through application of C2nd and cR, second-order convergence is ensured for the displacement field at both considered instants. As expected, the order of convergence for pressure and total energy becomes almost one order lower than that of displacement because these two physical quantities depend on the deformation gradient which is achieved by taking the gradient of displacement field.

Fig. 5 presents the snapshots of particles together with the pressure field at t = 3.74 s for the case of initially perturbed particle distributions and an initial particle spacing of d₀ = 0.025 m. Similar to the unperturbed case, the snapshot by TLSPH-C1st is characterized by hourglass instability and in a more intensified state due to the deliberate initial perturbation in particle positions. Through incorporation of C2nd scheme, a more stable particle distribution is achieved despite presence of noise in the reproduced pressure field. Through application of the cR scheme, the pressure field is stabilized and the results by TLSPH-C2nd-cR show a close qualitative similarity with the corresponding theoretical solution.

Fig. 6 depicts the graphs of the three error functions for accuracy and convergence confirmation for the initial particle spacings d₀ = 0.05 m, 0.025 m, 0.0125 m, 0.00625 m, 0.003125 m at t = 0.50 s and t = 3.74 s for the case of initially perturbed particle distributions. From the presented figure, even in presence of initial perturbation in particle positions, the TLSPH-C2nd-cR has provided accurate predictions characterized by second-order convergence. Due to the imposed perturbation at the beginning of calculations, for the considered fine particle spacings, the order of convergence corresponding to pressure is lower than that of the unperturbed case.

Download

You can download the full test case below:

References

• Gotoh, T., Sakoda, D., Khayyer, A., Lee, C.H., Gil, A.J., Gotoh, H., Bonet, J. (2025). An enhanced total Lagrangian SPH for non-linear and finite strain elastic structural dynamics. Comput Mech. https://link.springer.com/article/10.1007/s00466-024-02592-z

• Tang, X., Wu, D., Wang, Z., Haidn, O., Hu, X. (2024). An explicit multi-time stepping algorithm for multi-time scale coupling problems in SPH. Commun Comput Phys 36(5): 1219–1261. https://doi.org/10.4208/cicp.OA-2024-0037

• Brannon, R.M., Kamojjala, K., Sadeghirad, A. (2011). Establishing credibility of particle methods through verification testing. In: II international conference on particle-based methods—fundamentals and applications, (PARTICLES 2011). Barcelona. https://upcommons.upc.edu/handle/2117/189792 

• Kamojjala, K., Brannon, R., Sadeghirad, A., Guilkey, J. (2015). Verification tests in solid mechanics. Eng Comput 31:193–213. https://doi.org/10.1007/s00366-013-0342-x 

• Nishawala, V.V., Ruggirello, K.P. (2015). Evaluating the material point method in CTH using the method of manufactured solutions. Sandia National Lab. (SNL-NM), Albuquerque, New Mexico (United States). https://www.osti.gov/biblio/1458891 

• de Vaucorbeil, A., Nguyen, V.P., Hutchinson, C.R. (2020). A Total-Lagrangian Material Point Method for solid mechanics problems involving large deformations. Comput Methods Appl Mech Eng 360:112783. https://doi.org/10.1016/j.cma.2019.112783 

• Khayyer, A., Shimizu, Y., Gotoh, H., Nagashima, K. (2021). A coupled incompressible SPH-Hamiltonian SPH solver for hydroelastic FSI corresponding to composite structures. Appl Math Model 94:242-271. https://doi.org/10.1016/j.apm.2021.01.011