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

In military operations, sappers must often breach wooden structures. The formulas for determining the destructive explosive loads available in instructions and manuals used by sappers are simplified because they consider only a few variables, such as structure member diameter, whether the wood is dry or damp, or the wood species of the structure. In this study, the destructive explosive loads needed to breach pine, birch and oak members were computed via the finite element method. Static compression tests in three directions were conducted to define the orthotropic constitutive models of those wood species, and the results were used as an input to the numerical models. The damage model for wood considered different levels of energy criteria. The finite element analyses of contact explosion of TNT charges against cylindrical log beams were conducted for selected wood species, and destructive explosive loads were computed for different log diameters. Assuming different energy criteria, the results showed that the traditional approach in military instructions and manuals is higher than the values obtained from the numerical approach, i.e., standard manuals suggest using more explosive than may be needed.


Introduction
Despite the overwhelming advantages of more durable and simpler construction materials such as concrete or steel, there are applications and situations in which only wood may be used. One such application is military engineering, which must often use locally available materials due to speed of work and processing simplicity. Wood has complex properties, including anisotropic mechanical properties, hygroscopicity and a natural occurrence of various types of material defects. All of these must be considered when using wood for load-bearing structures.
Material properties for particular applications should be properly identified according to the nature of the load and the phenomena to be considered. Thus, when designing load-bearing structures from wood, the static properties of different wood types should be obtained from static strength tests of standardised samples. Such studies may be easily found in the literature (e.g. Keunecke et al. 2007;Zhou et al. 2017). In those papers, the elastic properties of wood materials were determined by ultrasonic waves. Dynamic tests are more challenging but are still often conducted. For instance, Bragov and Lomunov (1997) described dynamic tests of wood using the wheel method. This research methodology enabled them to determine the characteristics of the material at high strain rates. The study assumed that during plastic deformation, the samples at both ends would be in a state of equilibrium of forces. The theory of elastic wave propagation was used to determine the relationship between stress and strain. As a result of the research, diagrams of dynamic deformations of pine, birch and linden were obtained.
Apart from the static and dynamic properties of wood, other phenomena that may be studied are the damage and fracture of the material. This is especially important when considering military sapper engineering. Some papers deal with this aspect; however, they rarely test raw wood or wood-only structures. Raw woods were considered in Reiterer et al. (2002) and Stanzl (2008). However, the damage/ fracture described is not applicable to the approach proposed in this paper. Wood-based material was analysed in the publication by Hussein et al. (2020a), in which the impact of a blast wave on a probabilistic wall filled with sand and covered with oriented strand board (OSB) was analysed. In another study, by Hussein et al. (2020b), the construction of a wall was considered as an example of a cheap and simple structure protecting buildings against potential terrorist attacks. The study indicated that a sand-plate structure reduces the impact of the shock wave on a potentially endangered object. Blyberg et al. (2012) investigated the combination of various glass binders with wood and considered the degree of damage and types of damage to the adhesive joint. The conclusions indicated that acrylate provides the highest resistance to be damaged.
On the other hand, some damage models of wood may be found in the literature. The paper by Mackenzie-Helnwein et al. (2005) presents a model for softwood suitable for describing inelastic deformations in the plane and transverse to the crust surface. In contrast, the paper by Benvenuti et al. (2020) presents numerical constitutive models that take into account the complex and strongly nonlinear behavior of wood. Numerical models that reflect simulations of the mechanical behavior of a wooden structure, taking into account the dependence on moisture and temperature, are included in the paper by Schmidt and Kaliske (2009). On the other hand, Lukacevic et al. (2017) describe a new approach for predicting the mechanisms of wood failure taking into account microstructural features.
However, there is a lack of blast engineering data in the literature related to specific wood species. Sapper military literature, manuals and instructions, such as Szefostwo Wojsk Inżynieryjnych (1995), provide some rough approaches to determining the explosive strength of wooden structural members, but the results (as discussed in this paper) appear oversimplified. Military manuals or instructions do not consider the anisotropic properties of wood, the moisture of the particular sample and the defect characteristics of the wood. Additionally, the formulas recommended in military manuals represent an empirical approach that has since been applied to estimate the strength of similar wooden load-bearing structures. None of the above-mentioned studies or analyses directly reflect the impact of the explosive material on the wood, including the impact of shock waves or air or contact charges, nor material destruction after the explosion. Thus, a more up-to-date approach is needed.
The numerical approach applying the finite element method (FEM) can provide the answers to questions of what should be the destructive explosive load to breach a green wooden element involving particular wood species. To design a robust model to represent the material behaviour of particular wood species, both the static and dynamic properties should be available, and the resulting damage to the material must be determined. However, before determining the damage properties of a particular material, the correct modelling approach must be established. A sensitivity analysis of unknown damage parameters must first be conducted.
In this paper, the static properties of selected wood species were determined. Furthermore, the destructive explosive loads of those wood species were computed according to military sapper engineering and with a numerical approach through FEM. Birch, oak and pine were considered because those wood species are often used as construction members. This paper presents the modelling approach to represent the orthotropic behaviour of wood subjected to a contact blast load. Furthermore, the damage parameter influencing the FEM results was defined with the proposition of the field test to establish it for particular wood species through a mixed experimental-numerical approach.

