## 2 Materials and Methods
### Numerical Analysis
After some algebraic manipulation involving Equations (1)-(5), the relative displacement at damage onset becomes  
\[\delta_{0,\text{m}}=\frac{\delta_{0,\text{I}}\delta_{0,\text{II}}\sqrt{1+\beta ^{2}}}{\sqrt{\beta^{2}\delta_{0,\text{I}}^{2}+\delta_{0,\text{II}}^{2}}} \tag{6}\]  
Figure 4: Bilinear cohesive law in mixed\(-\)mode I + II implemented in ABAQUS. Separation (\(\delta\)) and traction (\(f\)). where \(\delta_{0,i}\) (\(i\) = I, II) are the relative displacements at damage onset under pure mode loading (Figure 4). During propagation, the traction-separation relationship can be directly defined as  
\[f_{\rm m}=(1-d_{\rm m})k\delta_{\rm m} \tag{7}\]  
where \(d_{\rm m}\) is the damage parameter defined from the linear softening relationship between \(\delta_{0,\rm m}\) and the ultimate relative displacement \(\delta_{\rm u,m}\) according to  
\[d_{\rm m}=\frac{\delta_{\rm u,m}(\delta_{\rm m}-\delta_{0,\rm m})}{\delta_{\rm m }(\delta_{\rm u,m}-\delta_{0,\rm m})} \tag{8}\]  
The ultimate relative displacement \(\delta_{\rm u,m}\) is calculated by means of the linear energetic criterion used to simulate damage propagation under mixed-mode I + II loading:  
\[\delta_{\text{m}}=\sqrt{\delta_{1}^{2}+\delta_{\text{II}}^{2}} \tag{1}\]0  
where \(G_{\rm ic}\) (\(i\) = I, II) are the critical fracture energy and \(G_{i}\) (\(i\) = I, II) are the energy components of the current mode mixity given by  
\[\delta_{\text{m}}=\sqrt{\delta_{1}^{2}+\delta_{\text{II}}^{2}} \tag{1}\]1  
with \(\delta_{0,\rm m,i}\) and \(\delta_{\rm u,m,i}\) representing the critical mode \(i\) components of the relative displacements at damage onset and at failure, respectively. From Equations (3) and (10), the strain energy release rate under mixed-mode loading can be written as  
\[\delta_{\text{m}}=\sqrt{\delta_{1}^{2}+\delta_{\text{II}}^{2}} \tag{1}\]2  
Combining Equations (11) and (9) yields  
\[\delta_{\text{m}}=\sqrt{\delta_{1}^{2}+\delta_{\text{II}}^{2}} \tag{1}\]3  
Knowing the critical strain energy release rate (\(G_{\rm Tc}\)), \(\delta_{\rm u,m}\) can be determined as follows:  
\[\delta_{\text{m}}=\sqrt{\delta_{1}^{2}+\delta_{\text{II}}^{2}} \tag{1}\]4  
which makes it possible to define the damage parameter (Equation (8)) that represents the progressive propagation of damage under mixed-mode I + II loading. It should be noted that this general formulation also addresses the pure mode loading cases (I or II) since they are particular cases of the mixed-mode loading scenario. The mean values of the fracture parameters of European beech defining the mode I and II cohesive laws in the TL crack propagation system, as well as the elastic constants of the material, were obtained from previous work by the authors [10; 17; 18] and are compiled in Tables 2 and 3. For the dowel steel, a modulus of elasticity of 210,000 N/mm\({}^{2}\) and a Poisson's ratio of 0.3 were considered. These material properties were used as input data in the numerical models. \begin{table}
\begin{tabular}{c c c c c c c c c c} \hline  & \(E_{\rm L}\) & \(E_{\rm R}\) & \(E_{\rm T}\) & \(\
u_{\rm RL}\) & \(\
u_{\rm TL}\) & \(\
u_{\rm RT}\) & \(G_{\rm LR}\) & \(G_{\rm LT}\) & \(G_{\rm RT}\) \\  & **(N/mm\({}^{2}\))** & **(N/mm\({}^{2}\))** & **(N/mm\({}^{2}\))** & **(-)** & **(-)** & **(-)** & **(N/mm\({}^{2}\))** & **(N/mm\({}^{2}\))** \\ \hline Average & 13,811 & 1590 & 832 & 0.51 & 0.44 & 0.32 & 1108 & 706 & 349 \\ St. Dev. & 1323 & 541 & 115 & 0.030 & 0.015 & 0.041 & 202 & 139 & 53 \\ \hline \end{tabular}
\end{table} 
Table 2: Average values and standard deviation of the elastic constants of European beech [10].