Numerical Modeling of Damage Caused by Seawater Exposure on Mechanical Strength in Fiber-Reinforced Polymer Composites

Fiber-reinforced polymer composites are frequently used in marine environments which may limit their durability. The development of accurate engineering tools capable of simulating the effect of seawater on material strength can improve design and reduce structural costs. This paper presents a numerical-based approach to predict the stress–strain response of fiber-reinforced polymer composites exposed to different seawater immersion times, ranging from 0 to 900 days. A three-dimensional numerical model has been implemented using a static implicit finite element analysis along with a user-defined material (UMAT) subroutine. Puck’s failure criterion was used for ultimate failure analysis of the laminates, while Fick’s first diffusion law was used to predict the seawater absorption rate. Overall, the simulated stress–strain curves were close to those obtained experimentally. Moreover, the model agreed well with the experimental data regarding the maximum stress and the strain at failure leading to maximum errors lower than 9% and 11%, respectively. Additionally, the simulated strain fields agreed well with the experimental results measured by digital image correlation. Finally, the proposed procedure was also used to identify the most critical surfaces to protect the mechanical components from marine environments.


Introduction
Fiber-reinforced polymers (FRP) possess very attractive physical and mechanical properties, such as good fatigue resistance, high specific stiffness, high specific strength, and easy construction relative to traditional materials [1]. This is reflected in the growth of their use in areas where light and strong structures are required, for example, in the aerospace and automobile industries. The commonly used FRP composite reinforcements are made from carbon fiber, glass fiber, basalt fiber and aramid fiber [2]. Glass fiberreinforced polymers (GFRP) are the most frequently used FRPs mainly due to their low cost and relatively good tensile strength [3]. On the other hand, basalt fiber reinforced polymers (BFRP) are more expensive but offer improved specifications such as higher strength and chemical resistance [4]. In the case of aramid fiber reinforced polymers (AFRP), they are not frequently used in structural applications due to low compressive strength and high cost. However, they are one of the best choices for ballistic protection applications [5]. Lastly, carbon fibers possess exceptional mechanical properties having high specific strength, low density, and also good electrical properties, namely excellent conductivity [6]. As a result, although carbon fiber-reinforced polymers (CFRP) are high-priced, they are also widely used in applications requiring high specific stiffness, enhanced fatigue properties and high specific strength, including aerospace, military and wind power applications [7,8]. In addition, due to their good corrosion resistance, FRPs have aroused great interest in multiple marine applications, namely in ship hulls, sonar domes, propellers, railings, submarine structures and offshore sectors [9][10][11][12]. However, although FRPs have been used for over 50 years in marine structures, their adoption frequently requires high design safety factors, since long-term durability predictions are generally based on experience [13]. This can be partly justified by the fact that the usage of composite materials in hostile environments can significantly reduce their mechanical behaviour.
Extensive research has been conducted to evaluate the effects of hostile environments on composite materials [14][15][16][17][18][19][20]. The absorption of water is detrimental to the mechanical properties of FRP composites, since it affects the resin matrix by plasticization and also deteriorates the matrix-fiber interface [21][22][23][24]. However, the degradation of the matrix resin due to water absorption is, in general, less relevant than the degradation of the matrix-fiber interface. In fact, several studies have demonstrated that the bonding between the matrix and fiber is the predominant factor causing the degradation of the composite [21,[25][26][27][28]. Furthermore, Apicella et al. [29] demonstrated that water can also hydrolyze the resin [30]. The same study showed that isophthalic resins revealed the lowest hydrolytic stability, while, on the contrary, vinyl ester resins revealed the highest stability after exposure.
Temperature changes also play an important role in mechanical behaviour of composite materials. Netravali et al. [31] studied the effect of water sorption at different temperatures on the reversibility of the plasticization and concluded that, at ambient temperature, the plasticization was mostly reversible. However, the samples exposed at 70 • C revealed a considerably lower glass transition temperature. In another study, Wong et al. [32] addressed the effect of multiple sorption cycles in a glassy epoxy resin. It was found that anomalous sorption cycles are verified not as a result of the development of cracks but also due to an alteration in resin chemical structure caused by an oxidation reaction. Ellis et al. [33] investigated the influence of water in epoxy resins and concluded that even a small amount of water uptake (e.g., 1-2 wt%) has a considerable plasticizing action, leading to degraded properties.
As described above, the damage induced by seawater can emerge in numerous forms. Analogously, water absorption in composite materials is also a complex phenomenon, which many causes can influence, such as the resin and curing agent, fabrication, fiber volume fraction, orientation of reinforcement, area of exposed surfaces, and temperature, among others [27,34,35]. Water penetrates the structure mainly via diffusion, through the resin matrix [36,37], and most water absorption is absorbed in the matrix. However, Tsenoglou et al. [38] evaluated the water penetration of an imperfectly bounded fibermatrix interface and found that, essentially in low quality composites, the micro-channels surrounding the fibers can be the predominant conduits of moisture diffusion. Furthermore, voids can be developed during the composite cure process, which also increases the water absorption during the FRP life cycle [39]. Weitsman et al. [40] analysed the microcracking influence on the water absorption behaviour in polymers and found that micro-cracks increase the absorption capacity and accelerate the diffusion process. Although microcracking is extensively linked with the increase in absorption capacity, it is still highly random and difficult to evaluate [40]. Pritchard and Speake [41] studied the kinetics of moisture absorption on GFRP for temperatures between 30 and 100 • C and found that the saturation point is reached more quickly at higher temperatures. However, the maximum water absorption remained unaltered. The same study showed that the residual mechanical properties were found to be functions of true absorbed water content.
Saadatmanesh et al. [42] studied the effect of different aggressive environments on the tensile properties of FRPs. An interesting outcome was that the carbon FRPs, when exposed to most environments, revealed excellent durability. Even after 27 months of exposure, the mechanical properties of the specimens were minimally affected by these environments. However, the exposure effects were significantly higher in the glass FRP specimens. After 20,000 h of exposure, i.e., about 830 days, the ultimate normalized strength of the carbonfiber specimens decreased by around 10%, while the glass-fiber ones decreased by 30 to 40%. In another study, Bian et al. [35] verified that the tensile strength of glass FRPs was reduced by 13% after only 42 days of immersion. Similarly, Gellert and Turley [43] conclude that the flexural strength continued to degrade for GFRPs as water uptake continued toward saturation. Moreover, the authors reported strength losses between 15 and 21% as well as a decrease in the interlaminar shear strength by 21% after 485 days of immersion.
Despite the detrimental effect of hostile environments on mechanical behaviour of FRPs being well documented in the scientific literature, the development of numerical tools capable of correlating the exposure time with the mechanical performance is much less common. Gholami et al. [44] conducted a micro-scale finite element analysis concerning the hydrothermal degradation of FRP. However, the scope of the study was limited to predicting the properties of the composite using representative volume elements. More recently, Feng et al. [45,46] addressed the effect of hydrothermal environments, via both experimental and numerical approaches, using a multiscale modelling method. Nevertheless, the effects of permeable and non-permeable surfaces on mechanical behaviour of immersed structures were not considered.
Therefore, although the development of the above-mentioned numerical tools are a major challenge, from a perspective of engineering design, such tools can help to improve the overall design of marine structures, reduce safety factors, partly reduce the expensive and time-consuming experimental testing campaigns, and speed up the time to market. With the advent of computer technology, advanced numerical methods have been a strong ally in the analysis of complex phenomena. Due to its flexibility, and because it is widespread in current engineering industry, the Finite Element Method (FEM) is the ideal choice to devise this tool.
In this context, this paper aims to develop a three-dimensional (3D) numerical model capable of predicting the effect of seawater on mechanical strength of glass fiber-reinforced polymers (GFRP) using the FEM. Puck's failure criterion, a failure theory for unidirectional fiber-reinforced polymer composites, was applied to define the ultimate failure of the laminates. Then, Fick's first diffusion law was used to predict the seawater absorption rate. Finally, a relation between the Fick's law and the diffusion distance was introduced to predict the seawater concentration into the material. The numerical model was implemented in the ABAQUS commercial software by developing a user-defined material (UMAT) subroutine. After that, the Digital Image Correlation (DIC) technique was used to validate the implementation of Puck's theory by providing full-field strain measurements. Finally, the predictive capabilities of the proposed approach were compared with experimental tensile test results obtained for different seawater immersion times ranging from 0 to 900 days.
The first part of this paper features the materials and the procedures used to validate the numerical models. The second part outlines Puck's failure criterion, the proposed numerical model, and their implementation in the subroutine. After that, the numerical simulations are compared to the experimental stress-strain test results performed for immersed and not immersed glass fiber-reinforced polymer composite laminates. Finally, the proposed procedure is used to improve structural design by applying impermeable coatings.

Material and Experimental Methods
Composite glass/epoxy laminates were used in this research. The laminates were produced using 1195P (195 g/m 2 ) glass fiber and Biresin ® (Sika Services AG, Zurich, Switzerland) CR122 epoxy resin with the Biresin ® CH122-3 hardener. The amount of hardener was 30% of the resin's total weight, as suggested by the manufacturer. The laminates were placed in a vacuum bag and compressed with a load of 2.5 kN for 12 h to preserve the original thickness. During the first 4 h, the vacuum bag remained connected to a vacuum pump to eliminate air bubbles. The post-cure was conducted in agreement with the manufacturer in an oven at 60 • C for 8 h. The specimens were cut from the laminate panels with dimensions of 165 mm × 22.5 mm × 2.3 mm. Then, a central hole measuring 5 mm in diameter was drilled using a twist drill bit operated at a suitable cutting speed to prevent delamination. After the manufacturing process, three specimens were immersed in seawater for 900 days, and another three were used as control samples in a dry condition, i.e., without immersion in seawater. The former group was weighted regularly during the immersion stage to evaluate the water absorption rate. In total, six tensile tests were conducted, three were used as control samples (i.e., not immersed in sea water), and the other three were used to evaluate the effects of seawater (i.e., immersed 900 days in seawater). The tensile tests were performed at room temperature, according to the ASTM D3039 standard, with a displacement rate of 2 mm/min. A Shimadzu Autograph AGS-X (Shimadzu Scientific Instruments, Kyoto, Japan) machine with a load capacity of 150 kN was used. The stressstrain response was acquired from a 25 mm-long mechanical extensometer connected directly to the specimen using the Trapezium X software (version 2.2).
The strain fields at the specimen surface were obtained using Digital Image Correlation (DIC) using a VIC-3D system (Correlated Solutions, Irmo, SC, USA). The specimens were previously prepared and painted with a white matte base coat, and a random speckle pattern of black paint was applied to create a pattern on the center zone of the specimen. To capture the strain field, two high-velocity cameras (Correlated Solutions, USA) were placed at approximately 30 • , symmetrical in relation to the specimens, at the maximum resolution of 1624 pixels × 1224 pixels. During the loading stage, images were captured using these two cameras, and strain fields were obtained with specialized software. By comparing successive images of the random speckle pattern acquired during the loading history on the surface of each specimen, DIC provided the full-field displacement and strain distributions.

Failure Criterion and Proposed Seawater Exposure Damage Model
This section is divided into three parts. Firstly, the Puck's failure criterion is introduced. Then, the continuum dagame model used to predict the interlaminar progressive failure is addressed. The section ends with a proposed approach to simulate the effect of seawater on laminate.

Puck's Failure Criterion
The initial failure theory was based on Puck's failure criterion [47]. Puck's criterion can identify the different fracture modes occurring in FRPs and quantify the repercussions of the lamina's damage on the laminate. Furthermore, it can also predict the fracture plane's orientation. Puck's failure theory differentiates two types of fracture: Inter Fiber Fracture (IFF) and Fiber Fracture (FF). The IFF criteria are identified using stresses of the action plane since it is not possible to formulate a physically-based fracture criterion using the stresses based on the natural axes of the unidirectional (UD) composite [48]. Therefore, UD lamina stresses in the lamina coordinate system (1,2,3) are translated into the fracture plane coordinate system (1,n,t) using tensor transformation in the 2-3 plane: where θ is the inclination of the fracture plane, σ 2 and σ 3 are normal stresses perpendicular to the fiber, and τ 21 , τ 31 and τ 23 are shear stresses. The verification of IFF requires the previous identification of the fracture plane. Hereby, it is mandatory first to detect the plane's fracture angle θ f p . The fracture plane is defined as the action plane with the maximum local stress exposure, and the local stress exposure, f E , can be used as a direct measure of the fracture risk since it grows linearly with the stresses [49]. Therefore, θ f p is searched numerically, calculating f E for all sections between θ = −90 • ; 90 • to find the fracture plane. The local stress exposure can be obtained as follows: where: where R A ⊥⊥ , R A ⊥ or R At ⊥ are the fracture resistances of the action plane to a single transverse shear stressing τ ⊥⊥ , longitudinal shear stressing τ ⊥ , or transverse tensile stressing σ t ⊥ , acting on the action plane, respectively. R ⊥ and R t ⊥ are the longitudinal shear strength, and the transverse shear strength of the lamina, respectively, while p t ⊥ , p c ⊥ p t ⊥⊥ , and p c ⊥⊥ are experimentally determined inclination parameters provided by Puck et al. [50] and summarised in Table 1. After finding the fracture plane, Puck's failure criterion can be verified with the following expressions: Fiber failure in tension: Matrix failure in compression (σ n < 0): Matrix failure in tension (σ n ≥ 0): where ν ⊥ is the major Poisson's ratio of the fiber, E is the lamina longitudinal modulus parallel to the fibers, while E f is the longitudinal fiber modulus and m σ f is the magnification factor for the transverse stress in the direction of the fiber. In this work, m σ f was set to 1.3, as recommended for GFRP by Lee et al. [49]. Note that 1-direction corresponds to the fiber orientation, 2-direction is transverse in-plane orientation relative to the fiber direction, and 3-direction stands for transverse out-of-plane orientation.

Constitutive Relation and Progressive Failure Theory
Numerous progressive failure theories have been proposed to predict defect initiation and propagation in common materials. However, most damage models were designed for isotropic ductile materials, and their implementation is complex for anisotropic brittle materials such as FRP composites. Gu et al. [51] proposed a continuum damage model to predict the interlaminar progressive failure of CFRP laminates. A bilinear relation was used to describe the damage accumulation process. A good agreement between the experimental and numerical results was verified. In the present study, which is inspired in the work of Lee et al. [52], a stress-based Puck failure criterion was adopted to predict the damage progression. This approach is simple and widely used for strength estimation and progressive failure analysis.
The constitutive relation, which relates stress states to strain states, is written as follows [53]: where ε I and γ ij represent three normal strains and three shear strains, respectively; I and τ ij denote three normal stresses and three shear stresses, respectively, and G ij denote the shear moduli. In addition, C ij can be defined with the Young's moduli and Poisson's ratios as will be established later in this section.
The constitutive relation parameters change based on the Element Weakening Method (EWM), if failure is verified [52]. The constitutive relation is modified by the damage parameters d f and d m (Equations (18) and (19)) associated with the fiber and matrix failure, respectively. The variables d f t , d f c , d mt , and d mc are related to fiber and matrix damage under tensile and comprehensive stress states, respectively. These variables assume the value of 1 if their associated failure mode occurs, as summarized in Table 2. All parameters for the constitutive relation, fiber and matrix damage variables, and material constants for the undamaged and damaged material stiffness matrix are written as follows: Table 2. Damage values for each failure mode.

Failure Mode Damage Value
Tensile fiber (20) where C 0 ij and C ij are the undamaged and damaged material stiffness components, respectively; E 1 , E 2 and E 3 are Young's moduli in the three laminae directions (x 1, x 2, x 3 ), respectively; ν ij are the Poisson's ratios that can be defined as ν ij = −ε j /ε i , with i = j, and G 12 , G 23 and G 31 are shear moduli in x 1 − x 2 , x 2 − x 3 and x 3 − x 1 planes, respectively. Lastly, s mt and d mt are the loss control factors for the shear stiffness caused by the matrix tensile and compressive failures, respectively. In this study, the used loss control factors values were: s mt = 0.9 and s mt = 0.5 as used in work by Seung Lee et al. [52].

Proposed Model for Seawater Damage
In the present study, a numerical 3D model for predicting the exposure damage to seawater was developed. The proposed approach states that the damage can be modelled by changing the lamina's resistances as a function of the water concentration. Here, it is assumed that the water concentration of the finite elements closer to the permeable surface, C CS , decreases linearly with the distance to the closer permeable surface, as the following equation suggests: where D CS is the Gaussian point's distance to the closer permeable surface, and D dist is the diffusion distance. The lamina's resistances parallel and perpendicular to fiber R t,c ⊥, are degraded as a function of the concentration: where k is an experimental fitting parameter. Concerning the increase in stiffness, Young's modulus increases as a function of the normalized Fick's law. Unlike the resistances, which can be a function of the concentration, Young's moduli cannot, since it would create a fictitious stress concentration. Therefore, Young's moduli are a function of the normalized Fick's law as demonstrated in the following equation: where E i is the Young's moduli in the direction i, E 0,i is the original Young's moduli in the direction i; N Fl is the normalized Fick's law, a function of the exposure time in days; and A, and B are fitting parameters. The variation of the diffusion distance, D dist , is based on Fick's law. Hence, a linear relation between the weight gain, in Fick's law, and the wet volume, W volume (D dist ), was assumed. The wet volume is the specimen volume affected by the diffusion and is a function of the diffusion distance. With this assumption, it is possible to find a function for the diffusion distance that forces the normalized wet volume to have the same behaviour as the normalized Fick's law. Thus, it can be found a function for the diffusion distance that varies with the immersion days and is based on Fick's law. When saturation is reached, the diffusion distance must be such that the wet volume equals the structure volume. Accordingly, the diffusion distance function is composed of two terms, as shown: where A(t) is a continuum function that varies between [0,1] with the immersion time in days, t, which must be modified for the normalized wet volume to have the same behaviour as the normalized Fick's law. C is the structural dimension, so that W volume (D dist = C)= structure volume. It should be noted that finding a function for the wet volume varying with the diffusion distance, W volume (D dist ) can be relatively tricky particularly for complex structures. However, if the wet volume increases linearly with the diffusion distance, then the diffusion distance can be a direct function of the normalized Fick's law, i.e., A(t)= Normalized Fick s law. Many studies have shown that the weight gain in composite materials can be estimated using Fickian models [37,41,54,55]. Therefore, in this study, a one-dimensional Fickian model was used. The moisture absorption contends that M t is time-dependent and can be estimated by the following equation [56,57]: where M ∞ is the maximum amount of absorption, t is time in days, h is the thickness, and D Coe f is the diffusion coefficient. The following equation allows obtaining the diffusion coefficient D Coe f , using two points at time t 1 and t 2 from the linear part of Fickian diffusion curve, i.e., where k can be obtained as suggested by Wang et al. [14]: where M 1 and M 2 are the percentage of water absorbed at time t 1 and t 2 , respectively. The parameters are summarised in Table 3, and the comparison between the fitted Fick's law and the experimental results is displayed in Figure 1. Overall, as can be seen in the figure, the proposed function is close to the experiments regardless of the immersion time. The maximum differences, which were found for immersion times of 640 days, were lower than 8.8%.

Finite-Element Based Implementation
The iterative procedure implemented in this study to simulate the effect of seawater exposure on mechanical properties in FGRPs was developed using the Abaqus commercial software package along with a UMAT subroutine written in Fortran. The algorithm framework is presented in Figure 2. As can be seen in the figure, at each time increment for each integration point, the information on material properties, strain tensor, strain incrementation and state variables are passed to the UMAT subroutine. After sending that information, if the iteration is the first one, the waterdamage routine is called. This condition reduces simulation time since this routine requires data that do not vary with the increment number. Concisely, the routine verifies if the integration point is within the zone damaged by the seawater, and if so, the tensile and compressive strengths parallel and perpendicular to fiber, R t,c ,⊥ , are degraded. Then, the constitutive matrix is evaluated, and following that, the stress tensor is calculated. In the next step, the SearchFP routine is called, which calculates the value of (see Section 3.1). After that, both the IFF and the FF failure criteria are verified (see Section 3.1). The respective degradation law is applied if a failure occurs. Since the constitutive matrix can change in the case of failure, the failure criteria must be verified again. Finally, the state variables are updated, and the updated constitutive matrix, the stress tensor, and the state variables are sent to Abaqus. Convergence is achieved if the residual force is lower than the tolerance value. If it is higher, the increment is repeated with a lower increment time to ensure easier convergence.
The finite element mesh of the tested geometry, see Figure 3, was generated with one element through-thickness for each lamina using 8-node linear hexahedral elements with reduced integration (C3D8R). The assembled model had 123,262 elements and 146,874 nodes. This element type uses a single integration point located at the centroid of the element which provides excellent benefits, since the number of evaluations during the simulation is substantially reduced. However, the difference between the number of nodes and the number of integration points can lead to mesh instability, known as hourglassing. Hence, the enhanced hourglass control option was activated for all simulations. Then, the constitutive matrix is evaluated, and following that, the stress tensor is calculated. In the next step, the SearchFP routine is called, which calculates the value of θ f p (see Section 3.1). After that, both the IFF and the FF failure criteria are verified (see Section 3.1). The respective degradation law is applied if a failure occurs. Since the constitutive matrix can change in the case of failure, the failure criteria must be verified again. Finally, the state variables are updated, and the updated constitutive matrix, the stress tensor, and the state variables are sent to Abaqus. Convergence is achieved if the residual force is lower than the tolerance value. If it is higher, the increment is repeated with a lower increment time to ensure easier convergence.
The finite element mesh of the tested geometry, see Figure 3, was generated with one element through-thickness for each lamina using 8-node linear hexahedral elements with reduced integration (C3D8R). The assembled model had 123,262 elements and 146,874 nodes. This element type uses a single integration point located at the centroid of the element which provides excellent benefits, since the number of evaluations during the simulation is substantially reduced. However, the difference between the number of nodes and the number of integration points can lead to mesh instability, known as hourglassing. Hence, the enhanced hourglass control option was activated for all simulations. Mesh convergence is another critical subject in any finite element analysis. In general, mesh convergence problems arise in regions with severe geometric discontinuities [58]. To prevent this concern, the element size around the hole was reduced to a 0.1 mm side, ensuring convergence in terms of stress and strain. In remote positions, the element size was larger to reduce the computational overhead. Regarding the type of finite element analysis, although a quasi-static analysis is particularly stable for large time steps, using lower time steps can lead to improved solutions. However, it can result in higher simulation times. Thus, each simulation was executed with a time step of 0.01 s, while the total time of 1 s was simulated. The material properties used in the simulations were those recommended by the WWFE [59].
The boundary conditions used to simulate the tensile tests were adopted using kinematic constraints in zones A and B, represented in Figure 3. In the former zone, all translations and rotations were restrained, while in the latter zone only the translation on the x-axis direction was allowed. The loading was applied on constraint B in the x-axis direction.

Results and Discussion
The first objective of the present study was to validate the implementation of Puck's failure theory. Thus, the strain fields near the hole region simulated with the proposed approach were compared with those obtained in the experimental tests for different loading levels [60]. The experimental results, as referred to above, were measured using Digital Image Correlation (DIC) [61,62]. The comparative analysis was carried out using the first principal strain. As depicted in Figure 4, the numerical simulations and experimental results agree well, regarding the strain fields. In fact, when comparing DIC results to the implemented Puck's failure model, the strain profile is similar for the three loading levels which validates the proposed approach.
Furthermore, although during the tensile tests loads were uniformly applied on the external surfaces of the specimen during the tensile tests, from the experimental results, it is possible to observe that the first principal strain field at the hole surface is not symmetrical. Puck's failure theory can predict this asymmetrical behaviour, which is likely to happen due to the laminae sequence, resulting in an asymmetrical matrix damage evolution. Figure 5 compares the maximum first principal strain with the matrix damage as load increases. At 3.0 kN or lower, no matrix damage occurs at the surface and no asymmetry is perceptible in the strain field. As the load increases, matrix damage starts to build up, and the strain asymmetry is lightly perceptible. At 5.0 kN or higher, a pronounced asymmetrical matrix damage pattern occurs, generating the uneven strain pattern that is verified in the experiments. After the validation of the strain fields near the hole region, which is based on the implemented Puck's failure theory, the tensile stress-strain tests for the control samples, i.e., the specimens not immersed in seawater, were also compared with the numerical simulations. The maximum load and displacement at failure for the tensile tests for the control and the immersed specimens are summarized in Tables 4 and 5, respectively. The results for the maximum load are, for all the tests, relatively close to the expected value; however, the maximum displacement results revealed some scatter. The stress-strain plots determined numerically and obtained in the experiments are shown in Figure 6. In all simulations, the final failure is assumed to occur when the fiber failure crosses all the specimen width, accompanied by a sudden escalation of the displacement.  The numerical model predicts the failure stress with an error of 8.3% and the strain at failure with an error of 10.7% compared with the average experimental results. Even though the model numerical model predicts the final values with reasonable precision, the stress-strain plot exhibits a nonlinear behaviour at a stress of approximately 80 MPa and then resumes the initial linear behaviour until the total failure. Although there is a slight nonlinearity on the experimental curve, it is noticed a small disagreement between the experimental data and the numerical behaviours. This nonlinear behaviour is caused by the first IFF (see Section 3.1), which is documented in the work of Knops [49].
The numerical model can also provide important insights into the fundamental mechanisms associated with the damage development. The experimental evidence, as already verified by Gonçalves [60] and Aguiar [63], suggests that the borderline of the hole is where the cracks start to spread, which is substantiated by the numerical simulations.   Regarding the fiber failure, as can be distinguished in Figure 8, it only takes place on the 0 • plies, since the other pliers' fibers are not being loaded. At about 7.0 kN, the fiber damage appears, and then, as the load increases, it propagates along a 90 • angle relative to the loading direction. At last, the test specimen is fully fractured under approximately 9.7 kN, and it is clear that the failed region induced by IFF is much wider than that induced by FF. Figure 8 also shows the uniaxial stress-strain curves obtained experimentally and simulated numerically for 900 days of immersion considering the methodology proposed to account for the damage caused by exposure to seawater. Overall, both curves are relatively close in the entire range. The numerical model predicts the failure stress with an error of 2% and the displacement at failure with an error of 4.9%. Furthermore, the initial nonlinear behaviour of the curve, already distinguished for the specimens not immersed in seawater, is also visible in this case. Similarly, it is also softened with the decrease in the ductility of the tested composite. As can be seen in Figure 9, the mechanical behaviour of the laminate is accurately predicted for applied stresses higher than 80 MPa, which is an interesting outcome and demonstrates the adequacy of the proposed methodology. Figure 9 illustrates the seawater concentration after 30, 60, 150, 400 and 900 immersion days. It can be noticed that, with the increase of the immersion days, the diffusion distance increases, causing the diffused regions to cover the specimen's volume to a greater extent. It is interesting to observe that the concentration is identical from day 400 to 900, although with more than twice the exposure time. This result is expected since the weight gain in the last section of Fick's law remains approximately constant.  Figure 10 displays the stress-strain plots simulated with the proposed numerical model for immersion days in the range 0 to 900 days. It is clear that most of the strength and ductility loss occurs in the first days of immersion. The differences between the cases corresponding to 300 and 900 days of immersion are not significant, particularly with regard to the maximum stress; regarding the strain at failure, although the results are not much different, there is a more relevant reduction with the increase of the immersion time. Nevertheless, since most of the water absorption occurs in the first months of immersion, see Figure 1, it is expected that most of the damage also develops during this period. The relationship between the maximum applied stress and the immersion time obtained by means of the numerical model developed in this study is exhibited in Figure 11. As inferred above, it can be noticed that most of the damage occurs in the first three months of immersion, and after approximately 250 days, no significant changes occur in the failure stress. Based on the numerical simulations, the dialectical relationship between the failure stress and the immersion time was established via nonlinear regression: where S is the failure stress, and ID is the number of immersion days. Figure 11 also plots the experimental values of the failure stress obtained for 0 days and 900 days of immersion.
In this representation, it can also be observed that the simulated results agree well with the experiments. Another important application of the proposed methodology is the design improvement of engineering components subjected to hostile environmental media. For example, we can consider different surfaces as permeable, or impermeable, and simulate the effect of these changes on structural response. This can be extremely important to increase the long-term durability because surface protective coatings degenerate over time and need to be renewed. Since diffusion phenomena can happen on different surfaces and can have different effects on mechanical strength, the numerical model developed in this study is a valuable tool to identify a priori the most critical surfaces allowing an optimized design to be achieved in a simple and fast manner, but also allows to simulate the structural behaviour and quantify eventual benefits associated with different combinations of impermeable and permeable surfaces. Figure 12 compares the seawater concentration considering all permeable surfaces (Figure 12a), impermeable lateral surfaces (Figure 12b), and an impermeable hole surface ( Figure 12c). As can be seen in the figure, the protection against damage by seawater of the hole surface significantly reduces the seawater concentration in the central region around the hole (Figure 12c). On the contrary, in the case of the hole surface not being protected and the lateral surfaces being simultaneously protected (Figure 12b), the seawater concentration is significantly higher around the hole region even though it is necessary to protect a much wider area. When no protection is used, high seawater concentrations are found not only at the lateral surfaces but also at the hole region. The different seawater concentration distributions reported above are likely to affect the mechanical behaviour of the laminate. Figure 13 compares the uniaxial stress-strain response considering the non-permeability of the lateral surfaces and the non-permeability of the hole surface for 50, 150 and 900 immersion days. Comparing the model predictions between the simultaneously impermeable lateral surfaces (165 mm × 2.3 mm) and the impermeable hole surface, the difference in failure stress is very significant for the three immersion periods. For 50, 150 and 900 immersion days, the numerical simulations indicate that the specimen can sustain 10.4, 16.2 and 11.8% more load before failure if the hole surface is protected, respectively. This is a coherent outcome since, with an impermeable hole surface, no reduction in mechanical properties occurs near the hole where the fiber damage initiates (see Figure 8). As a result, the fiber damage only initiates at higher stresses enabling the specimen to support more load. In terms of structural integrity, when statically loading the specimen, protecting the hole surface is more efficient than protecting the lateral surfaces, even though the area of the lateral surfaces is more than 21 times the area of the hole surface. Thus, an efficient way to improve the structural strength in a marine environment for structural components made of GFRPs and having holes would be the application of protective coatings on the hole surface.

Conclusions
This paper presented a 3D numerical procedure for predicting the effect of seawater exposure on mechanical behaviour of glass fiber-reinforced polymers. The following conclusions can be drawn:

•
The linear relation between the weight gain and the wet volume was adequate to model the diffusion process. The maximum errors between the fitted function and the experimental results were lower than 9%. • The first principal strain fields near the hole surface simulated with the proposed methodology were quite similar to those measured in the experiments for the different loading cases with digital image correlation; • The model predicted both the maximum stress and the strain at failure with good accuracy, irrespective of the immersion time. In the absence of seawater, the results were more accurate; • The stress-strain plots exhibited a nonlinear behaviour. However, this non-realistic behaviour only occurred in the early stage of the loading process. After a certain stress level, the numerical simulations and the experimental results behaved similarly; • A power relationship between the stress at failure and the immersion time was found. Most of the damage caused by the exposure to seawater occurred in the first three months; • The model can help designers to identify the areas where the application of impermeable coatings can improve the resistance to seawater exposure. The cases analysed showed gains in maximum stress higher than 10%. Informed Consent Statement: Not applicable.

Data Availability Statement:
The raw data presented in this study are available on request from the corresponding author.

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