Experimental and Numerical Analysis of Multiple Low-Velocity Impact Damages in a Glass Fibered Composite Structure

Glass fiber-reinforced polymer structures (GFRPS) are widely used in civil and mechanical fields due to their light weight and corrosion resistance. However, these structures are prone to damage with very-low-energy impacts. The reliability of such structures is of prime importance before their installation and usage. This study aimed to identify, visualize, localize, and verify multiple barely visible impact damage (BVID) in a GFRPS using a combination of guided waves (GW)-based online structural health monitoring (SHM) and thermal strain-based nondestructive testing (NDT) approaches. Global NDT techniques like the use of a laser Doppler vibrometer (LDV) and digital image correlation (DIC) were used in the experimental analysis. The effectiveness of the experimental LDV-GW process was also checked numerically with the spectral element method (SEM). A threshold-based baseline free SHM approach to effectively localize the damages was proposed along with quick DIC verification of composite structure with thermal loading based on short-pulse heating as an excitation source. This study analyzed combined experimental- and numerical-based SHM-NDT methods in characterizing the multiple BVIDs located in a GFRPS.


Introduction
Composites are light weight, have high strength, and are used in many automotive, civil, aerospace, etc. applications. The GFRPS laminates have greater strength along the direction of fiber [1] GFPRS are widely used in retrofitting structures in construction industries, and, nowadays, these glass fibers are combined to prepare reinforcement bars for structural strengthening [2]. They are highly preferred due to their high tensile strength, low density, and corrosion resistance [3]. Even though several advantages exist, they also have disadvantages like delamination defects, impact cracks [4], internal fiber matrix damage [5], etc. If such damage fails to be noticed, it could lead to severe structural disintegration [6]. BVID and impact cracks in the GFRPS are of most common that occurs even due to a low-energy impact force and are mostly not visible for internal visual identification.
The engineering industry offers many methods for analyzing the defects in the structures. The NDT analysis of glass fibers using a combination of X-ray tomography, ultrasonic imaging, and optical proliferometry was proposed by Miqoi et al. [7] and on carbon fibers using digital image correlation (DIC) and scanning electron microscopy by Sarasini et al. [8], in addition to experimental compression after impact tests by Talreja et al. [9]. SHM finds its way in detecting such impact damage by performing continuous health monitoring of structures 1. SHM reduces the laborious human hours of checking the whole structure for baseline free elliptical threshold method (ETM) was proposed to show the effective online localization of the damages. Thus the proposed SHM-GW technique can be used as the first step of the quick analysis to locate the damages followed by NDT-DIC verification based on SHM results. The SHM combined with the NDT method can drastically reduce the time to analyze larger and complex structures.

Methodological Flow Chart
The flowchart (Figure 1) briefly explains the methodological process involved in analyzing the GFRPS using combined NDT and SHM methods. The process involves damage-identification, visualization, localization, and verification. A numerical method based on SEM code was developed and compared with the experimental LDV signals to analyze the damage. Since all three damages existed before, the numerical simulation was utilized to study each damage scenario separately. The RRMS method was applied to the signal data obtained from both LDV and SEM. The SHM-based baseline free elliptical threshold method (ETM) was proposed to show the effective online localization of the damages. Thus the proposed SHM-GW technique can be used as the first step of the quick analysis to locate the damages followed by NDT-DIC verification based on SHM results. The SHM combined with the NDT method can drastically reduce the time to analyze larger and complex structures.

Methodological Flow Chart
The flowchart (Figure 1) briefly explains the methodological process involved in analyzing the GFRPS using combined NDT and SHM methods. The process involves damage-identification, visualization, localization, and verification.

GFRP Material Studied
The 0.2 cm thick GFRPS has 12 layups stacked in [0/90]3s format, which are bonded together by HEXCEL-212Na adhesive. The material is procured from the G.ANGELONI group [30] (Venice, Italy). These materials are largely used in the mechanical and construction industries. The GFRPS is of dimensions 50 × 50 cm 2 in length and breadth. Multiple BVIDs of Ø 0.9 cm (BVID1) and Ø 1.2 cm (BVID2 and BVID3) are made with an impact force of 20, 25, and 30J. This is done by dropping a steel ball attached to a steel beam from varying heights guided via a pipe. All the 4 PZTs (Ø 1 cm and 0.1cm thin) are

