Effect of moisture on the nonlinear viscoelastic fracture behavior of polymer nanocompsites: a finite deformation phase-field model

The mechanisms underlying damage in high-performance polymer nanocomposites are remarkably affected by hygrothermal conditions. In this study, we develop a phase-field formulation to investigate the influence of hygrothermal conditions on the nonlinear viscoelastic fracture behavior of epoxy resins and their nanocomposites at finite deformation. For this, the Helmholtz free energy, capturing the effect of temperature and moisture and nanoparticle contents, is defined based on an additive decomposition of the energy into an equilibrium, a non-equilibrium, and a volumetric contribution with different definitions under tensile and compressive loading. The coupled displacement phase-field problem is solved using a quasi-Newton monolithic algorithm and a staggered solution scheme. Numerical examples show that the monolithic algorithm is more efficient. Simulations are performed to investigate the effect of temperature, deformation rate, and moisture content on the force–displacement response of boehmite nanoparticle/epoxy samples in benchmark numerical problems. Comparing numerical predictions and experimental data for compact-tension tests shows good agreement at different nanoparticle contents. Also, the model’s capability to predict fracture patterns is evaluated using simulations of single-edge notched nanocomposite plates under tensile and shear loading.

The thermal expansion and moisture swelling coefficients Absolute temperature np , w w BNP volume and moisture weight fractions a w , b w Parameters to define the effect of moisture The Cauchy stress D The symmetric Eulerian rate of the deformation tensor L The Eulerian gradient of velocity The Helmholtz free specific energy ± Positive/negative parts of the free specific energy eq , neq , vol Equilibrium, non-equilibrium, and volumetric parts of the free specific energy The phase-field parameter g( ) The energetic degradation function V e ,R e The pure deformation and rotation obtained from the polar decomposition of F ē D i ,W i The rate of stretching and spin D 0 i The objective rate of inelastic deformation

