Neural Networks for Predicting Fracture Toughness of Individual Wood Samples

Strain energy release rate (GIc) of Pinus radiata in the TL opening mode was determined using the compliance crack length relationship. A total of 123 specimens consisting of four sizes of specimen with each size having four different crack lengths were tested. For each specimen, grain and ring angles, density and moisture content were measured. Video imaging, was used to measure crack length during propagation. Since cracks extended in stages, full compliance-crack length relationship was developed for each specimen based on their initial and subsequent crack lengths. No significant differences in GIc, between initial and subsequent crack lengths were found for the smaller specimens by paired sample t-tests, but differences were significant for the largest specimen size. The Average fracture toughness was calculated from GIc and it was 215 kPa.m0.5. Three artificial neural networks were developed to predict the: 1) force required to propagate a crack, 2) crack extension, and 3) fracture toughness of an individual specimen. Each was successful, producing respective R2 of 0.870, 0.865, and 0.621 on validation data. A sensitivity analysis of the networks revealed that the crack length was the most influential with 21% contribution followed by grain angle with 14% contribution for predicting the applied force. This was followed by volume and physical properties. For predicting the crack extension, density had the greatest contribution (20%) followed by previous crack length and force contributing 16% equally. Fracture toughness was dominated by the dimensional parameters of the specimen contributing (42%) followed by anisotropy and physical properties.


Introduction
Wood contains inherent flaws such as knots, holes, notches, and pre-existing cracks some of which are intentional. These can cause failure to occur below the normal failure stress of unflawed wood due to the stress concentrations. Therefore, there is considerable interest in the development of non-destructive test methods for accurately assessing strength of wood (Porter 1964).
The American Society of Testing and Materials (ASTM 2005) provide an empirical relationship between lumber strength and knot size. However, this ignores the effect that pre-existing cracks and other such flaws could have on timber strength. It can be used to predict, with reasonable accuracy, the strength distribution of graded timber but gives no indication as to the strength of an individual member. It has now become accepted that fracture mechanics can and should play a major role in these timber engineering problems (Valentin et al. 1991).
The following study employs the energy balance-based fracture theory and video imaging to determine strain energy release rate that allows to obtain fracture toughness of wood. Furthermore, the study attempts to use neural networks to model fracture properties to obtain an estimate for individual specimens by incorporating their variation in physical and geometric properties.

Crack Propagation and Stress Concentration
Fracture Mechanics relate parameters such as material properties, flaw geometry, and applied loads to the resulting stress concentrations surrounding a crack tip. The linear-elastic fracture mechanics (LEFM) is the basis for most fracture studies today. Wood is an anisotropic, inhomogeneous and non-linear elastic material. For these reasons, the determination of stress/strain relationships throughout the material is difficult. In fact, Samarasinghe and Kulasiri (2000a and b) demonstrate, graphically and mathematically, displacement profiles of stressed wood surfaces and near-tip regions of cracked bodies obtained from digital image processing and these clearly illustrate the nature of heterogeneity of local displacement profiles. It is for these reasons that Griffith's (1920) fracture theory based on energy required to fracture has so much appeal.
According to LEFM, the critical stress σ cr required for growth of a crack of length 2a in a plate under tension σ can be expressed in the form of: where E is the Young's modulus and G c is the critical strain energy release rate (twice the specific surface energy) and k c is the critical stress intensity factor or fracture toughness. For structural members, three stress-intensity factors, k 1 , k 2 , and k 3, are identified by their respective modes of crack extension (Mode I, II, III). The stress intensity factor is an indication of the level of stress surrounding a crack for each respective crack mode. Fracture toughenss, k c , and critical strain energy release rate G c , are two related fracture strength parameters of a material. Thus, the mode I (opening) fracture toughness, k 1c, is related to the critical strain energy release rate G 1c as: Although the two fracture strength parameters are related, the more fundamental energy release rate concept based on total energy in a body can provide insight into fracture phenomena (Porter 1964) along with realistic values for the fracture strength, especially for heterogeneous materials.