GFRP Material Studied
The 0.2 cm thick GFRPS has 12 layups stacked in [0/90] 3s format, which are bonded together by HEXCEL-212Na adhesive. The material is procured from the G.ANGELONI group [30] (Venice, Italy). These materials are largely used in the mechanical and construction industries. The GFRPS is of dimensions 50 × 50 cm 2 in length and breadth. Multiple BVIDs of Ø 0.9 cm (BVID1) and Ø 1.2 cm (BVID2 and BVID3) are made with an impact force of 20, 25, and 30J. This is done by dropping a steel ball attached to a steel beam from varying heights guided via a pipe. All the 4 PZTs (Ø 1 cm and 0.1 cm thin) are made from SONOX material P502 [31] and are attached to the GFRPS using cyanoacrylate glue. The PZTs (S1-S4) and BVIDs location coordinates are presented in Figure 2. made from SONOX material P502 [31] and are attached to the GFRPS using cyanoacrylate glue. The PZTs (S1-S4) and BVIDs location coordinates are presented in Figure 2.

LDV Setup
A single head (1D) LDV was chosen, which analyzed the structures using out-ofplane wavefield components. Polytec PSV 400 was the LDV type used. LDV helps to visualize the GW propagation and also to acquire the GW full-wave fields. The GW propagation was done by exciting the PZT. The retro-reflective film was applied only to an area of 26 cm × 50 cm (blue dotted lines), as shown in Figure 2, to reduce the LDV calculation time. The LDV setup consisted of the signal generator, the power amplifier, and the oscilloscope. A 5-cycle sinusoidal tone burst generated with a Hanning window was applied to the PZT. Frequencies of 50 kHz, 100 kHz, and 200 kHz were used in this analysis. The LDV measurement at each grid point was time-averaged 5 times with an overall sampling period of 1024. A 16 Vpp voltage with a gain of 20 was applied to PZT via the generator. All four symmetrical positioned PZTs were excited via this process.

SEM Numerical Approach
The simulation was conducted with the algorithm developed in Matlab software (version 2020a) based on the implementation found in [32]. The GFRP structure was modelled according to the Mindlin-Reissner first-order shear deformation theory [33], and the laminate theory [34] was used to determine the effective engineering constants shown in Table 1. The details of the determination of the effective constants are given in Appendix A. The model of PZT transducers was omitted in the analysis to simplify calculations. Instead, the out-of-plane forces were applied to the nodes in the place of the actuator. The damage was modulated in two ways. First, Young's modulus of the structure was reduced by 50% in the circular of the damaged area. Second, three cracks were introduced using node separation [35] at the location marked by the dashed line in Figure 3. The LDV experimental setup values, as mentioned earlier, were used in studying the numerical GFRP models.

LDV Setup
A single head (1D) LDV was chosen, which analyzed the structures using out-of-plane wavefield components. Polytec PSV 400 was the LDV type used. LDV helps to visualize the GW propagation and also to acquire the GW full-wave fields. The GW propagation was done by exciting the PZT. The retro-reflective film was applied only to an area of 26 cm × 50 cm (blue dotted lines), as shown in Figure 2, to reduce the LDV calculation time. The LDV setup consisted of the signal generator, the power amplifier, and the oscilloscope. A 5-cycle sinusoidal tone burst generated with a Hanning window was applied to the PZT. Frequencies of 50 kHz, 100 kHz, and 200 kHz were used in this analysis. The LDV measurement at each grid point was time-averaged 5 times with an overall sampling period of 1024. A 16 V pp voltage with a gain of 20 was applied to PZT via the generator. All four symmetrical positioned PZTs were excited via this process.

SEM Numerical Approach
The simulation was conducted with the algorithm developed in Matlab software (version 2020a) based on the implementation found in [32]. The GFRP structure was modelled according to the Mindlin-Reissner first-order shear deformation theory [33], and the laminate theory [34] was used to determine the effective engineering constants shown in Table 1. The details of the determination of the effective constants are given in Appendix A. The model of PZT transducers was omitted in the analysis to simplify calculations. Instead, the out-of-plane forces were applied to the nodes in the place of the actuator. The damage was modulated in two ways. First, Young's modulus of the structure was reduced by 50% in the circular of the damaged area. Second, three cracks were introduced using node separation [35] at the location marked by the dashed line in Figure 3. The LDV experimental setup values, as mentioned earlier, were used in studying the numerical GFRP models. where E represents the Young's modulus, G represents the shear modulus,  represents the Poisson's ratio,  is the density, and V is the volume fraction. The four model cases ( Figure 3) were made numerically as the structure obtained already had three BVIDs in the initial state itself. Thus, SEM models added an advantage to visualize how the GW propagation would be with the individual BVIDs.