Introduction
One of the main challenges in today's engineering is structural weight reduction to provide higher performance and functionality for specific applications. Therefore, besides structural optimization, research is also focused on developing new materials with enhanced thermo-mechanical properties at low weight. One of these new classes of materials is polymer nanocomposites, where the functionalities of polymers, including low weight and high ductility, are combined with the nanoparticles' unique features [37,57]. Recently, boehmite nanoparticle (BNP) reinforced epoxy composites have been considered as one of the most promising composites in lightweight structures due to their high strength-to-weight ratio [31]. BNP/epoxy nanocomposites exhibit remarkably improved mechanical properties, including strength and fracture toughness, compared with neat epoxies [7,31].
The material innovation demands reliable models to predict the effect of external conditions (e.g., loading rate, temperature, and moisture) and microstructural parameters (e.g., nanoparticle/matrix interactions) on the damage and fracture behavior of the nanocomposites. Continuing research activity on polymers and their composites has led to a variety of phenomenological or physically motivated constitutive models [48,50] to elucidate their nonlinear rate-and temperature-dependent behavior. Boyce et al. [12] developed a constitutive model based on a composite-type formulation considering the microstructure of semicrystalline polymers. In the model, the soft amorphous and stiff crystalline phases are treated as the matrix and fillers, respectively. Later, based on the model, Qi and Boyce [53] proposed a viscoelastic-viscoplastic constitutive model to capture the nonlinear, rate-dependent, and softening behavior of thermoplastic polyurethanes. Li et al. [36] introduced a physically based viscoelastic constitutive model for elastomers at large deformation, where elastomers are assumed to be crosslinked networks with superimposed free chains. Nguyen et al. [48] developed and experimentally calibrated a ratedependent damage constitutive model for epoxy resins to study the nonlinear behavior of amorphous glassy polymers. Based on the definition of Helmholtz free energy, N'Guyen et al. [49] derived a thermodynamical framework for the thermo-chemo-mechanical couplings in polymer materials at finite deformation. Predicting the nonlinear stress-strain response of polymer nanocomposites, which also contain nano-scale additives, is more challenging due to the heterogeneous distribution of agglomerated nanoparticles in the matrix and complex interactions between the matrix and nanoparticles. Fankhänel et al. [20] presented an atomistically informed finite-element (FE) model to study the material properties of BNP/epoxy nanocomposites. Within the multiscale model, the interphase properties between BNPs and an epoxy matrix were first characterized using molecular simulations. The interphase properties were then upscaled to the continuum scale, where the effective material properties were homogenized using FE simulations of representative volume elements of the nanocomposite. Arash et al. [6,7] proposed a combined simulation-experiment framework to calibrate a viscoelastic damage model for BNP/epoxy nanocomposites at finite deformation. The experimentalnumerical validation proves the predictive capability of the model to capture the main features of the stress-strain relationship of the nanocomposites, including the nonlinear hyperelastic, rate-dependent, and softening behavior. The experimental-numerical validation proves the capability of the model to capture the main features of the stress-strain relationship of the nanocomposites, including the nonlinear hyperelastic, rate-dependent, and softening behavior. In seeking a robust parameter identification procedure, Unger et al. [59,60] extended the multiscale approach to characterize the thermo-viscoelastic damage behavior of BNP/epoxy nanocomposites. These studies show that the proposed simulation-based framework allows significantly reducing the number of experimental tests required for identifying the material parameters used in modeling polymer nanocomposites without significant loss in accuracy.
As explained in [15,48], the evolution of damage in polymer materials at the small scales, such as microvoids and microcracks coalesce, leads to the birth of cracks at the macro-scale. The progressive evolution of damage tends to the localization of deformation into a narrow zone accompanied by a softening behavior. Therefore, the modeling of fracture in nanocomposites requires an accurate prediction of damage initiation and propagation in the materials. However, FE models based on local continuum description of damage [6,7] suffer from an inherent mesh dependence for strain softening problems [22,47], resulting from illposedness of the boundary value problem due to the loss of ellipticity in statics or hyperbolicity in dynamics. To ensure the well-posedness of the boundary value problems, various computational techniques have been proposed in the literature. The so-called regularized solutions for damage and failure in materials, allowing interactions between neighboring material points, are among the most successful. The solutions link the continuum damage mechanics to fracture mechanics through the coupling of a diffusion-type equation of nonlocal variables with the momentum balance equation. A nonlocal variable, such as damage, is related to the spatial average of the field of the variable in a certain neighborhood of a given point [9]. In these methods, such as the gradient-enhanced damage model and its variants [8,51,52,61], and the phase-field model (PFM) [3,21,44], a sharp crack is approximated by a diffuse damage band thanks to introducing a length-scale parameter that controls interactions between material points.
Generalizing Griffith's theory, PFMs have emerged as variational fracture models to adequately predict the crack initiation, propagation, and branching [1,2,54,65,67]. In these models, a fracture can be revisited as the minimization of the potential energy consisting of the stored bulk energy, the work of external forces, and the surface energy. In addition, an auxiliary variable, the so-called phase-field parameter, describes a smooth transition from an intact material to a fully broken state. PFMs have been used to study brittle fracture [11,43], quasi-brittle fracture [19,63], and ductile fracture [1,42]. Among others, Msekh et al. [45,46] developed PFMs for clay/epoxy nanocomposites, where clay nanoparticles were modeled as linear elastic materials embedded in a hyperelastic matrix. The models were also used to investigate the surface energy dissipation in the nanocomposites during fracture. Recently, Goswami et al. [25,26] proposed a neural network algorithm for phasefield modeling fracture in brittle materials by minimizing the variational energy of a system. The simulation results show that the crack path predicted by the proposed approach is in agreement with those reported in the literature. Furthermore, some PFMs have been developed to study the rate-dependent fracture of solids [14,38,55,65]. Shen et al. [55] derived phase-field formulations for fracture of viscoelastic materials at small deformation. Loew et al. [38] calibrated a ratedependent PFM for rubbers, where the material parameters were extracted using uniaxial tensile and double-edge tensile tests. Moreover, digital image correlation determined the length-scale parameter by measuring local strain at the crack tip. Yin and Kaliske [65] integrated a viscoelastic model into a PFM to study the rate-dependent fracture behavior of elastomers. In this model, different from those proposed in [55] and [65], the crack driving force does not include the viscous energy dissipation. Instead, the driving force is defined by the elastic strain potential given by the equilibrium and non-equilibrium networks of the viscoelastic model. Brighenti [14] proposed a PFM for elastomers using a statistical physics-based micromechanical model, which captures the rearrangement of a polymer network over time. Due to the non-convexity of the total potential energy functional with respect to the kinematic variables in PFMs, hindering convergence and robustness in Newton method-based monolithic methods, different numerical strategies have been suggested to overcome the drawback [18,23,28,40,56]. Among them, staggered solution schemes [11,43] have been shown to be robust to solve the coupled damage-displacement governing equations, but they offer computationally expensive solutions.. Recently, to remedy the computational overheads, a quasi-Newton method-based monolithic algorithm has been proposed in [34,64] as a robust and numerically efficient solution scheme for phasefield fracture modeling.
As for thermoset nanocomposites, nanoparticle-polymer matrix interactions and ambient conditions (e.g., temperature and moisture) affect the evolution of damage [17,30,60]. In this contribution, to address the open questions, we derive a phase-field formulation for investigating the effect of hygrothermal conditions on the rate-dependent fracture behavior of BNP/epoxy nanocomposites at finite deformation. The corresponding free specific energy is defined on the basis of a volumetric-deviatoric decomposition of the total deformation gradient, where the volumetric part is split into positive/negative components. Modified Guth-Gold and Kitagawa models are adopted to study the influence of the nanoparticle and moisture contents and temperature on the crack propagation in the polymer nanocomposites. The coupled governing equations are solved using a quasi-Newton monolithic algorithm and a staggered solution scheme. Numerical examples confirm that the monolithic algorithm yields identical results to the staggered solution with higher efficiency. The effect of temperature, deformation rate, and moisture content on the force-displacement behavior of the nanocomposite samples are studied using the numerical simulation of benchmark problems. Also, the proposed model is validated through the comparison of numerical predictions obtained from the modeling of compact-tension (CT) tests with experimental data.
This work is organized as follows. Sect. 2 presents a nonlinear viscoelastic model describing the temperature-and moisture-dependent behavior of polymer nanocomposites at finite deformation. The governing equations of the PFM and the corresponding discretized equations are provided in Sect. 3. In Sect. 4, the proposed PFM is validated using the numerical simulations, and the effect of temperature, deformation rate, and moisture on the fracture behavior of the nanocomposites is investigated. Finally, Sect. 5 summarizes the findings.