Orthotropic and damage properties of wood
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.

Determining the elastic properties of wood
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 × 30 × 30 mm, whereas the longitudinal direction tests used samples of 20 × 20 × 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 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 2 . For the radial and tangential directions, the area was 600 mm 2 . 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. 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).

Defining the modelling of wood-orthotropic constitutive law with damage energy criteria
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).
(1)  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 ij and strains 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: RT = 0.683 , RL = 0.054 and TL = 0.033 . Other Poisson's ratios, namely, TR , LR and 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.

Determining the destructive explosive loads of wood according to the finite element method
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), Fig. 3 Overall view of longitudinal section of a pine log with perpendicular supports presenting strain energy output at t = 0.2 ms after contact detonation of 400 g of TNT (element mesh size equal to 2 mm) Fig. 4 Magnified longitudinal sections of pine models presenting maximal stored strain energy at t = 2 ms and after contact detonation of 400 g of TNT with two selected mesh sizes: a 5 mm and b 2 mm (black cross represents the explosive) 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.

Elastic properties of selected types of wood
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

Destructive explosive loads of selected wood species according to a military approach
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.

Destructive explosive loads of selected wood species according to the finite element method
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,  1 3 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.

Conclusion
In this study, birch, oak and pine samples were considered for computing the destructive explosive loads of wooden members for sapper engineering purposes. Destructive load curves from the traditional military calculations and advanced numerical estimations were compared for three wood species and different diameters of wooden beams. The finite element method was employed using an orthotropic constitutive law of wood with a damage model for numerical Fig. 9 Destructive mass of TNT for a birch, b oak and c pine members according to military approach (as calculated) and numerical estimation; see Tables 4 and 5 and Eqs. (1-2) 1 3 estimations. The wooden material model was fed with the experimental static data from compression of wood samples, tested in three directions.
As a result of the work, the following issues may be identified: -Experimental orthotropic properties of the selected woods were presented, as were their relationships in three directions of the samples to be used in the numerical models -Traditional military calculations greatly simplify the numerical calculations of destructive loads of wooden members and do not consider many properties of the structure or the material -Moisture has an important influence on the strength of wood, although this influence was not explored in the paper -The influence of damage energy criteria on the destructive loads of wooden beam members in the numerical estimations is significant and was shown in this paper -There is a need for more practical tests in the field using wooden beam samples subjected to explosive charges and designed to determine the damage energy criteria for particular wood species. Due to the current lack of such knowledge, the aim of future research will be to determine the damage energy criteria for selected wood species.