The influence of fracture energy on wooden structural members due to contact explosion

Because wood is one of the oldest construction materials, there are plentiful references to the mechanical properties of different types of wood. Not only its elastic properties have been described. For instance, elastic properties of mass timber panels were determined by modal testing used on full-size solid wood (Zhou et al. 2017). Moreover, common yew shear moduli were determined with ultrasonic waves (Keunecke et al. 2007). Wood density as well as its change during radial growth during the ripening period of oak and black pine was studied in Krajnc et al. (2021). Even the opposite, that is, the effect of reducing the density of wood while increasing the number of trees in forests, has been studied in Pretzsch et al. (2018). The fracture energy was also studied, for instance for spruce in Stanzl-Tschegg et al. (1994) or for wood cracks considering moisture levels (Vasic and Stanzl-Tschegg 2007). Noteworthy also are studies of wood subjected to a tensile force in the longitudinal direction (Navi et al. 1995), detailed descriptions of wood defects due to its structure, such as the influence of knots on the properties of wood (Wood Handbook 2010), 3D modelling of knots (Lukacevic et al. 2019), and descriptions and processes that affect the formation of cracks and deformations in statically loaded timber (van Blokland et al. 2020).

Furthermore, information was also obtained on the fracture behaviour of birch and spruce in the radial–tangential direction of propagation (Tukiainen and Hughes 2013). Moreover, a vacuum drying system for green hardwood parts (Chen and Lamb 2004) as well as acoustic wave propagation through a cross-section of green wood (Dikrallah et al. 2010) were studied. The methodology of modelling multi-layer structures (Gliniorz et al. 2002) and glued–laminated timber designs (Stürzenbecher et al. 2010), such as the mechanical properties of innovative laminated beams made of bambus wood (Fang et al. 2015) and analyses and measurements of wood-based panels (Olmati et al. 2014; Zhou et al. 2016), have also been addressed in the literature.

The present research study specifies, among other things, the elastic properties of birch, oak and pine. Table 1 shows the elastic properties of selected wood species acquired from the literature. Note large discrepancies between the values for different wood species, but also for the same genus. For this reason, the orthotropic moduli of elasticity were determined in this study.