Constitutive model for nanoparticle/ epoxy
The stress response of a nanoparticle/epoxy system shown in Fig. 1 can be decomposed into an equilibrium part and a viscous part to capture the nonlinear and rate-dependent behavior of the material at finite deformation [7]. In the following section, a constitutive law, which models the nanocomposites as homogeneous continua by neglecting complex interactions between nanoparticles and the epoxy matrix, is presented. The approach is suitable for modeling the nanocomposites at the macro-scale.

Kinematics
The total deformation gradient is multiplicatively split into a volumetric and a deviatoric part to define the kinematics as [49] where J = det [F] and F are the volume variation and the isochoric deformation gradient, respectively. Here, it is assumed that the volume variation can be decomposed into three terms: the mechanical compressibility J m , the thermal dilatation J , and the moisture-induced swelling J w where and The proposed form for J w to define the swelling is based on the simulation results presented in [16]. In the above equations, and w are, respectively, the thermal expansion and moisture swelling coefficients, is the absolute temperature, 0 is the temperature of the reference configuration, and w w is the moisture content. The deviatoric part of the deformation gradient can also be decomposed into elastic and inelastic parts by introducing an intermediate configuration as follows [35]: Accordingly, the total and elastic left Cauchy-Green deformation tensors, related to the reference configuration, are given by and