DIC Setup
As a verification procedure, 3D DIC was proposed; so, the measuring setup consisted of two USB3.0 Baumer 12.3 Mpx cameras (Baumer GmbH, Radeberg, Germany) with VS-1220HV lenses, combined in a Q400 system by DANTEC DYNAMICS GmbH. The strain field changes resulted from heating the GFRP-plate using a halogen lamp HEDLER HF 65 (HEDLER GmbH, Hesse, Germany); the exposition time was 2 min. The size of the observed area was a compromise between the highest possible accuracy and the widest observed area. The area observed by each camera was about 340 mm in width and 249 mm high, which considering the camera resolution of 4096 × 3000 px gave the spatial resolution equal to 83 μm/px. The calibration was made using the Al-15-BMB_9×9 calibration target. During the evaluation of pictures, the facet size was established to 17 pixels. In the presented examination, the mean accuracy of the displacement measurement was: in the horizontal direction, 0.00006 mm; in the vertical direction, 0.00007 mm; and, perpendicular to the observed plane, 0.002 mm. The whole specimen was not observed. The measuring stand is shown in Figure 4a, and the scheme of the experiment is shown in Figure 4b, top view. The four model cases ( Figure 3) were made numerically as the structure obtained already had three BVIDs in the initial state itself. Thus, SEM models added an advantage to visualize how the GW propagation would be with the individual BVIDs.

DIC Setup
As a verification procedure, 3D DIC was proposed; so, the measuring setup consisted of two USB3.0 Baumer 12.3 Mpx cameras (Baumer GmbH, Radeberg, Germany) with VS-1220HV lenses, combined in a Q400 system by DANTEC DYNAMICS GmbH. The strain field changes resulted from heating the GFRP-plate using a halogen lamp HEDLER HF 65 (HEDLER GmbH, Hesse, Germany); the exposition time was 2 min. The size of the observed area was a compromise between the highest possible accuracy and the widest observed area. The area observed by each camera was about 340 mm in width and 249 mm high, which considering the camera resolution of 4096 × 3000 px gave the spatial resolution equal to 83 µm/px. The calibration was made using the Al-15-BMB_9×9 calibration target. During the evaluation of pictures, the facet size was established to 17 pixels. In the presented examination, the mean accuracy of the displacement measurement was: in the horizontal direction, 0.00006 mm; in the vertical direction, 0.00007 mm; and, perpendicular to the observed plane, 0.002 mm. The whole specimen was not observed. The measuring stand is shown in Figure 4a, and the scheme of the experiment is shown in Figure 4b, top view.
target. During the evaluation of pictures, the facet size was established to 17 pixels. In the presented examination, the mean accuracy of the displacement measurement was: in the horizontal direction, 0.00006 mm; in the vertical direction, 0.00007 mm; and, perpendicular to the observed plane, 0.002 mm. The whole specimen was not observed. The measuring stand is shown in Figure 4a, and the scheme of the experiment is shown in Figure 4b

Damage Analysis Methods
The damage identification and visualization were done using RRMS for both LDV and SEM data (Section 3.1), localized by the developed SHM-based ETM method (Section 3.2) and verified by the proposed DIC thermal speckle-based variation method (Section 3.3).

RRMS-Based Damage Analysis
The visualization method using RRMS [18] for four PZTs excitations were analyzed. The RMS provides the energy distribution of the signals. However, if more energy gets accumulated in a certain region, it shows a higher value of RMS. It results in blurring of the damage zone sometimes. The RRMS function Equation (1) was added to counteract such cases, which enhanced the results obtained.
where k i -ith sample amplitude, and N-the number of registered time samples. RWF is the radially weighted factored value for each scanning point, and "a" is the power of the weight factor. In this function, the excitation points were the origin point (x 0 , y 0 ), and the radial distance was measured for all the points (x i , y in ), as shown in Equation (2).
With the help of the obtained Equation (2), a combination of LDV and SEM models were analyzed as mentioned in Table 2.