Energy Release Rate and Fracture Toughness of Wood
Prior to Porter (1964), very few fundamental studies had been performed on the manner in which fracture occurs in wood. Atack et al. (1961) measured the surface energy of green wood. Wu (1963) developed mixed mode fracture criteria using orthotropic balsa wood plates. Porter (1964) concluded that the combination of an energy release rate method and experimentation allowed considerable insight into the mechanisms by which wood fails. Since then, there has been a great deal of research in the field, with publications of both fracture toughness values and development of the related experimental procedure and theory (Patton-Mallory andCramer 1987, Conrad et al. 2003). Gustafsson (1988) demonstrated that simple design expressions can be obtained by balancing the change in a system's strain energy with the energy required to form a fracture surface. Boström (1992) combined fracture energy observations with a "Fictitious Crack Model" to accurately predict the behaviour of ASTM compact tension specimens and to model formation of drying cracks in lumber.
Generally, energy methods show great promise for situations where the fracture process zone near the tip is not small compared with the length of the crack. This opens the door for accurate prediction of fracture behaviour in various wood structural components including connections and to make reliable interpretation of results from "standardised" fracture property tests (Smith and Chui 1994).
The experimental technique adopted in the present study to determine the strain energy release rate is an adaptation of that employed by Porter (1964) and Triboult et al. (1984). It requires the determination of the change in stiffness of the specimen as a crack propagates. For a plate with a crack of length a and thickness B stressed with a load P, the strain energy release rate is given by where C is the compliance, C = Δ/P, where Δ is the crack opening displacement.
Thus, the strain energy release rate can be determined with known compliance-crack length relationship, the applied load, and specimen width. (The work described in our paper uses this approach with a great deal of attention being paid to the characteristics of compliance.) For orthotropic material such as wood, the mode-I critical stress intensity K IC is related to G IC by (Sih et al. 1965): where S* is the inverse of an apparent elastic modulus given by The relationship in Eq. 4 is valid for orthotropy, plain strain conditions and linear elastic behaviour (Sih et al. 1965). Triboult et al. (1984) undertook to investigate the validity of these assumptions and concluded that they were in fact reasonable for wood. However, most wood used in service deviates from orthotropy as indicated by non-zero grain and ring angles that can alter the fracture properties. In fact, Samarasinghe andKulasiri (1998, 2004) demonstrate the predominant influence of grain angle on Mode-I fracture toughness of New Zealand Pinus radiata. Furthermore, wood density within a species varies with a number of factors including location in a tree, location within the geographical range of the species, site conditions, and genetic sources. However, the designer of a wood product has no control or knowledge of where the particular tree was cut, in what part of the tree it originated, or necessarily even the particular species.
All wood properties are affected by the moisture content. For instance, when testing over the possible range of moisture contents, double edge-notched specimens display ultimate stress intensities that differ, for some woods, by up to 30% (Ashby et al. 1985). These effects further complicate the prediction of wood properties for a particular specimen. Thus, the fracture resistance that is considered to be a material property has high variability and a complex relationship with the environment and geometry. Therefore, approaches that incorporate measurable physical and geometric properties for predicting fracture strength are urgently needed and these could significantly improve the efficiency of material use. The present study is a step towards this direction.

Objectives
The goal of the study is two-fold: 1) To determine mode-I (TL) fracture strength of Pinus radiata using energy release rate measurements obtained from crack length extension measurements made by using video imaging; and 2) To use neural networks to predict fracture toughness, peak load, and crack extension for individual specimens with known geometric and material properties. This involved following sub-objectives: 1) Determine energy release rate in crack extension through compliance measurements obtained for various initial crack lengths using video imaging and investigate in detail characteristics of energy release rate; 2) Develop neural network models for predicting fracture toughness for individual specimens from its easily measurable physical and geometric properties and study the influence of these properties on fracture toughness; and 3) Develop separate neural network models to predict maximum load required for crack propagation for individual specimens and the resulting crack extension from their physical and geometric properties.

Specimen Geometry and Characteristics
Specimens were prepared as shown in Fig. 1 from kiln-dried knot free New Zealand grown radiata pine (Pinus radiata) boards obtained locally from a sawmill in Christchurch. There were four specimen sizes and each size contained eight sets of four differing initial crack lengths, as indicated in Table 1. There were a total of 128 specimens but tests were successfully conducted only on 123 test specimens. Cracks were made using a sharp knife-edge especially prepared for the purpose. Specimen size was varied in order to ascertain geometrical dependency of G IC and K IC . Crack length was varied for each size in order to obtain a compliance value (C = Δ/P) for each crack length and use these to establish the compliance crack length relationship (Ca) for each specimen size. The partial derivative (∂C / ∂a) of this relationship was to be used in Eq. 3 to obtain strain energy release rate. The prepared specimens were left in the laboratory for over six months in order for them to become accustomed to the local climatic conditions yielding an equilibrium moisture content of around 8%. They were stored in a plastic sheath until required for the test.
Before testing, all the dimensions, including original crack length, were measured to ±0.5 mm accuracy. The grain angle was measured from the end of the crack with an estimated accuracy of 1 degree using a protractor. The load direction was not exactly in either radial (R) or tangential (T) directions. Therefore, the load direction in the RT plane was determined by measuring on the cross section the angle (ring angle) the tangent to a l w h Fig. 1. Geometry and configuration of specimens tested for the determination of strain energy release rate. a late wood ring makes with the original crack. If the ring angle is zero, load direction is radial and if it is 90º, load is in the tangential direction.