Clausius-Duhem inequality
To derive the constitutive law coupled to the phase-field evolution within the kinematic framework, we proceed with the Clausius-Duhem inequality for an isothermal process. Combining the first and second principles of thermodynamics, the inequality in an Eulerian configuration reads Here, D is the dissipation, is the Cauchy stress, D = 1 2 L + L T the symmetric Eulerian rate of the deformation tensor, and ̇ is the material time derivative of the Helmholtz free specific energy. L =Ḟ.F −1 represents the Eulerian gradient of velocity. The Helmholtz free specific energy is defined by where cracks are characterized by a phase-field parameter varying between 0 and 1. = 0 denotes an intact material and = 1 represents a fully cracked material. The material time derivative of the free energy is then given by [8] where The pure deformation V e in Eqs. (11) and (12) is obtained from the polar decomposition of F e =V e .R e , and the objective rate of inelastic deformation D 0 i is given by where The velocity gradient of the relaxed configuration L i can be additively decomposed into a symmetric tensor D i and a skew-symmetric tensor W i , so-called the rate of stretching and spin, respectively, such that L i =D i +W i . Since the intermediate configuration can be taken in different ways, a convenient form for the spin of the relaxed state without loss of generality is W i = 0 [10,13]. By substituting Eqs. (11)-(13) into Eq. (10), we arrive at Finally, introducing relation (16) into (8), we obtained the dissipation inequality To proceed further, we make two assumptions based on Germain's work [24]: D m and D are independently positive, and the dissipation results from the thermodynamic flux or force related to the internal variable F i . Based on the assumptions, the thermodynamic force related to the flux D has to remain null regardless of the thermodynamic process. Therefore, the the stress can be split into equilibrium eq , non-equilibrium neq , and volumetric vol terms as follows: where superscript D denotes the deviatoric operator. The following terms remain in the dissipation: Assuming that two terms in Eq. (19) . Equation (20) is a natural choice for D 0 i to make the first term of Eq. (19) a positive-definite quadratic form. Using Eq. (14), it can be remarked that and the viscoelastic flow ̇i is defined by the Argon model [62] where k b is the Boltzmann constant. The model is characterized by a pre-exponential factor ̇0 , the activation energy H , and the athermal yield stress 0 . Recently, Unger et al. [60] showed that the Argon model leads to a good agreement with experimental data in a wide range of temperatures. Substituting Eq. (21) into Eq. (15) gives (19), the evolution of damage is an energy dissipation process (i.e., −̇≥ 0 ).
Here, the midpoint method is used to numerically obtain the inelastic deformation gradient at the end of a time increment, that is To calculate the elastic deformation gradient at the midpoint, it is required to find the total deformation gradient at the midpoint. This is done by taking the average of the deformation gradient at the start and end of the increment

Phenomenological viscoelastic model coupled with a phase-field description
Following the additive decomposition of the free energy proposed in [1], the overall free energy of the material can be decomposed into an equilibrium eq , a non-equilibrium neq , and a volumetric part vol as: and The Heaviside step function is defined as The energetic degradation function g( ) captures the evolution of the strain energy versus the phase-field parameter and satisfies the following conditions: The conditions prescribe a monotonic decreasing behavior during the fracture evolution. To prevent crack propagation under compression, the volumetric strain energy does not change when J < 1 in Eqs. (28) and (29).
Here, the equilibrium eq and non-equilibrium neq parts of the free energy are defined by the neo-Hookean hyperelastic model as and where I 1 (⋅) = tr[⋅] is the first invariant of the tensor. The material parameters eq and neq depend on temperature, BNP volume fraction np , and water content w w Assuming that BNPs are well-dispersed rigid particles in the epoxy matrix, the Guth-Gold model is adopted by which the effective stiffness of particle-filled solids is obtained by [27]. The amplification factor X is typically a function of fillers' volume fraction and distribution. So far, some attempts of various levels of sophistication have been conducted to incorporate the effect of particle/matrix interactions on the effective modulus of polymer composites. Most of these models suggest a polynomial series expansion for the amplification factor. Here, a modified Guth-Gold model is proposed to account for uniformly distributed nanoparticles and moisture content as follows: In Eq. (36), the effect of moisture content on the material behavior is included based on the experimental data reported in the literature [16]. Accordingly, moisture absorption results in volumetric swelling, which in turn leads to an epoxy network with higher structural mobility and lower stiffness. However, the modulus is partially recovered beyond a certain moisture content. This effect can be explained by type II bound water [66]. Compared with the mobile water (no hydrogen bond) and type I bound water (forming one hydrogen bond), type II bound water molecules form two hydrogen bonds with epoxy chains. As a result, more type II water binds with the polymer chains during the water absorption process, resulting in higher cross-linking and stiffness partially retrieving. In this study, we choose a widely used model and modified the model to capture the effect of BNPs and water molecules as well-dispersed particles in an epoxy system.The modified amplification factor proposed in this work is a first step to capture the stress-strain behavior of BNP/epoxy nanocomposites under hygrothermal conditions. To consider the effect of temperature on the material properties, a modified Kitagawa model proposed by Unger et al. [60] is utilized with the following equations: The volumetric part of the free energy vol is also defined by where t he bulk modulus is assumed to be