Various aspects of the fracture mechanics of wood have been studied in the literature. Novel measurement techniques were used for the observations, such as micro-scale computed tomography (Forsberg et al. 2008), because the fracture phenomenon is induced by the wood microstructure (Carlsson and Isaksson 2020). Moreover, the fracture energy of glued joints may be determined, as in Ammann and Niemz (2015), who investigated the influence of humidity on the crack propagation rate for European beech wood. The study showed no relationship, and all tests indicated similar values. Further, the cracking of wood and its composites is influenced by many parameters; Stanzl-Tschegg and Navi (2009) further investigated the mechanisms involved during cracking. In Ostapska and Malo (2020), the fracture parameters of Norway spruce were measured and described using a digital image correlation system. According to the authors` knowledge, no study has examined the fracture energy of birch, oak and pine woods or its influence on the blast properties of those wood species.

In this study, compression tests of birch (Betula verrucosa Ehrh), oak (Quercus robur L.) and pine (Pinus sylvestris L.) wood were conducted to determine the moduli of each species of wood in three directions. Green wood samples were considered, because the sappers work in field conditions, where the structures to breach are built from green wood material (from felling, built ad hoc), i.e., from locally available trees rather than seasoned wood typical for engineering structures. The samples were tested in three main directions of wood: radial, tangential and longitudinal. Polish norms PN-D-04102 (1979) and PN-D-04229 (1997) were used to obtain the mechanical properties of those wood species. The first norm regulates the experimental procedure for compression wood testing along the grain (longitudinal direction), while the second norm standardises compression wood testing across the grain (radial and tangential directions).

Cuboidal samples were used in all tests, regardless of the direction. The samples for radial and tangential direction tests were \(20\times 30\times 30\) mm, whereas the longitudinal direction tests used samples of \(20\times 20\times 30\) mm. First, the slices from wood logs were cut by a chainsaw and a circular table saw; see Fig. 1a. The slices were selected to avoid or limit wood defects, which would affect the experimental results (the samples closer than 5 cm from the defect boundary, considered in the slice plane, were eliminated). The slices were cut further by a circular table saw to obtain the wood strips; see Fig. 1b–c. Finally, the circular table saw was used to obtain cuboidal samples; see Fig. 1d. Three direction samples were cut from one strip: radial, tangential and longitudinal. They were marked as R, T and L in Fig. 1e. Some samples were excluded from the tests due to defects or damage sustained during their acquisition.

All tests were conducted on a typical compression press. In the present study, the INSTRON 5982 press was used with a force capacity of \(100\) kN; see Fig. 2a. The speed of the displacement-controlled test was \(1.5\) mm/min (mean data acquisition sampling was approximately 16 Hz). In the tests, time, force and displacement were each recorded, as was the strain in an axial direction by the video extensometer; see Fig. 2b. The final number of measured force–displacement test curves for birch was 28, 28 and 28; for oak it was 45, 36 and 54; and for pine it was 20, 24 and 24 for the radial, tangential and longitudinal directions, respectively.

Post-processing computations were performed for each of the 287 samples tested. First, the test output data were smoothed to compute the elastic moduli. The smoothing technique used was an author's “marching” algorithm, in which 20 points (10 before and 10 after the point of interest) were used to average a single point on the smoothed curve. Second, each force vs displacement curve was divided into 400 points and interpolated for those displacement values. Next, the stress vs strain curves were computed from those data. The basic mechanical approach was used here, in which \(L\) was the sample height. As mentioned above, the dimensions for the samples were different depending on the wood direction considered. Due to different in-plane dimensions, the area of the wood samples for the longitudinal direction was \(400\) mm². For the radial and tangential directions, the area was \(600\) mm². The height of the samples was the same for all three directions: \(L=30\) mm. Then, due to the nonlinear character of the plot, multiple linear regressions were computed for the marching collection of m points (\(m=20\)). From each linear regression, the elastic moduli were determined by receiving the slope ratio for the data considered. The maximal value of the moduli obtained (for the almost linear part of the stress vs strain plot) was interpreted as the sought elastic moduli of the sample in the particular wood direction.

Wooden members are destroyed with blasting agents of normal strength. Charges can be applied to logs, beams, and groups or bundles of piles, both externally and internally. For calculations, several virtual diameters of round wooden beams were assumed: 15 cm, 20 cm, 25 cm, 30 cm and 35 cm. Those diameters were assumed for three different wood species, considering the humidity. The destructive explosive loads may be computed from the following formulas:where \(W\) is the mass of a charge, \(K\) is a tabular coefficient (see Table 2) representing wood type and humidity at the same time, and \(D\) is the wood diameter in cm (measured in the plane of a puncture).

A realistic approach to modelling the constitutive behaviour of wood is required to determine the numerical estimation of the destructive explosive loads of wooden structural elements. This approach has been adopted based on the understanding that static wood may be efficiently modelled by orthotropic elasticity. It directly results from the structure of wood; see Fig. 1c. The wood is stiffer along the grain. Both the radial and tangential properties of wood perpendicular to the grain are different, and wood is weaker in both directions. Thus, in this paper, an elastic orthotropy was used in numerical modelling, (Dudziak et al. 2016).

The constitutive law for wood was determined by specifying the terms in the elastic stiffness matrix, \({D}_{ijkl}\). The elastic stiffness matrix \({D}_{ijkl}\) defines the relationship between stresses \({\sigma }_{ij}\) and strains \({\varepsilon }_{kl}\):

The elastic stiffness matrix \({D}_{ijkl}\) takes the following form: where

The cylindrical system of coordinates was used in the present case, in which Cartesian axes for classic orthotropy—1, 2 and 3—were replaced with R: radial, T: tangential and L: longitudinal directions, respectively.

The unknown elastic parameters of the woods (birch, oak and pine) used in this study were in part determined experimentally, and others were assumed according to the literature. The experimental determination of elastic moduli \({E}_{R}\), \({E}_{T}\) and \({E}_{L}\) is presented in the Results section. Poisson's ratios and shear modulus are required for orthotropic elasticity, apart from elastic moduli. Poisson’s ratios were assumed according to Wilczyński and Gogolin (1989), and the mean values were computed from the data reported in the study: \({\nu }_{RT}=0.683\), \({\nu }_{RL}=0.054\) and \({\nu }_{TL}=0.033\). Other Poisson’s ratios, namely, \({\nu }_{TR}\), \({\nu }_{LR}\) and \({\nu }_{LT}\), were computed according to the symmetry:

where \(\{i,j\}\) are \(R\), \(T\) or \(L\). Shear moduli were computed according to the classical formula of orthotropic elasticity, namely:

In the paper, apart from elastic orthotropy, the fracture behaviour of wood was modelled by applying the element deletion through strain energy criteria:

where \(U\) is the strain energy, and \({U}^{c}\) is the strain energy criteria considered in this study as the unknown material parameter of particular wood species. Numerically, Eq. (15) means that if the strain energy \(U\) is greater than \({U}^{c}\), the element is deleted from finite element computations. Here, the strain energy \(U\) is computed at a level of Gauss points as follows:

In the present model, the element deletion is active from the moment of the blast overpressure peak, that is, \(1.23674\) ms from detonation; the maximal strain energy until this moment is stored, and if it is greater than \({U}^{c}\), the particular element is deleted at the moment of the blast overpressure peak or later.

In general, the elastic stiffness matrix is moisture dependent, but the moisture aspect has not been considered and included in the modelling.

A detailed, three-dimensional finite element model was built to model the contact detonation of wooden members for different wood species (pine, oak and birch) to determine the destructive explosive loads. This destruction is assumed to vanish 80% of the cross-section in the middle of the wood log beam. The dynamic analyses according to explicit FEM code were performed in Abaqus Unified FEA. The geometry used represented the horizontal log of 3 m with three diameters (15, 25 and 35 cm) placed on two perpendicular support logs. Thus, the overall log sample could be considered as a simply supported wooden beam; see Fig. 3—the example for a pine log of 25 cm diameter. The contour plot demonstrates the stored strain energy at time moment \(t = 0.2\) ms. The explosive was assumed to be placed in the middle of the log (point A). The material parameters for a surface explosion of the CONWEP model represent the behaviour of trinitrotoluene (TNT) explosive (Sielicki and Gajewski 2019; Sielicki and Łodygowski 2019; Al-Rifaie et al. 2021). The Explicit Abaqus FEA was used including VUMAT user subroutine package, while the explosion was modelled by Kingery-Bulmash surface blast approximation implemented in CONWEP library. The transversal anisotropic behaviour of wood was fully programmed using FORTRAN coding (VUMAT) including damage criteria (material deletion). In CONWEP feature, the shape of the charge is not included. However, it represents the load for ideally spherical charge with equivalent explosive mass. The constitutive law to model the mechanical behaviour of selected wood species used in FEM computations was described in the previous subsection.

The mesh sensitivity study was performed to determine the reliable size of the mesh. According to previous experience, the mesh sizes assumed were \(1.5\) mm, \(2\) mm and \(5\) mm. The number of finite elements for the biggest logs (\(35\) cm diameter) was 68.1 million, \(36.08\) million and \(2.31\) million elements, respectively for \(1.5\) mm, \(2\) mm and \(5\) mm mesh. The maximal increments were computed by the authors according to the stable time increment size in the explicit FEM code. Thus, the final values were fixed manually to 1e-8, what is lower than stable time increment size, 1.3e-6, 1.7e-6 and 4.4e-6, respectively. The differences between the results for 1.5 and 2 mm models were negligible regarding the destruction of the cross-section for particular TNT charge (here 400 g was used). Thus, 2 mm case may be used as the model with reliable finite element mesh. In Fig. 4, the exemplary results for those models are presented for a pine log (\(25\) cm diameter) with the energy criterion \({U}^{crit}\) assumed to be \(3500\) J. The contour plots demonstrate the maximal stored strain energy. The explosive mass was set to \(400\) g of TNT. The figure shows the significant difference in element deletion between models with different mesh sizes. For the case with larger elements, the deletion was more significant; the difference is not negligible. Thus, in the final computations, a \(2\) mm element was assumed. The finite element of C3D8 was used (according to ABAQUS coding). The C3D8 element is a general purpose linear brick element with full integration scheme.

Due to numerous unknowns in a single FEM model and multiple forward explicit FE analyses to be conducted (see Sect. 3.3), the computations were performed using the Polish high-performance computing infrastructure at Poznan Supercomputing and Networking Center (PCSS), namely, PL-GRID. PL-GRID is a nationwide supercomputer centre for scientific research. This study used an EAGLE supercomputer with 48 CPUs 2.6 GHz Intel Xeon E5-2697 per single analysis. The approximate number of single forward FE analyses performed, considering testing and nominal analyses, was 2160. The wall clock hours of parallel computations were approximately 8640 h in total, using 192 CPUs available (4 analyses in parallel with 48 CPUs in maximum). The main reason for conducting so many analyses was to vary the energy criteria, wood species and log diameters.

Figure 5 presents two strain energy field examples for different energy criteria assumed for pine (25 cm diameter), namely, 3500 J and 6000 J. The contour plots demonstrate the maximal stored strain energy. The difference is not negligible, the sensitivity of the model regarding energy criteria is evident. The longitudinal damage (along the grain) characteristic for wood was obtained; see the arrows.

Based on the measured plots of force vs displacement of birch, oak and pine, the stress vs strain relations were determined according to basic equations of strength of materials. All plot results are included in Figs. 6, 7 and 8; the samples were shifted on the vertical (strain) axis to ensure appropriate visibility (* different shifts were used for each subplot). Each curve represents a single test. The black section for each curve illustrates the linear part of the plot from which the elastic modulus was computed (see Sect. 2.2) by finding the section for which the plot is the closest to the linear regression. In some cases, the phenomenological plastic behaviour may be observed, for example, pine in radial direction. However, because the contact detonation is an extremely dynamic process, it was not included in the modelling of wood.

The statistics of the test output data are presented in Table 3. The mean elastic moduli were computed for each wood direction and species. Moreover, the minimal and maximal moduli are shown, as well as the standard deviation. Note that a low data dispersion was observed. Those mean values of elastic moduli computed were used in the orthotropic constitutive law to model wood; see the Materials and Methods section.

One of the most important aspects determining the wood properties is its moisture. Here, the information on moisture has been included in Table 3 (first column). It should be emphasized that due to the use of green wood samples, which are characterized by high moisture, the obtained values of the modulus of elasticity were low here; to compare with the literature see Table 1, especially in longitudinal direction.

The computations of destructive explosive loads of birch, oak and pine logs according to a military approach were determined in this chapter. Logs of 15 cm, 20 cm, 25 cm, 30 cm and 35 cm diameters were considered. The military approach used was presented in Sect. 2.3. For birch and oak, the results are presented in Table 4; for the pine, the results are presented in Table 5. The plots representing the final masses of the charges obtained are presented in Fig. 9. According to military engineering, the destructive explosive loads for birch and oak are calculated in the same way.

The analysis of the results shows that 25% more explosives are needed to destroy damp wood. Notice that the military approach is an approximation used for a rapid assessment, which allows for quick calculations during field operations.

This paper aims to compare the results of destructive explosive loads for wooden elements according to the traditional military approach with a numerical modelling approach. The finite element method model used for the numerical part was described in the Materials and Methods section. An important question to address regarding the numerical modelling approach adopted is: what damage energy criteria should be used for particular wood species? Selected reasonable values were assumed because these were not explicitly available in the literature.

The numerical model was computed for many different damage energy criteria. The trial and error method was used for first attempts to determine charges for the numerical model, i.e. the charge mass was assumed and its effectiveness was verified by conducting one forward FEM analysis. Only the selected values were presented in the paper; see Fig. 9. For birch, 3500 J and 6000 J were used; for oak, 4500 J, 7000 J, 15,000 J and 25,000 J were used; for pine, 3500 J, 4500 J, and 6000 J were used. For a single energy criterion (i), particular wood species (ii) and log diameter (iii), multiple forward explicit FEM analyses were conducted to find a minimum mass of TNT charge that will destroy a wooden member. The result for a particular selection of (i), (ii) and (iii) is represented by a single square marker on the numerical trend curve; see Fig. 9. The curve was created by a cubic interpolation for three values computed via the FEM approach described above (three square markers for each curve).

As observed in Fig. 9, the military estimation values (dot markers) are usually much higher than the numerical estimation. This discrepancy is especially clear for log diameters of 25 cm and less; for larger log diameters, the numerical estimation gives values nearer to the military estimation. For birch, the numerical estimation is always lower; for larger oak logs, the numerical estimation is slightly higher than the military approach; for larger pine logs, the numerical estimation is almost equal to the military one.

The above conclusions cannot be applied to sapper military engineering unless the damage energy criteria for particular wood species can be determined. Depending on the damage energy criteria, the differences vary and may be significant. For instance, for a 35-cm-diameter oak log, the mass of the charge for energy criteria 25,000 J is more than eight times greater than for 4500 J (3200 g vs 380 g).

The numerical results lead to the conclusion that it is necessary to determine damage energy criteria for particular wood species. Only then may the numerical modelling be used as a robust tool for estimating destructive loads of wooden elements in contact explosions. A future study will be devoted to determining the damage energy criteria for selected wood species in field tests to confirm or reject the preliminary hypothesis that the traditional simplified military approach is severely overestimating the necessary sapper destructive load.