## 2 Materials and Methods
### Numerical Model  
A structural-mechanical model of the four-point bending test was simulated using the commercial software LS-DYNA. The 8 layers of the woven fabric composite were modeled with Belytschko-Tsai 4-node shell elements, with one integration point through the thickness. For each layer of the composite, a layer of shell elements was considered. Solid elements of cohesive material were placed between the shell layers, in order to represent the resin among layers. Solid elements of cohesive consider a specific formulation of LS-DYNA, which allows the transfer of moments to the shells. In particular, nodes of solid elements coincide with those of the shells. The material card \"MAT_LAMINATED_COMPOSITE_FABRIC, which is specific for fabric composite, was adopted for the composite layers. This material model assumes a slightly modified version of Hashin criteria for failure [18]. However, in this paper, the focus is on the elastic field, and parameters related to the post-failure behavior of the material will not be addressed. The cohesive material model \"MAT_COHESIVE_TH, with properties of the resin, was assigned to the cohesive elements. Finally, the span noses were modeled as rigid cylinder walls. The upper cylinders, shown in Figure 1b, move downwards with a prescribed motion law. The upper and lower cylinders were placed at the initial distance of 1 mm from the laminate in order to avoid penetrations and numerical instabilities. Contact between laminate and cylinders was based on a penalty formulation with a penalty factor which adaptively guarantees the numerical stability [24]. This formulation accounts for thickness offsets of shells. The accuracy of the numerical model was firstly verified with respect to the four point-bending tests of the two composite plates. Then, the same numerical model was used in comparison with the experimental results of the progressively damaged plate. In this case, the structural-mechanical model had to account for the variation of Young's modulus on the plate, as pointed out through the \(DI_{d}\) technique. The residual elastic properties of the damaged plate are mapped through the \(DI_{d}\) technique. The Young's modulus varies along the x-direction according to the location of measurement, with the damage mostly located in correspondence of internal rolls. It seems reasonable to assume that the Young's modulus varies continuously along the longitudinal direction. Here, a polynomial curve is considered for each region included between two consecutive measure locations. For each polynomial curve, the two Young's moduli at the extremities, \(E_{k}\) and \(E_{k+1}\), with k the considered location, are known. In addition, the two derivatives at the extremities can be estimated. In order to evaluate the derivates, the difference quotient for each region has been calculated as  
\[\frac{\Delta E_{k}}{\Delta x_{k}}=\frac{E_{k}-E_{k-1}}{x_{k}-x_{k-1}} \tag{1}\]  
in accordance with its definition. The difference quotient is a measure of the average rate of change of the function, here the Young's modulus, over the interval \(\Delta x_{k}\). By considering two consecutive difference quotients, \(\frac{\Delta E_{k}}{\Delta x_{k}}\) and \(\frac{\Delta E_{k+1}}{\Delta x_{k+1}}\), the derivative at the k-th location has been estimated as follows:  
1. when the product \(\frac{\Delta E_{k}}{\Delta x_{k}}\frac{\Delta E_{k+1}}{\Delta x_{k+1}}\) returns a negative value, the derivate at the k-th location is assumed equal to zero. This can be justified by considering that a change in the sign of the difference quotient implies a change in the derivative of the function, as well;
2. when the sign of the product \(\frac{\Delta E_{k}}{\Delta x_{k}}\frac{\Delta E_{k+1}}{\Delta x_{k+1}}\) is positive, the derivate at the k-th location is assumed equal to average value \(\frac{\frac{\Delta E_{k}}{\Delta x_{k}}+\frac{\Delta E_{k+1}}{\Delta x_{k+1}}}{2}\);
3. the derivates at the extremities, \(k=1\) and \(k=7\), are assumed equal to zero.