Phase-field model at finite deformation
To evaluate the predictive capability of the proposed material model, we use the model to develop a PFM formulation for nanoparticle/polymer composites. This section presents a variational phase-field formulation for fracture at finite deformation. In the following, to show the procedure of analysis, we derive the continuum mechanics incremental scheme and FE equations.

Problem field description
The strong form of the boundary value problem in spatial description for the coupling between displacement u and phase-field variable can be written as with the following boundary conditions: where l 0 is the length scale that controls the width of the diffuse crack, = t ∪ u , b represents the vector of body forces, n is the outward unit normal vector on the boundary of the body Ω , t is the traction force, and u d represents the prescribed displacements at the boundary u . To take into account the effect of BNP content on the fracture evolution in the nanocomposites, the energy release rate is taken to be G c = X np , w w G 0 c . Following Miehe et al. [44], H is the local history field of maximum positive reference energy defined by: which prevents the healing of cracks when the source term + 0 decreases. Here, a monotonically decreasing degradation function g( ) , satisfying conditions presented in Eq. (31), is chosen [29] where k is a small positive parameter introduced for ensuring for the stability of the solution [44].
It should be noted that the propagation of cracks in the benchmark tests investigated in the present study takes a few minutes to a few days depending on the displacement rate. While moisture diffusion into a specimen and the degradation of its mechanical properties due to water (43) .
absorption takes some weeks to some months [16]. It implies that moisture diffusion and crack propagation can be considered decoupled phenomena physically.

Finite-element formulation
To derive the governing equations in weak form, the weighted residual approach is used. Thus, Eq. (43) is multiplied by weight functions and integrated over as and Using the divergence theorem and imposing the boundary conditions .n = t and ∇ .n = 0 presented in Eq. (44), the weak form of the governing equations can be written as follows: Employing the Bubnov-Galerkin method, the displacement, the phase-field field, and the corresponding weight functions are discretized in each element where the shape function matrices N u and N interpolate the nodal values u and , respectively, and B u and B are the gradient operators for the displacement and the nonlocal equivalent strain, respectively. The same shape functions interpolate the nodal values of the weight functions u and . Substituting the relations into the weak formulation of the governing equations yields which have to hold for any choice of u and . The discretized equations can therefore be rewritten as The equations can then be expressed in terms of external and internal nodal forces as where

Consistent incremental-iterative scheme
By linearizing Eq. (56) at iteration i + 1 with respect to the previous iteration i, a consistent tangent stiffness is obtained as follows: The linearized equations are finally summarized as follows: where Generally, two approaches have been suggested to solve the phase-field-displacement system of equations: simultaneously solving and (monolithic algorithm) and sequentially solving and as coupled staggered fields. However, it is known that monolithic schemes poorly perform in solving Eq. 63, since the energy functional is non-convex with respect to the unknowns [11]. This non-convexity causes the Jacobian matrix in the Newton method to become indefinite, preventing convergence and robustness in monolithic solutions. On the other hand, although staggered algorithms are relatively robust, the time increment must be adequately small to prevent deviating from the equilibrium path. To make a trade-off between robustness and efficiency, a monolithic solution scheme based on the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm has been proposed in the literature to solve the system of coupled equations [34,64]. In this work, the BFGS method presented in the following section is considered as the primary solution scheme.

The Broyden-Fletcher-Goldfarb-Shanno algorithm
In quasi-Newton methods, the stiffness matrix is replaced by an approximation of the stiffness ̃ . The approximated stiffness matrix satisfies the quasi-Newton equation as follows: where The quasi-Newton method is an algorithmic secant method based on the series of successive approximations to the solution, which finds the root of nonlinear equations using the current and previous iteration steps. In the BFGS quasi-Newton method, the approximated stiffness matrix is updated by its predecessor and a correction matrix of rank 2 It can be shown that a symmetric and positive-definite initial guess ̃ 0 leads to the symmetric and positive-definite updated one [39]. Therefore, since the stiffness matrix in Eq. 63 is not necessarily symmetric and positive-definite, an uncoupled stiffness matrix is taken as the initial guess In the above equation, uu 0 is symmetric and positivedefinite. Also, 0 is symmetric and positive-definite for an appropriate length-scale parameter. Accordingly, the approximated stiffness matrix given by the BFGS algorithm is symmetric and positive-definite. It should be noted that although the off-diagonal inter-field coupling terms have been dropped in Eq. 73, the approximation 72 couples the and fields. This is not the case for Newton-based monolithic algorithms adopted for solving weakly coupled problems.
The BFGS algorithm along with a line search method helps to prevent divergence of equilibrium iterations resulting from the inexact stiffness matrix. Therefore, the solution is updated as follows: where the multiplier s is chosen in the way that component of in the search direction is zero with a tolerance Here, convergence is assumed to be achieved if both the residual and solution correction controls are met, that is where g max is the largest residual in the balance equations for the displacement field ( = u ) and the phase-field ( = ), g is the overall time-averaged flux for the displacement and phase-field obtained during the current time step including the current increment, c max is the largest correction to the field variable given by the current iteration, and z max is the largest incremental change to the corresponding solution variable in the current time increment. Here, the tolerances are taken to be g = 0.01 and c = 0.01 . The implementation of the BFGS algorithm is outlined in Table 1.