Experimental Procedure
All testing was undertaken using a Sintech 30/D testing machine run by "Testworks" material testing and analysis software as shown in Fig. 2. Each specimen was secured to the machine by means of pins that passed through its width on either side of the crack and the other end was simply supported. The load was applied to these pins through a universal joint, which ensured lateral restraint was not a contributing factor to the results. Each specimen was painted black in order to ensure that the camera captured the full extent of the cracking. As the crack propagated, the white flesh of the wood was exposed and contrasted against its black exterior.
For the purpose of acquiring crack extension data, a Hitachi VM-S7280E colour TV camera was set up on a tripod to film the specimens during the testing. The crack extension images taken at appropriate times were transferred in real time to a computer. A timer was located in view of the camera to synchronise the video recording with the time recorded by the "Testworks" software in order to obtain the loads at various crack events from load-time graphs.
Following procedure was followed for each test. After securing the specimen, the camera was focused on the specimen, and it was ensured that the line of sight was perpendicular to the specimen. Prior to testing, an image of a graph paper set against the specimen was taken to calibrate the pixel dimensions to actual dimensions. Then the camera, timer and the testing machine were started simultaneously and the test was run until the crack propagated to the end of the specimen. The load-displacement-time relationships were also recorded during loading. After the test, load/ time graph was used to select appropriate crack length images.
After testing, several pieces were cut from near the center of each specimen and moisture content was determined using the oven-dry method and dry density was obtained from the volume determined by immersion in water of dry specimens coated with wax.

Crack Length Measurement
The method of energy release rate requires the measurement of crack opening displacement (C.O.D.) and crack length. The C.O.D. between the pins was obtained from the images, testing machine and from an electronic calliper. From repeated measurements, a perfect linear relationship between these measurements was found with an R 2 value of 0.97. The crack length was remotely measured at pertinent points during propagation from the desired digitized individual images of 512 × 512 pixels. The length of the crack was measured in pixels and converted to mm by using the pixel/mm relationship established previously.

Determining Compliance/Crack Length Relationship
In order to determine the strain energy release rate G IC , it was necessary to find a suitable expression for the relationship describing the partial derivative of the change in compliance as the crack extended: ∂C/∂a. This requires determination of compliance C (= Δ/P) for various crack lengths (a) and relating C to a. The original idea was to test specimens with different crack lengths separately and use them to obtain the whole compliance curve but during the testing it was found that a single specimen can be used to obtain compliance for several crack lengths. This was because the initial crack propagated in quite distinct stages making it possible to measure several definite crack lengths on each specimen and the corresponding load required to extend them, without even having to unload and reload the specimen after each crack extension. This was equivalent to testing several specimens with different initial crack lengths. This is novel and significant because it not only produces realistic compliance curves but also allows to determine several measurements of strain energy release rate and therefore, fracture toughness, from a single specimen. This approach is different to that described by Conrad et al. (2003) which involved unloading and reloading the specimen after each crack extension. At each crack extension, compliance was calculated and the values for a single specimen that had 5 crack extensions is shown in Fig. 3 and an exponential function shown in Eq. 7 and depicted in Fig. 3 accurately modelled this data. In Eq. 7, A and B are constants and a is the crack length.
All of the specimens that produced sufficient crack length data (3 crack propagations and above) were graphed in this manner, and a similar trend was seen in each case: all had an R 2 value of above 0.9. Compliance-crack length relationship for four specimens with similar geometry (Category 2 in Table 1) but different initial crack lengths (IC) is shown in Fig. 4 from which an exponential approximation shown in the figure was determined to be reliable with an R 2 value of 0.9164. The individual compliance crack length curves in Fig. 4 overlap and illustrate that the initial length of the crack is of little or no consequence to the determination of the compliance/crack length relationship. It was also attempted to develop a single relationship for all 32 specimens of similar geometry (Category 2) as shown in the plot in Fig. 5.
However, the divergence observed for the larger crack lengths (initial and subsequent) in  . The exponential nature of the compliance -crack length curve illustrated using data from geometrically similar specimens (Category 2 in Table 1) with varying initial crack lengths.
produced a regression with a low R 2 value of 0.546, rendering this grouping of little use in analysis as it would be a poor approximation to the compliance/crack length relationship for individual specimens. This divergence is attributed to the lack of homogeneity between specimens. In order to overcome this, the individual compliance/crack length relationship created for each specimen in the form of an exponential function was used for the subsequent analysis. The mean and standard deviation of constants A and B for each geometry are shown in Table 2 and these were found through linear regression. By differentiating Eq. 7 with respect to crack length, rate of change of compliance was obtained and substituted into Eq. 3 to obtain strain energy release rate for each crack length on each specimen. Since several crack length measurements were made on each specimen, there were many more measurements of G Ic than there were specimens.

Determining Strain Energy Release Rate
Strain energy release rate was determined for the initial as well as the subsequent stable crack lengths. Fig. 6 shows the trend observed for the average initial G IC values when plotted against its associated initial crack length for the largest specimen geometry (see Table 1). It shows an increasing trend of G Ic with the initial crack length.
When all the measured G IC values are grouped for each geometry and plotted against crack length however, the trends become difficult to observe as shown in Fig. 7 which represents the largest specimen geometry. Both individual specimen G IC from initial and subsequent crack lengths displayed such scatter in the data. A similar scatter in G IC was observed by Triboulot et al. (1984) in their study using only the initial crack lengths but in their study they did not reveal the existence, in this scatter, of much more coherent and meaningful relationships in the case of individual speci-   Fig. 6. The average initial strain energy release rate for initial crack lengths for the largest specimen geometry.  Table 1).
mens, as clearly demonstrated in this study. The average G Ic value for the data set is represented by the hyphenated line in Fig. 7. Average properties are normally used currently in practical decision making and the figure highlights the deviation of individual specimen values from the mean due to the heterogeneity of properties of the specimens. Fig. 7 demonstrates the highly variable nature of the strain energy release rate when all the values are combined for just a single geometrical size.
If values for all four sizes were combined, even more scattered data plot would result.