Threshold-Based Elliptical Method in Damage Analysis
The localization of the damages was done using the SHM-based strategy using the ETM. It works on the methodology of convergence and interaction of ellipses. The method requires the time of arrival (TOA) of the wave peaks in the damage calculation [36,37]. The time of arrival (theoretical) (TOA T ) was obtained between GW paths to the arbitrary grid points, the sensor, and the actuator. After applying the Hilbert transformation (HT) envelope to the signals, amplitude peaks were obtained. The obtained TOA T was then tried to match up with the HT peak values (TOA E ).
A higher value in the mesh grid was obtained when the TOA T matched with an HT value. A fine mesh of 0.1 mm spacing was used in the calculation. Since most reflections occurred at the damage spot or near the damage spot, that region was highlighted in the grid. A threshold value (90%) was applied to the algorithm to separate the higher values in the grid during the run. The threshold value was selected based on the experience and the references from the threshold-based probabilistic study [38]. The damage localization approach is explained in a flowchart, as shown in Figure 5. (G1, G2) are arbitrary grid point coordinates; (A1, A2, S1, S2) are actuator and sensor coordinates; VA, VS represents the group velocities obtained from the actuator to damage and the sensor to damage, respectively; A T represents the amplification factor; TOA E represents the time of arrival obtained from the experiment; D (x, y) represents the damage index; AS represents the total number of actuator sensor paths studied; τ is the decay factor value of 5 μs, and N represents threshold value used in the calculation.

Thermal Speckle-Based Verification Method
The examined GFRPS was coated by the typical DIC speckle pattern (registered by two cameras as shown in Figure 6), so the damages were invisible. The specimen was observed under the heating process realized by the light source. Since the damage location was initially localized with SHM, the lamps were focused in the localized area. The reference picture was made at stable room temperature, and then the whole heating took two minutes. During this time, each camera registered two new photos every one second (sampling frequency equal to 2 Hz). For the analysis of picture sequences, the commercial ISTRA 4D software (version 4.6) was used.

RRMS Models-Visualization of Damages
RRMS studies were carried out with the LDV and SEM data to visualize the damage. A constant power/weight factor value of a =0.5 was used for all the calculations.  (G 1, G 2 ) are arbitrary grid point coordinates; (A 1 , A 2, S 1, S 2 ) are actuator and sensor coordinates; V A , V S represents the group velocities obtained from the actuator to damage and the sensor to damage, respectively; A T represents the amplification factor; TOA E represents the time of arrival obtained from the experiment; D (x, y) represents the damage index; AS represents the total number of actuator sensor paths studied; τ is the decay factor value of 5 µs, and N represents threshold value used in the calculation.

Thermal Speckle-Based Verification Method
The examined GFRPS was coated by the typical DIC speckle pattern (registered by two cameras as shown in Figure 6), so the damages were invisible. The specimen was observed under the heating process realized by the light source. Since the damage location was initially localized with SHM, the lamps were focused in the localized area. The reference picture was made at stable room temperature, and then the whole heating took two minutes. During this time, each camera registered two new photos every one second (sampling frequency equal to 2 Hz). For the analysis of picture sequences, the commercial ISTRA 4D software (version 4.6) was used.
14, x FOR PEER REVIEW 7 of 14 (G1, G2) are arbitrary grid point coordinates; (A1, A2, S1, S2) are actuator and sensor coordinates; VA, VS represents the group velocities obtained from the actuator to damage and the sensor to damage, respectively; A T represents the amplification factor; TOA E represents the time of arrival obtained from the experiment; D (x, y) represents the damage index; AS represents the total number of actuator sensor paths studied; τ is the decay factor value of 5 μs, and N represents threshold value used in the calculation.

Thermal Speckle-Based Verification Method
The examined GFRPS was coated by the typical DIC speckle pattern (registered by two cameras as shown in Figure 6), so the damages were invisible. The specimen was observed under the heating process realized by the light source. Since the damage location was initially localized with SHM, the lamps were focused in the localized area. The reference picture was made at stable room temperature, and then the whole heating took two minutes. During this time, each camera registered two new photos every one second (sampling frequency equal to 2 Hz). For the analysis of picture sequences, the commercial ISTRA 4D software (version 4.6) was used.

RRMS Models-Visualization of Damages
RRMS studies were carried out with the LDV and SEM data to visualize the damage. A constant power/weight factor value of a =0.5 was used for all the calculations. Table 2's