Consistent tangent moduli based on the Jaumann-Zaremba stress rate
To integrate the viscoelastic model into the incremental-iterative FE framework, the tangent modulus tensor needs to be explicitly specified. However, a closed-form calculation of the tangent tensor is not a straightforward task. Here, we use an efficient numerical approximation of the tangent moduli proposed by Sun et al. [58]. In this approach, by perturbing the deformation gradient, the tangent moduli for the Jaumann rate of the Cauchy stress are accurately approximated by a forward difference of the Cauchy stresses. The Jaumann rate of the Cauchy stress can be expressed as

The linearized incremental form of Eq. (78) is then obtained from
To numerically calculate components of ℂ J , Eq. (79) is perturbed by applying small perturbations to components of D and W tensors. Here, W ij and D ij tensors with perturbed (i,j) components are expressed as and where the corresponding perturbed F ij is obtained from perturbing its (i,j) component as [41] (78) where is a small perturbation parameter. By substituting Eq. (82) into Eqs. (80) and (81), we have It is noteworthy that D has six independent components due to its symmetry. Therefore, the choice of (i,j) would be (1,1), (2,2) (83) (85)

Fracture experiments and numerical simulations
In the following section, the potential of the quasi-Newton method-based monolithic algorithm in achieving convergence and reducing the computation time of the rate-dependent fracture of solids is first investigated. The numerical results of CT tests of BNP/epoxy samples are then compared with the experimental data to identify the material parameters required for the proposed PFM. Next, the effect of hygrothermal conditions and deformation rate on the fracture behavior of the polymer nanocomposites is studied. Finally, the model's capability in predicting fracture patterns is qualitatively assessed using single-edge notched tests.

Experiments
To prepare the nanocomposite specimens, BNPs are first mixed with an epoxy resin. The epoxy system chosen in this study is a commercially available amine-cured epoxy system from Olin Epoxy, namely the epoxy resin of type AIRSTONE 880E and the hardener of type AIRSTONE 886H with an epoxy-hardener mixing ratio of 100:31 by weight. A dispersion process is then conducted to break up nanoparticle agglomerates under mechanical shear loading using a three-roll mill (Exakt, 80E). The dispersion process ensures that nanoparticles are well dispersed in the epoxy matrix. The mixture is next de-gassed at room temperature and poured into a mold for curing nanocomposite plates. The curing process is performed at 353 K for 5 h. The mass density and glass-transition temperature of the epoxy are 1.2 g/cc and 355 K, respectively. The plates are cut into CT specimens with a width of 35 mm and a nominal thickness of 5 mm using a CNC milling machine according to DIN EN ISO 13586. A pre-crack with a length of 16 mm is introduced into the notch root by forcing a razor blade. The pre-crack ensures that a crack is formed ahead of the razor blade tip. The specifications of the specimen are illustrated in Fig. 2. Fracture toughness experiments are carried out at a deformation rate of 10 mm/min and 296 K. Since experimental results may be affected by uncertainties resulting from the manufacturing process, each force-displacement response is obtained from five CT tests to secure the experimental data statistically.