Statistical Comparison of G Ic for Intro duced (Initial) Cracks and Subsequent Naturally Occurring Cracks
Since it was possible to determine G Ic from the initial introduced crack length as well as from subsequent extended crack lengths on each specimen, it is possible to test if these two sets of values are similar. The difference between the initial and extended cracks is that the initial crack was artificially created to resemble a natural crack and extended cracks were naturally occurring cracks due to a propagating crack. There were 3 to 5 crack extensions per specimen and thus as many G Ic values per specimen. No difference was apparent between G Ic for initial and subsequent crack lengths; therefore, statistical tests of significance were conducted and for this purpose, all G Ic values for each geometry were grouped into initial (i.e.1st), 2nd, 3rd, 4th, and 5th crack length samples. Although these subsequent crack lengths were not identical, they were grouped on the basis of the stage of crack extension. A paired sample t-test (two-tailed) that compares the means of two groups for a single variable was undertaken in SPSS 10 for Windows to determine whether there is a significant difference between the observed strain energy release rate for initial and subsequent crack lengths. It assumes that the observations for each paired group are normally distributed and to test the distribution of the data groups for normality, a skewness test was undertaken, also in SPSS 10 for Windows. Assumption of normality is violated if the ratio of skewness to its standard error is less than -2, or greater than +2. Out of 19 groups only 3 were skewed.
A paired samples t-test (two-tailed) was computed for those 16 data groups that were deemed normally distributed. A 95% confidence interval was selected, and so a t-test returning a p-value greater than 0.05 determines that there is no significant difference between groups' means. Tests were conducted comparing G Ic(initial) with G Ic for each extended crack (i.e., 2nd, 3rd ,4th, etc.) as well as comparing G IC for extended cracks with each other. There was no significant difference either between the G Ic for initial and any of the subsequent crack lengths, or between the G Ic for subsequent crack lengths for the two smaller geometric categories in Table 1 indicating that the artificial cracks and natural cracks have similar behaviour for these two specimen geometries. In these two cases, there were a total of 10 comparison tests.
For the third largest geometry, there were 6 comparisons tests and all but one were insignificant. A significant difference was found only between G Ic(initial) and G IC for the 2nd crack length. However, since G Ic for initial and 3rd (and 4th) crack lengths were insignificant, this result is inconclusive. More important for this case is that 5 out of 6 tests indicate that there were no significant difference between the initial and subsequent crack lengths and that subsequent natural cracks display similar behaviour. As for the largest size category, there were 10 comparison tests and significant differences were found in four tests between the G Ic for the initial crack length and G Ic for all the subsequent (2nd), 3rd, 4th, and 5th crack lengths. However, there was no significant difference between the G Ic for any of the four subsequent crack lengths indicating that all the natural cracks have similar behaviour in this geometry as well. These results indicate that for the largest specimen geometry, the made crack and natural cracks have different characteristics.
Comparisons were also made between the mean of the initial crack lengths' energy strain release rates and the combined mean G Ic of all subsequent crack lengths for each size category, and again their differences were insignificant for the two smaller sizes. There was a significant difference for the third largest category; however; this is due to the significant difference between just one pair of G Ic (for initial and 2nd crack length) out of six pairs as already discussed. There was a difference between the mean G Ic for initial crack length and the mean value for all subsequent crack lengths for the largest size category. The results indicate that an introduced crack may be a good approximation to a naturally occurring crack for a range of smaller specimen sizes but not for larger specimens. These mean G Ic results for the four size categories are presented in the next section. Fig. 8 shows the relationship of G Ic to volume. The top figure is for average initial G Ic and the bottom figure is for the average G Ic for the subsequent (moving) crack. Fig. 8a and b show an increasing trend of G Ic with volume. G Ic 's volumetric dependence appears to have more effect on the specimens with a smaller volume, as indicated by the logarithmic nature of the G Ic vs. volume trends. The two figures also highlight the validity of the assumption that the G Ic of a pre-fabricated crack (i.e., initial cracks) is a good experimental approximation to that for a naturally occurring (subsequent) crack. For instance, the average G Ic values for two smaller sizes are similar and divergence just sets in for the 3rd category and becomes prominent for the largest category.