RRMS Models-Visualization of Damages
RRMS studies were carried out with the LDV and SEM data to visualize the damage. A constant power/weight factor value of a = 0.5 was used for all the calculations. Table 2's case results are shown in Figures 7 and 8. Initially, cases 4 and 5 ( Figure 7) were considered for the visualization study as they represent the real damage scenarios. After analyzing all the frequencies, it was found that both LDV_RRMS and SEM_RRMS 200 kHz plots identified all the BVIDs more prominently with rounded red circles (Figure 7). Table 2  After checking the results from S1, similar BVID identification and visualization checks were done for other symmetrically located sensors, as shown in Figure 9, using the RRMS-based weightage formulation. All sensors identified the damages clearly both from the experimental and the numerical data.  After analyzing all the frequencies, it was found that both LDV_RRMS and SEM_RRMS 200 kHz plots identified all the BVIDs more prominently with rounded red circles (Figure 7). Table 2's other cases (1)(2)(3) were also checked only with 200 kHz. Thus, only 200 kHz was chosen as the identified frequency for other RRMS and SHM results (Section 4.2). The SEM_RRMS plots, as shown in Figure 8, identified the inner numerically made BVIDs (1, 2, and 3) separately and proved that the methodology can be used for single damages. After checking the results from S1, similar BVID identification and visualization checks were done for other symmetrically located sensors, as shown in Figure 9, using the RRMS-based weightage formulation. All sensors identified the damages clearly both from the experimental and the numerical data. After analyzing all the frequencies, it was found that both LDV_RRMS and SEM_RRMS 200 kHz plots identified all the BVIDs more prominently with rounded red circles (Figure 7). Table 2's other cases (1)(2)(3) were also checked only with 200 kHz. Thus, only 200 kHz was chosen as the identified frequency for other RRMS and SHM results (Section 4.2). The SEM_RRMS plots, as shown in Figure 8, identified the inner numerically made BVIDs (1, 2, and 3) separately and proved that the methodology can be used for single damages.
After checking the results from S1, similar BVID identification and visualization checks were done for other symmetrically located sensors, as shown in Figure 9, using the RRMS-based weightage formulation. All sensors identified the damages clearly both from the experimental and the numerical data.

ETM-Localization of Damages
The amplification and threshold factor-added elliptical-based algorithm was applied to Table 2 in damage localization. The damage points were identified based on elliptical intersections. Equidistant nodal points (P1-P10) were taken from the LDV area scan results with excitation from S1 (Figure 10a) for the process. Since S2 and S3 are in a straight line to the detected damages, they were not considered for taking the SHM study. After analyzing S1 and S4, the S1 RRMS maps detected the damage better than its symmetrical counterpart S4. A frequency range of 200 kHz was considered for this study, as mentioned earlier. The full wavefield data of the cases showed that the A0 GW mode identified the damage cases clearly, which was chosen to study further. The A0 mode velocity values  After checking the results from S1, similar BVID identification and visualization checks were done for other symmetrically located sensors, as shown in Figure 9, using the RRMS-based weightage formulation. All sensors identified the damages clearly both from the experimental and the numerical data.

ETM-Localization of Damages
The amplification and threshold factor-added elliptical-based algorithm was applied to Table 2 in damage localization. The damage points were identified based on elliptical intersections. Equidistant nodal points (P1-P10) were taken from the LDV area scan results with excitation from S1 (Figure 10a) for the process. Since S2 and S3 are in a straight line to the detected damages, they were not considered for taking the SHM study. After analyzing S1 and S4, the S1 RRMS maps detected the damage better than its symmetrical counterpart S4. A frequency range of 200 kHz was considered for this study, as mentioned earlier. The full wavefield data of the cases showed that the A0 GW mode identified the damage cases clearly, which was chosen to study further. The A0 mode velocity values obtained for the LDV by scanning points at 90° and 0° were 2231 m/s, 1768 m/s, 2135 m/s, and 1677 m/s for the SEM calculations. The velocity profile was obtained by fitting the velocity values to an elliptical function and determining the 360° velocity profile.
The damage localization plots are shown in Figure 10b-f. The maps were plotted for all five cases (4 SEM and 1 LDV), and higher-energy regions were isolated with threshold values. The region shows that higher energy was nearer to the created BVIDs (circled in black). The error estimation is shown in Tables 3 and 4. The values obtained were closer to the BVID regions and are in an error marginal range of >1 cm difference maximum. The damage localization plots are shown in Figure 10b-f. The maps were plotted for all five cases (4 SEM and 1 LDV), and higher-energy regions were isolated with threshold values. The region shows that higher energy was nearer to the created BVIDs (circled in black). The error estimation is shown in Tables 3 and 4. The values obtained were closer to the BVID regions and are in an error marginal range of >1 cm difference maximum. velocity values to an elliptical function and determining the 360° velocity profile. The damage localization plots are shown in Figure 10b-f. The maps were plotted for all five cases (4 SEM and 1 LDV), and higher-energy regions were isolated with threshold values. The region shows that higher energy was nearer to the created BVIDs (circled in black). The error estimation is shown in Tables 3 and 4. The values obtained were closer to the BVID regions and are in an error marginal range of >1 cm difference maximum.