Simulations
We first study the crack growth in a CT specimen. The specimen specifications and dimensions in mm are illustrated in Fig. 2. Considering the symmetry at the mid-length of the specimen, the FE analysis of half of the specimen using symmetric boundary conditions would provide a complete solution of the full model with less computational cost.
Loading and boundary conditions of the reduced model are  Table 2. The model is discretized with four-noded quadrilateral (Q4) elements. Figure 3b shows the corresponding mesh with 6130 Q4 elements. Due to stress and strain concentrations at the bottom-right part of the model, the mesh is refined toward the part, so that the characteristic element size of the fine-level discretization is eight times smaller than the phase-field length scale. In the following CT simulations, two times the vertical displacement of the point placed at the external force is used to measure the displacement.
The resulting force-displacement curves of the CT tests obtained from both the quasi-Newton monolithic and staggered solutions are shown in Fig. 4a. The load is applied to a neat epoxy sample at a constant deformation rate of ̇u = 10 mm/min, = 296 K and w = 0 . In the monolithic solution, the load is applied with a constant displacement increment of u = 5 × 10 −3 mm. It can be seen that the staggered solution is sensitive to the increment size, and two orders of magnitude smaller displacement increment (i.e., u = 5 × 10 −5 mm) is required to reproduce the monolithic solution. The cumulative number of iterations versus displacement is also presented in Fig. 4b. The simulation results show that recovering the monolithic result using a staggered solution requires a displacement increment of 5 × 10 −5 mm and a total number of iterations that is one order of magnitude larger. The significant reduction in the number of iterations justifies the computational efficiency and applicability of the monolithic solution for studying rate-dependent fracture problems in polymerbased materials. Next, CT simulation tests are performed to asses the predictive capability of the proposed PFM model in capturing  [59] Modified Kitagawa parameters (Eqs. (37) and (38)) (K -1 ) 0.01093 [60] ref (  We then study the effect of temperature on the fracture behavior of BNP/epoxy nanocomposites. Figure 6 shows the force-displacement response in CT simulation tests of specimens with 15 wt% of BNPs. Here, the deformation rate is kept constant at ̇u = 10 mm/min, while temperature varies from = 296 to 346 K. The simulation results show that the peak force decreases from around 98 to 51 N by increasing temperature from 296 to 346 K. Furthermore, the displacement at the fracture initiation rises from 0.27 to 0.53 mm by increasing temperature from 296 to 346 K. The simulation results can be explained as follows. First, the shear and bulk modulus associated with the equilibrium, non-equilibrium, and volumetric responses decrease by increasing temperature (see Eqs. (37) and (38)), leading to a less stiff material. Second, the nonlinear viscoelastic flow defined by the Argon model in Eq. (22) is temperature-dependent, so that the viscoelastic flow increases with an increase in temperature. As a result, a larger portion of the strain energy is dissipated at higher temperatures, leading to a smaller peak force and larger displacement at the fracture initiation. It is noteworthy that the effect of temperature on the energy release rate has not been considered in the simulations. Although this assumption may be acceptable in the studied range of temperatures according to experimental data [32], further studies are required to investigate the variation of energy release rate with temperature. In Fig. 7, CT simulation tests are conducted to investigate the effect of the deformation rate ( ̇u ) on the fracture response of BNP (15 wt%)/epoxy nanocomposites. Here, temperature is set to be 346 K, and large deformations up to 10 mm at different deformation rates, varying from 0.001 to 1 mm/ min, are applied to the specimens. From the figure, the displacement at the peak force decreases from 1.31 to 0.53 mm by increasing the deformation rates from 0.001 to 1 mm/ min. In analogy with the effect of temperature, the viscous effect is smaller, and the dashpot behaves more like a solid at higher deformation rates. Therefore, a smaller portion of the strain energy is dissipated, and displacement at the fracture initiation and peak force become smaller. The evolution of damage in a BNP(15 %wt)/epoxy sample at the deformation rate of 0.01 mm/min and imposed displacements of 1, 2, 5, and 10 mm is illustrated in Figs. 8a-d.
In the following, we study the effect of moisture on the fracture behavior of the nanocomposites. Figures 9a and  9b show the force-displacement response in CT simulation tests of specimens with 15 wt% of BNPs at 296 and 346 K, respectively. In the simulations, the deformation rate is kept constant at ̇u = 10 mm/min, while the moisture content varies from 0 to 3 wt%. Here, 3 wt% is considered as the saturation moisture content that can be contained in the polymer nanocomposites. The influence of moisture content on the material behavior of the nanocomposites is captured using the modified Guth-Gold model presented in Eq. (36). The simulation results show that the peak force, respectively, decreases from around 98 to 87 N and from 5 0 to 45 N at 296 and 346 K by increasing the moisture content from 1 to 2 wt%. The peak force, however, rises by further increasing the moisture to 3 wt%. For instance, the peak force increases from around 87 to 91 N at 296 K by increasing the moisture content from 2 to 3wt%. This observation can be interpreted as follows. The epoxy's moduli decrease  with moisture absorption when the moisture content is less than 1.7 wt%. However, when the moisture content exceeds 1.7 wt%, moduli regain some of their reduction as discussed in Eq. (36). Therefore, the peak force increases by increasing the moisture content from 2 to 3wt%.
Finally, the capability of the proposed model to predict fracture patterns is evaluated using the well-known singleedge notched tensile and shear tests of BNP(15 %wt)/epoxy samples. The geometry and boundary conditions are shown in Fig. 10a. A horizontal notch is placed at middle height from the left outer surface to the center of the specimen. The bottom side of the specimen is fixed, while the top side is moved. Both tensile and shear loads are applied at the deformation rate of ̇u = 10 mm/min with constant displacement increments of 10 −4 mm. The simulations are performed at 296 K under plane strain conditions with the material parameters listed in Table 2. Meshes are refined in areas where cracks are expected to propagate. Accordingly, 12509 elements and an effective element size of 0.003 mm in the central strip of the specimen are generated for the tensile test, and 21045 elements with refined meshes in the lower right diagonal strip of the specimen are used for the shear test. Also, the length-scale parameter is set to be l 0 = 0.015 mm. The predicted fracture patterns for the two cases are shown in Figs. 10b, b. It can be seen that the crack path is horizontal for the tensile case, while there is a curved crack path for the pure shear case. The crack patterns are in agreement with those presented in the literature [43]. To study the effect of moisture on the fracture behavior of the nanocomposites, the tensile and shear tests are performed in both dry and wet  Figures 10a, b present the resulting force-displacement curves at the moisture content of 0 and 2 wt%. According to the simulation results, the peak force obtained from the tensile and shear tests, respectively, decreases from around 43 to 38 N and from 18 to 16 N by increasing the moisture content to 2 wt% (Fig. 11).

Summary and conclusion
A phase-field fracture model has been proposed to investigate the effect of hydrothermal effects on the fracture behavior of BNPs/epoxy nanocomposites at finite deformation. To explore the impact of nanoparticle and moisture contents on the rate-, temperature-dependent fracture evolution in the polymer nanocomposites, the PFM has been coupled to a nonlinear viscoelastic constitutive model in a thermodynamically consistent way. For this, the Helmholtz free energy, which describes the nanocomposites' rate, temperature, and moisture-dependent behavior, is additively decomposed into an equilibrium, a non-equilibrium, and a volumetric contribution with positive/negative components. Within this framework, modified versions of the Guth-Gold and Kitagawa models have been adopted to capture the role played by hydrothermal conditions on the fracture evolution in the materials. To improve the the computational efficiency of the PFM, a quasi-Newton monolithic algorithm proposed in [34,64] has been employed to solve the coupled governing equations. The BFGS algorithm applied to nonlinear viscoelastic fracture problems confirms its higher efficiency in comparison with a staggered solution scheme. Benchmark examples show that recovering the monolithic solution using the staggered solution requires about ten times the cumulative number of iterations is about ten times more number of iterations.
To further evaluate the capability of the proposed PFM, numerical predictions should be compared with the experimental data at different hydrothermal conditions and strain rates in the future. Also, the effect of nonuniform dispersion of nanoparticles and moisture content and varying temperature profiles across a specimen on the fracture behavior of polymer nanocomposites needs to be investigated in future studies. Furthermore, interactions between water molecules, nanoparticles, and an epoxy matrix would cause changes in the material properties such as the viscosity and energy release rate. However, due to the complex interactions at small scales, the mechanisms leading to these possible variations are not clear. To gain a deep understanding of the microstructure's effect on the macroscopic properties, the phase-field modeling can be coupled to molecular models to characterize the polymer nanocomposites' material behavior [4,5]. Furthermore, nanoparticles tend to form agglomerates in an epoxy matrix, resulting in insufficient dispersal [33]. It leads to degrading material properties due to relatively inferior interfacial interactions between nanoparticles and the matrix. The effect of surface modification of BNPs on the the fracture behavior of the polymer nanocomposites should be investigated in the future studies.