Geometric Investigation of G Ic
Average strain energy release rate can be used to obtain fracture toughness from Eq. 4 with appropriate material properties in Eqs. 5 and 6. Since elastic moduli for individual specimens were not available, average material properties obtained for the same material in a previous study (Bandara et al. 1999) were used and these are shown in Table 3. The average strain energy release rate and fracture toughness values for each specimen size is shown in Table 4 for initial crack lengths (G Icinit and K Icinit ) and for all crack lengths (G IcAll and K IcAll) . These results are what can be expected for this material (King et al. 1999, Samarasinghe andKulasiri 2004).

Physical and Geometric Variables and Fracture
There are multiple variables involved in wood fracture dynamics. These relate to anisotropy and physical and geometric properties of wood and characteristics of cracks. For example, Samarasinghe and Kulasiri (2004)    (angle between the crack plane and the direction perpendicular to a growth ring in the RT plane on the cross section) on K Ic of Pinus radiata. They found a positive relationship between the two with correlation coefficient r of 0.75. This is a result of anisotropy due to geometric axes not coinciding with the material axes as assumed in pure orthotropy. Fonselius and Riipola (1992) also found a significant effect of anisotropy in their study in which fracture toughness in Mode-I of standard compact tension specimens showed a strong positive correlation with the orientation angle of load in the RT plane, measured as the angle between a tangent to a growth ring and original crack plane. They also found density to be the next most important variable for K Ic and developed linear regression models to predict fracture toughness from the orientation angle and density. The size effect was insignificant in their study for thickness ranging from 12 to 50 mm. They also found that it is insensitive to constant moisture conditions corresponding to relative humidity ranging from 50% to 80%. They suggest that it is not the moisture content but the moisture gradient that influences the fracture toughness. Further confirming the effect of anisotropy, Smith and Chui (1994) found from SEN (single edge notched) beam specimens of red pine that G Ic increases with orientation angle increasing from 0° up to 60° beyond which it decreases sharply. Porter (1964) in his study on strain energy release rate G Ic of western pine found that it is independent of crack length, specimen height, length and width, where the latter varied from 3 mm to 25.4 mm. This led to his assertion that G Ic is truly a material property. However, different studies have seen that different size effects, such as volume, height, or height*length, predominate (Bandara et al. 1999). Porter (1964) found that G Ic increases with moisture content increasing from 0% to about 18%. Similar trend was found for G Ic for red pine from single edge notched beam specimens by Smith and Chui (1994). They found that G Ic increases with moisture content from 7% to 18% and then decreases from 18-24%. However, this effect cannot be directly translated to fracture toughness as the latter is related to G Ic through Young's modulus which also varies with moisture content. After an extensive review of literature on fracture of solid wood, Conrad et al. (2003) state that the current knowledge of the subject indicates that fracture toughness increases with increasing density and strain rate and decreases with the number and size of defects, and that it is maximum around 6-8% moisture content. King et al. (1999) report that opening mode fracture toughness in K IcTL (crack propagating in L direction) for wet and dry New Zealand Pinus radiata does not differ significantly for SEN bending specimens. In all the other crack planes except TR, opening mode fracture toughness for dry wood is larger. Petterson and Bodig (1982) from an extensive study of fracture in mode-I for a number of softwood species in a range of moisture content up to 30% developed a formula to compute fracture toughness from known moisture content and specific gravity of green wood. In their study, fracture toughness is linearly related to green specific gravity and exponentially and negatively related to moisture content up to fiber saturation point. These studies have incorporated a range of species but have not investigated the relations in natural varying properties in a batch of kiln dried wood of the same species with the general assumption that the variation is too great to yield any significant trends.
The moisture content and density variation among specimens from the same kiln dried batch can be large enough to make their properties differ significantly but currently, there are no methods to integrate these values of individual specimens into a prediction of G Ic or K Ic of individual specimens. These results point to the fact that a complete picture of the combined effect of the physical, geometric, and crack properties on wood fracture is still lacking in description. Table 5 shows the mean and the standard deviation for the physical and geometric variables measured in our study and these variations must explain the scatter in the results such as G Ic values observed in Fig. 7. The relationship shown in Fig. 9 was found between G Ic and the ring angle and this trend is compatible with the findings of Samarasinghe andKulasiri (1999, 2004). The volume and crack length effects have already been shown in Figs 6, 7, and 8. For the other measured variables, the scatter was too great to discern obvious trends.
Past methods for determining fracture toughness (including those outlined in the previous section) have involved large spreads of data and relatively inconclusive and non-applicable results.
In an attempt to reduce this large spread of results and produce a good predictive tool for fracture toughness for individual specimens of wood, research was undertaken to include some of the afore-mentioned specimen parameters. However, there is little agreement, if any, between researchers on the relative contribution of these variables. Judging by the values in Table 5, specimen geometry, density and anisotropy that show evidence of spread must be the factors contributing to the scatter in the data.
It is evident that producing a suitable model to predict fracture toughness is a very difficult task. As discussed previously, fracture takes place in wood when all the necessary parameters reach a threshold simultaneously. The greatest difficulties arise, however, not from trying to determine these thresholds, but from determining their interrelationships. It is this task the current study attempts to accomplish using Artificial Neural Networks, a powerful nonlinear data modelling method.