Thermal Speckle-Based Verification of Damages
The collected data were carefully analyzed, considering the displacements in three directions and strain fields. In Figure 11a-c, the horizontal displacement of points in the first, 120th, and 240th steps (after two minutes of heating) were shown, respectively. In the figures, the horizontal displacement field was quasi even in each step and corresponds to the boundary conditions (clamping on the right side). In such

Thermal Speckle-Based Verification of Damages
The collected data were carefully analyzed, considering the displacements in three directions and strain fields. In Figure 11a-c, the horizontal displacement of points in the first, 120th, and 240th steps (after two minutes of heating) were shown, respectively. In the figures, the horizontal displacement field was quasi even in each step and corresponds to the boundary conditions (clamping on the right side). In such a short period of heating, there were no visible differences between point displacements in damaged and undamaged regions.

Conclusions
In this study, a combined SHM-GW and NDT-DIC method was proposed in studying BVIDs presented in the GFRPS. A frequency range of 200 kHz was analyzed for the GWbased studies. SEM numerical models were compared with the experimental data.
▪ SEM helped to model BVID separately and together (experimental data) to validate the SHM methodology. ▪ The RRMS study provided results about the damage locations even near the higherenergy zones. ▪ S1 provided good analyses of the results, and scanning points were taken to analyze the damage zones. ▪ The threshold-based proposed ETM algorithm predicted the location of the damages. ▪ The error range was less than 1%, as shown in the difference (cm) of Table 4 in all the analyzed cases. ▪ DIC-based diagrams of variation in displacements differentiated the damage region and thus verified it. ▪ Combined online monitoring and verification of the structure was proposed.
Therefore, it is promising to use the proposed SHM-GW and DIC-NDT application in monitoring mechanical, civil structures to verify structural integrity. Future research studies involve implementing similar threshold-based localization and verification approaches in analyzing different structures for certification and reliability purposes.  However, the estimated standard deviation of the displacement in the horizontal direction for every data point (shown in Figure 11d) allows to differentiate the failure areas (damages); the standard deviation for damages was greater than for healthy regions in the neighbourhood. Similar observations were made for the other directions.

Conclusions
In this study, a combined SHM-GW and NDT-DIC method was proposed in studying BVIDs presented in the GFRPS. A frequency range of 200 kHz was analyzed for the GW-based studies. SEM numerical models were compared with the experimental data.
I SEM helped to model BVID separately and together (experimental data) to validate the SHM methodology.

I
The RRMS study provided results about the damage locations even near the higherenergy zones. I S1 provided good analyses of the results, and scanning points were taken to analyze the damage zones.

I
The threshold-based proposed ETM algorithm predicted the location of the damages.

I
The error range was less than 1%, as shown in the difference (cm) of Table 4 in all the analyzed cases. Therefore, it is promising to use the proposed SHM-GW and DIC-NDT application in monitoring mechanical, civil structures to verify structural integrity. Future research studies involve implementing similar threshold-based localization and verification approaches in analyzing different structures for certification and reliability purposes.

Data Availability Statement:
The data presented in this study are available on request from the first, corresponding authors.

Acknowledgments:
The authors also acknowledge Task-CI for the use of computational resources. This publication is based upon work from COST Action CA18123 (ODIN), supported by COST (European Cooperation in Science and Technology-Brussels, Belgium).

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The elastic constants of the unidirectional fibre composite in the form of the Hahn functional are gathered in Table A1. The Hahn function is given in the form: where V f , V m are the volume fraction of the fibers and matrix, respectively. The bulk modulus (K T ) for fibers and the matrix is equal to K f = E f /2 1 − υ f , K m = E m /2(1 − υ m ), respectively, and the expressions for η are given as follow: The transverse moduli and υ 23 of the composite are given: . According to the rules, the mixture mass density of the composite is equal ρ = ρ f V f + ρ m V m .