Artificial Neural Networks
This emerging technology is rooted in various disciplines and the present author provides a comprehensive introduction for scientists and engineers in "Neural Networks for Applied Sciences and Engineering" (Samarasinghe 2006). One feature among many of Neural Networks' capabilities is their ability to find complex nonlinear interrelationship among many variables that produce an outcome. The concept of artificial neural networks were inspired by biological neural networks that consist of many neurons (nerve cells) that process information in the brain. A biological neural network is a nonlinear, highly parallel system characterised by robustness and fault tolerance (Samarasinghe 2006). The features which an Artificial Neural Network (ANN) attempts to capture are: learning by adapting to changes in the surrounding environment (trends); handling imprecise, fuzzy, noisy, and stochastic information (variability); and generalising from known tasks or examples to unknown ones (robustness).
An artificial neural network (Fig. 10) consists of many artificial neurons that process their inputs and send the output to one or many neurons con-  nected to them until the information propagation is complete and the network produces an output. They are trained to learn the relationships in data by repeatedly presenting examples (input-desired output data) to it and adjusting internal parameters until the predicted outcome is as close to the desired outcome as possible. The layout of the most popular neural network, called Multilayer Perceptron (MLP), is shown in Fig. 10 where neurons are organized into layers where the input layer represents input variables and the output neuron represents the output variable. The middle layer consists of so called hidden neurons that make it possible to model highly nonlinear relationships between inputs and output. In this network, data pass from the input layer to the output layer. A weight (w ij ) is the weight on input i feeding neuron j and indicate the relative influence of the corresponding input i in producing the outcome. Information processing in an individual hidden neuron involves first summing (∑) the weighted inputs. The hidden neuron j further processes this weighted sum of inputs to produce a nonlinear output using a nonlinear activation function (σ). The output is weighted by the hidden-output weights w j and sent to the output neuron where a similar process is repeated and an output is produced.
The neuron activation functions can be linear, logistic, Gaussian, Sine etc. as shown in Fig. 11 and it is these that give a network its powerful nonlinear mapping capabilities accomplished through local mapping done in individual hidden neurons.
Mathematically, the input-output relationship of the networks can be expressed as: which represents the neural network model. The network goes through an iterative learning process using the input-output data to determine weights w ij and w j and the final weights represent the coefficients of the model. As in regression, mean square error between actual and predicted is used as the measure of the goodness-of-fit of the model. The model must then be calibrated and validated. From a calibrated model, various issues related to the problem, such as sensitivity of output to inputs, relative contribution of inputs to the output, stability of coefficients etc., can be assessed. This type of network has been used in our study presented in this paper to predict fracture properties for individual specimens from their geometric and physical properties. Seibi and Al-Alawi (1997) successfully applied neural networks to predict fracture strength of steel.

Implementing Artificial Neural Net works for Modelling Fracture in Wood
Three ANNs were developed in this study with the intention of: 1) Predicting the applied load required to initiate crack propagation in an individual specimen, 2) Predicting the distance over which a crack will propagate in an individual specimen, and 3) Predicting the fracture toughness of an individual specimen.
The predictions for individual specimens are made from their easily measurable physical and geometric properties (Table 5), such as dimensions, density, crack length, load orientation and grain angle. The corresponding output was either Peak load, Crack extension, or Fracture toughness, for the three networks. The property variations considered are natural variations that can be found in a batch of kiln-dried wood. The networks were developed on NeuroShell 2 (1997).

Data Preparation
Separate input and output variables were defined as appropriate for each ANN and these are presented in Sections 6.3.4, 6.3.5 and 6.3.6 where actual networks are developed. Inputs were scaled between -1 and 1 to maintain a similar range for all variables so that the sheer magnitude of one variable does not mask a significant influence of another variable. Then the data for each network was divided into three sets: the training set, test set, and validation set, consisting of 70%, 15%, and 15% of the data, respectively, as shown in Table 6. All sets where chosen randomly, and then the individual data records were randomly ordered in the data sets to eliminate any bias.
It is the training set that was used to adapt network weights. The test set was run through the network periodically during training to ensure maximum generalization by avoiding overfitting. Overfitting sets in when the networks tries to model noise as well. During training, the error on both training and test sets decrease until overfitting sets in at which point, the error on the test set starts increasing while training error continues to decrease. The network weights are optimum at this point and the training is stopped. Finally, the validation dataset is used to test the trained model.

Training Networks
Training commenced with setting initial weights to randomly values extracted from a uniform distribution and termination criteria to mean square error (MSE) set to 0.0002. Training is a quite an involved process whereby many aspects of a network are modified. These are number of hidden neurons, type of nonlinear activation function in neurons, testing the capability of different learning methods available, retraining the network with different initial weights in order to avoid local minima in the error surface. All of these were done in this study. For each ANN, various network configurations and training methods were tested until the best network configuration (i.e., number of hidden neurons and activation functions) was found through validation. Multiple Linear Regression was also conducted for each case for comparison.

Validation
Validation was undertaken by running the previously withheld validation data set through the network, and calculating R 2 -the coefficient of multiple determination. For neural networks, it compares the accuracy of the network model to the accuracy of a trivial benchmark model that only predicts the mean of the sample. A perfect fit would result in an R 2 value of 1; a very good fit, near 1; and a very poor fit, less than 0. If the neural network model predictions are as accurate as could be predicted by the mean of the sample, then the R 2 value would be 0. This point is particularly important for wood because in the absence of methods to predict the required properties for individual specimens in a batch of kiln dried wood, mean is usually used as an estimate.
Thus, if a model yields R 2 values larger than zero in predicting the value of a certain property for an individual specimen, then the model prediction is better than using the group average to represent this property for that individual specimen.

Predicting the Force Required for Crack Extension
Of the 123 specimens fractured, it was possible to determine 403 separate data records in which an applied load produced crack propagation. The output and input variables for this network are listed below.
This network predicts the peak load required for crack extension of an individual specimen from its material and geometric properties and initial crack length. The architectural specifics for the best network are as follows: 3 Layer MLP (Fig. 10); 9 Inputs; 23 Hidden neurons; and 1 Output neuron. The 23 hidden layer neurons had Gaussian Complement transfer functions. This is the inverse of the Gaussian function (Fig. 11) and is sensitive to extreme values of inputs. The output layer has one neuron with a logistic transfer function.
The network predictions for the validation set, presented in conjunction with the actual output data form the same set, are shown in Fig. 12 which indicates good correspondence. Various other measures, such as errors, correlation, and plots between the actual and predicted output, were assessed to further validate the results for the network (This extra validation procedure was used for the other two networks as well). The network produced an R 2 value of 0.865. Multiple linear regression model resulted in R 2 0.80 which is also high but still less than that produced by the best neural network. This indicates that the prediction of crack extension force is a relatively linear problem.

Predicting Crack Extension Length
This section addresses prediction of crack extension length for a given peak force applied to an individual specimen. Of the 123 specimens that fractured, it was possible to determine 288 complete data records in which there was both a previous crack length, a (n-1) , and a new crack length, a n . The output and input variables for this network are presented below: Output: Crack lengtha n (mm) Inputs: Specimen length (mm), Specimen height (mm), Specimen width (mm), Grain angle (degrees). Ring angle (degrees), Moisture content (%), Density (kg/m 3 ), previous crack lengtha (n-1) (mm), Applied force -P n (kN) This network predicts the new crack length of an individual cracked specimen from its physical and geometric properties, initial crack length, and peak load. The optimal network architecture found is listed below and it is similar to the one developed previously (i.e. 3 Layer M.L.P, 9 Inputs, 23 Hidden neurons with Gaussian Compliment, 1 Output node with Logistic function); however, the network weights would be different. It produced an R 2 value of 0.870 when the validation set was run through it. The Multiple linear regression produced an R 2 value of 0.588 on the validation data which is much lower than that given by the best neural network. This indicates that crack extension is a highly nonlinear problem. The network predictions for the validation set, presented in conjunction with the actual output data form the same set, are shown in Fig. 13. It shows a very good correspondence between the actual and predicted crack lengths.

Predicting Fracture Toughness of an Individual Specimen
For modelling fracture toughness, accurate fracture toughness for each specimen must be available. Using G Ic with values for the required material properties in Eqs. 4-6 would produce accurate values for fracture toughness of individual specimens. Since all these properties were not available for each specimen, the fracture toughness for this exercise was obtained from the standard formula (Tada et al. 1985) in Eq. 9 to improve the accuracy of prediction by the neural network. It basically uses the left hand side of Eq. 1. The objective of this modelling exercise is to see whether the variation of individual specimen characteristics can account for variation in fracture strength and therefore justifies this approach to determination of K Ic . The fracture toughness of each specimen at each point of crack propagation was determined. Of the 123 specimens fractured, it was possible to describe 216 data records in which a full set of input variables could be related to a fracture toughness, K Ic , given by where, P max is load at crack extension, a is crack length, w is specimen width, h is specimen height and ν is poisson ratio in LT plane (0.45). The bracketed expression on the right hand side is a geometry correction factor and Eq. 9 applies to plane strain conditions. The optimal network architecture found was: 3 layers; 8 Inputs; one hidden layer with two sets of 7 neurons with one set having a hyperbolic tangent (tanh) and the other Gaussian transfer functions; and Output node with Logistic activation. This produced an R 2 value of 0.621 when the validation set was run through it. The network predictions for the validation set, presented in conjunction with the actual output data from the same set, are shown in Fig. 14. Note that the units used are MPa.m 0.5 . The Multiple Linear Regression results were much poorer with an R 2 value of 0.20 indicating that the prediction of fracture toughness is a highly nonlinear problem. The above results for accuracy of prediction of fracture toughness is significant because the mean value, that is normally assumed in practice, corresponds to an R 2 value of zero (i.e., inputs do not influence the output). The linear regression state that inputs account for 20% of variation of K Ic and it still cannot explain 80% of the scatter in K Ic . The neural network has captured 62% of variation in K Ic as being explained by the inputs when it predicts K Ic from inputs not seen during training. This indicates that inputs do affect K Ic and this effect is nonlinear.

Determining Influence of Geometric and Physical Variables on Fracture Properties
ANNs develop a nonlinear mapping function between input and output variables. It can be used to determine the sensitivity of the output to certain variables or groupings of variables. This is of great interest in the case of fracture toughness, as it allows comparison of the influence of several physical and geometric variables. A sensitivity analysis was conducted on each of the networks to ascertain the relative contributions of each input variable on NeuroShell 2. Tables 7 through 9 show the relative contributions in the three networks investigated above. The left-hand column lists the variable, and the right lists its relative contribution to the output as a percentage of all the other variables' contributions. It gives the rank order of the inputs in terms of the significance of a variable in predicting the output. Table 7 displays the relative contribution of the input variables to predicting peak load for crack extension. Table 8 does likewise for the network  predicting the corresponding crack extension, and  Table 9 for fracture toughness. Table 7 shows that the peak load required for crack extension is largely governed by geometric properties with crack length being the most predominant which is intuitive. The grain angle is the second most influential variables confirming the results by Samarasinghe and Kulasiri (2004). These two variables together contribute 35% to the prediction. Another 30% is contributed by the volume. The physical properties -density and moisture content -together contribute 20%. For the range of data for the ring angle in this study, its influence is about 9%. As for the volume, height is the most prominent followed by width and length.
In predicting crack extension, density, initial crack length and applied force paly a major role contributing 52% to the prediction (Table 8). Den-sity is the most influential variable. The volume has about 26% influence, with the height again being the most predominant followed by width and length. The anisotropy contributes 20%.
As for fracture toughness, Table 9 indicates that volume is predominant (42%) with height as the most influential variable followed by width and length. This latter pattern was consistent in all three networks. The anisotropy is the second most important factor contributing 25% and crack length 12%. Physical properties make a 20% contribution. Strain energy release rate was determined experimentally through compliance measured for cracked beams of four size categories. The compliance/crack length relationship was found to be non-linear and of an exponential form. Though a good nonlinear regression could be developed for this relationship in the case of an individual specimen, the highly inhomogeneous nature of wood precluded the development of a reliable general regression curve for each geometric category. Consequently, a compliance/crack length relationship was developed for each individual specimen. The strain energy release rate, G Ic was computed for each point of crack propagation in the specimens. The values for G Ic had a large spread; however, the logical methods followed in the study clearly demonstrated that there is more certainty of behaviour in the case of individual specimens but this is masked by the heterogeneity among specimens. Paired Sample t-tests on G Ic for initial and subsequent lengths indicated that the initial crack, which was artificially made, may be a good representation of the naturally occurring crack front in smaller specimens but not necessarily so as specimens gets larger. The average value of fracture toughness over all specimens for initial and all crack lengths was found to be 211 and 218 kPa.m 0.5 , respectively, and these agree with what is expected for this material.

Neural Networks for Predicting Peak Load, Crack Extension and Fracture Toughness and Assessing Contribution of Variables
Three Neural Networks were successfully developed to predict the nature of crack propagation in New Zealand-grown Pinus radiata. When conducting tests with validation data, an R 2 value of 0.865 indicated that the force required to propagate a crack was predicted accurately. Likewise, with an R 2 value of 0.870, the length to which a given crack would extend could be predicted. The fracture toughness, which was obtained through the theory of Linear Elastic Fracture Mechanics, was predicted with an R 2 value of 0.621. This is considered to be highly significant because as discussed earlier, an average value for fracture toughness is normally used; this would by definition have an R 2 of 0.0. Hence, it can be concluded that an ANN, which has been properly developed to predict fracture toughness, can give a better prediction of fracture toughness for an individual piece of wood than can a method that assumes an average value across all specimens. A sensitivity analysis was run for all the variables in each network to determine their relative contribution to the output parameter of that network. Crack length and grain angle predominate the force required for crack propagation, whereas, density and previous crack length along with force predominate crack extension. The latter implies that the length of the initial crack is relevant as was also indicated in Fig. 6. Fracture toughness was mostly determined by dimensional parameters of the specimen implying a volume effect. As for the volume effect on fracture, height is the most prominent followed by width and length in all three networks. Anisotropy (ring and grain angle) was important for all predictions.
These findings are of significance to applications in the timber and building industries. The large spread of fracture strength values, due to the highly inhomogeneous nature of wood, precludes the use of an average value unless an appropriately large safety factor is also adopted; the result is that compensating for the ensuing uncertainty often involves excess material usage (over-sizing). A method of grading timber that used neural networks, or gave consideration to the effect of individual member characteristics, would predict fracture strength for individual members and may enable efficient use of wood in future.