A nonlinear viscoelastic–viscoplastic model for adhesives

We consider the nonlinear viscoelastic–viscoplastic behavior of adhesives. We develop a one-dimensional nonlinear model by combining Schapery’s nonlinear single integral model and Perzyna’s viscoplastic model. The viscoplastic strain was solved iteratively using the von Mises yield criterion and nonlinear kinematic hardening. The combined viscoelastic–viscoplastic model was solved using Newton’s iteration and implemented into a finite element model. The model was calibrated using creep-recovery data from bulk adhesives and verified from the cyclic behavior of the bulk adhesives. The finite element model results agreed with experimental creep and cyclic responses, including recoverable and permanent strain after load removal. Although the contribution of the viscoplastic strain was small, both viscoplastic and viscoelastic components of strain response were required to describe the adhesive creep and cyclic response.


Introduction
Adhesively bonded joints are widely employed in the aerospace and aviation industries due to their advantages of light weight, stress distribution over a wider area, reduction of stress concentration, and increased resistance to fatigue and corrosion, as compared with conventional joining methods. With the increasing use of adhesive joints, understanding their resistance to damage, fatigue, and mechanical behavior is important. As most adhesives are polymeric, their mechanical behaviors are time-dependent and may involve nonlinearity and viscoplasticity (Groth 1990;Tuttle et al. 1994;Pasricha et al. 1996;Cognard et al. 2006;Lemme 2016;Morin et al. 2015). Since adhesive joints often involve complex geometry and loading, the finite element method (FEM) is a common tool in their design. However, the application of these analyses faces two challenges: (1) the accuracy of the constitutive model used to describe the physical phenomena and (2) the discretization and solution algorithms.
There have been several studies using finite element analysis to describe the timedependent behavior of adhesive systems, among which Schapery's general single integral constitutive law (Schapery 1969), from the concepts of irreversible thermodynamics, is commonly used for nonlinear viscoelastic materials. Henriksen (1984) extended Schapery's model for multiaxial stress and implemented it using finite elements. The Poisson effect was incorporated in the model, and the strain was decomposed into an instantaneous and a hereditary component which was updated by a recursive method. Results from this approach found good agreement with experimental creep data in uniaxial tension. Following Henriksen's work, Roy and Reddy (1986) developed a finite element formulation of the multiaxial Schapery model based on an updated Lagrangian formulation of motion. The computational model was validated from analytical and experimental work (Peretz and Weitsman 1983;Romanko et al. 1982) on adhesive joints. The model was extended to account for moisture effects (Roy and Reddy 1987). Lai and Bakker (1996) modified Henriksen's model to incorporate the effects of temperature and physical aging through the definition of a reduced time. They allowed different tensile and compressible response by separating the strain into deviatoric and hydrostatic components and defining nonlinear viscoelastic parameters. Henriksen (1984) and Lai and Bakker (1996) used a similar recursive algorithm with nonlinear viscoelastic parameters, which limited the time step size. Based on their work, Haj-Ali and Muliana (2004) developed a new iterative procedure with a stress correction algorithm to reach the converged stress and strain at each time step.
Perzyna's viscoplastic model (Perzyna 1971) is often used for time-dependent plastic deformation of polymers and adhesives. Corigliano and Ricci (2001) applied Perzyna's viscoplastic model to interfaces for time-dependent delamination processes in composite materials. Giambanco and Fileccia Scimemi (2006) modified the model for decohesion processes in adhesively bonded joints and proposed a numerical implementation. Following Perzyna's work, Pandey et al. (1999), Pandey and Narasimhan (2001) developed a uniaxial elastoviscoplastic model using Perzyna's model with a pressure-sensitive modified von Mises yield function. They discretized the equilibrium equation into an incremental form for finite element implementation. Nguyen et al. (2016) combined the linear generalized Maxwell viscoelastic model and Perzyna's model incorporating a modified Drucker-Prager yield function and kinematic hardening. They also introduced a pressure-sensitive failure criterion, which captures the rate-dependent mechanical behavior and failure of amorphous glassy polymers. Whereas the foregoing work considered linear viscoelasticity and viscoplasticity, many adhesives are nonlinear, particularly at high stress levels and temperatures. To our knowledge, there is no model that combines nonlinear viscoelasticity with Perzyna's model for adhesives.
The cumulative deformation under cyclic loading for polymers and adhesives is known to be complicated and important for failure prediction. It can be modeled using nonlinear viscoelastic models as demonstrated by Xia et al. (2005) and Tscharnuter et al. (2012) and using elastoviscoplastic models as proposed by Motta et al. (2019), Yu et al. (2001), Hassan et al. (2011), andAsmaz (2015). However, the cyclic behavior of adhesives involving both viscoelasticity and viscoplasticity has received less attention, particularly for solutions involving finite elements (Tuttle et al. 1994;Nguyen et al. 2016;Rocha et al. 2019). We consider the nonlinear viscoelastic and viscoplastic response of a standard and toughened adhesive under uniaxial creep and cyclic loading. The approach combines Schapery's single integral viscoelastic model and Perzyna's viscoplastic model. An iterative solution for the viscoplastic strain and an algorithm for the viscoelastic-viscoplastic system are proposed and implemented using FEM. The model is compared with experimental results of uniaxial creep and cyclic tests.

Nonlinear viscoelastic-viscoplastic model
In this study, we decompose the total strain of an adhesive into viscoelastic and viscoplastic components with the assumption that the viscoelastic and viscoplastic components are uncoupled (Frank and Brockman 2001;Dufour et al. 2016;Ilioni et al. 2018): where the superscripts tot, ve, and vp represent total, viscoelastic, and viscoplastic components, respectively. The viscoelastic strain is recoverable, whereas the viscoplastic strain is unrecoverable.

Viscoelastic component
Schapery's single-integral model (Schapery 1969) was employed to describe the viscoelastic strain as where ψ t is the effective time given by a, g 0 , g 1, and g 2 are nonlinear parameters dependent on the current stress σ t and temperature. (For linear viscoelastic materials, g t 0 = g t 1 = g t 2 = a t = 1.) Temperature effects were not considered. D 0 is the instantaneous compliance, whereas D ψt is a transient compliance expressed as a seven-term Prony series or where D n and λ n are material coefficients. Substituting Eq. (4) into Eq.
(2) yields where q t n is the hereditary integral

Viscoplastic component
The viscoplastic component was accounted for using Perzyna's model (Perzyna 1971) as it has flexibility to accommodate different yield functions. The model in tensor form iṡ where N is a material constant, η is a viscosity parameter, g is the potential energy function, and an associated flow rule, that is, g = f , was used; · is the McCauley bracket defined as where σ 0 y is the initial yield stress, and f is the yield function. The von Mises criterion with nonlinear kinematic hardening was considered as it was shown to agree with the plastic response of adhesives (Mohapatra 2018(Mohapatra , 2014. The yield function is defined as where σ e is a modified effective stress, S ij is the deviatoric stress tensor, and α ij is the deviatoric back-stress. For uniaxial tension, Substituting Eqs. (10) and (11) into Eq. (9), we obtain an expression for the yield surface: The evolution law for the back-stress α ij is defined from the Armstrong-Frederick model (Frederick and Armstrong 2007) where c is the initial kinematic hardening modulus, k determines the rate at which the kinematic hardening modulus decreases with increasing plastic deformation, and ε vp e is the effective plastic strain. Equation (13) can be written for uniaxial tension as By integrating Eq. (14) and substituting into Eq. (12) the back stress can be expressed as The related yield function f is written as 3 Numerical implementation and parameter calibration

Increments and derivatives
In this study, we consider uniaxial stress using the methodology presented above. Considering a time increment t , the strain history ε t at time t is obtained from the prior strain ε t− t and the change in strain over the time increment ε t as Accordingly, Eq.
(1) can be presented as and stress can be expressed as Haj-Ali and Muliana (2004) proposed a numerical integration scheme to obtain the viscoelastic strain increment for Schapery's model, in which ε ve,t can be expressed as where Assuming that the time increment t is small, we obtain Then the stress increment can be obtained by substituting Eq. (23) into Eq. (20): Employing the backward Euler integration method (Butcher 2016), the uniaxial viscoplastic strain from Eq. (7) can be presented in incremental form as From Eqs. (16) and (26) we can see that ε vp,t e and f are interactive. To calculate ∂f ∂σ t and, consequently, obtain ε vp,t e , we employed Newton's iteration (Abadue 1986). A flowchart of the viscoplastic strain calculation is shown in Fig. 1. At the beginning of the iteration, we initialize by Substituting Eq. (27) into Eq. (26), we can write the viscoplastic strain increment as Then the derivative of ε vp,t e with respect to σ t can be obtained from Eq. (28): The derivative ∂f ∂σ t can be obtained from Eq. (16): Combining Eqs. (29) and (30), we have The difference between ∂f ∂σ t and γ 0 is defined as To minimize R, the γ 0 is corrected at the end of every iteration according to The superscripts K and K + 1 represent the Kth and (K + 1)th iterations, respectively. The derivative of R with respect to γ is derived from Eq. (32): Once the iteration converged, namely R < 10 −7 , we get γ and the viscoplastic strain increment from Eq. (28).

Solution algorithm
The viscoelastic-viscoplastic model was implemented into commercial FEA software (Abaqus 2016) with user's material subroutine (UMAT), where the stress and the Jacobian matrix are required for each time increment. As shown before, both viscoelastic and viscoplastic strain increments, as well as the nonlinear viscoelastic parameters, are related to the current stress. Thus, Newton's iteration was employed for the stress calculation. The flowchart is shown in Fig. 2. At the beginning of each time increment t , we assumed a trial stress increment σ tr containing only the viscoelastic stress from Eq. (24) with the parameters passed from the previous increment. The nonlinear viscoelastic parameters and the viscoelastic strain increment ε ve,t were recalculated from the trial stress. To correct the trial stress for subsequent steps, the derivative of the viscoelastic strain increment with respect to the stress was needed, which is described further.
For linear viscoelasticity, For nonlinear viscoelasticity, g t 0 , g t 1 , g t 2 , and a t are functions of σ t . In addition, ψ t = t a t . Thus Eq. where Substituting the trial stress into Eq. (16) to obtain the yield function f , the iteration of the viscoplastic strain (Fig. 1) is activated when f > 0. Considering the stress correction, the derivative of the viscoplastic strain increment with respect to the stress is obtained from Eq. (25): Substituting Eq. (31) into Eq. (37) and neglecting the high-order term ∂ 2 f ∂σ 2 , we obtain Note that when the applied stress is below the yield stress, that is, f ≤ 0, The total strain increment is obtained from both components given by To obtain the current stress and the viscoelastic-viscoplastic strain increment, we used an iterative algorithm with a strain difference. The strain difference was defined as where ε ext,t is calculated from the stiffness matrix and external force and passed into the UMAT (Abaqus 2016). The stiffness matrix is updated at each increment from the Jacobian in the UMAT. The convergence was defined by ξ < 10 −7 . If convergence has not occurred, then stress correlation is needed. The stress at the (K + 1)th iteration was updated by where the superscripts K and K + 1 represent the Kth and (K + 1)th iterations, respectively. In Eq. (42), which is also the inverse of the Jacobian matrix.

Parameter calibration
In this work, we considered the responses of a "standard" (Cytec FM300-2) and "toughened" (Hysol EA9696) film adhesive, both of which are thermoset and epoxy-based. Nonlinear viscoelasticity has been observed for both systems above 50% UTS, whereas permanent deformation occurred above 80% of UTS (Lemme 2016). Considering a linear viscoelastic response with constant stress σ 0 , Eq. (5) reduces to whereas the associated recovery is The nonlinear parameters g 0 , g 1 , g 2, and a were fit to linear functions of stress of the form Accordingly, the instantaneous strain at 50% and 80% of UTS was used to obtain g 0 from Eq. (45). The recovery strain was used to find g 2 and a at 50% and 80% of UTS from Eq. (46). Finally, the creep response was used to obtain g 1 using Eq. (45). It should be noted that permanent strain was subtracted from the above experimental response at 80% UTS before fitting the viscoelastic parameters listed in Table 1.
Referring to the viscoplastic component parameters c and k, they were determined from uniaxial tension tests (Mohapatra 2018) using Eq. (16) such that The material constants η and N in Eq. (7) were fit to the creep-recovery response at 80% of UTS. The viscoplastic parameters are listed in Table 2.

Comparison of the measured and predicted creep and recovery response
Using the UMAT with Newton-Raphson iteration, the nonlinear model was compared with experiment (three replicates for each case) for creep-recovery of two adhesives at three levels of stress. The FE model was 152-mm long with 100 two-node truss elements (T2D2), where one end of the truss was fixed, and the other end was under a concentrated force. A comparison of the simulated and experimental results is shown in Fig. 3. The viscoplastic component was activated and combined with the viscoelastic component at 80% UTS, whereas the responses to other two levels of stress are only viscoelastic (full strain recovery).
Good agreement of the model with the experimental data for the toughened adhesive was obtained, both in creep and recovery. Good agreement was also observed for creep with the standard adhesive, as well as the viscoelastic recovery. The predicted viscoplastic strain at 80% UTS was 14% lower than experiment. The difference amounted to 70 µε, which is within the experimental repeatability of the extensometer for this duration (30 h).
The contributions of the viscoelastic strain and viscoplastic strain at 80% UTS are compared in Fig. 4. The figure shows that both components increase with time, whereas only the viscoelastic strain recovered after load removal. The viscoplastic strain is unrecoverable as expected, contributing 4.9% and 2.3% to the adhesive strain for the toughened and standard adhesives, respectively.

Cyclic behavior and discussion
To verify the model capacity of predicting the time-dependent response of adhesives, we compared the model with the cyclic behavior of bulk adhesives. The coupons were loaded with a saw tooth wave at 0.5 Hz in uniaxial tension. The load was cycled between the maximum and 10% of the maximum load for one thousand and ten thousand cycles, after which the load was removed, and the coupons were allowed to recover. The different cyclic durations allowed comparison of the viscoplastic component during recovery. The maximum load was 20%, 50%, and 80% of UTS. The peak strain of each cycle was captured during the loading stage. More detail about the experiments can be found elsewhere (Lemme 2016). Figure 5 compares the peak cyclic strain between the simulated and experimental results. The agreement is generally favorable. The permanent strain increased as the cycles increased from 1000 to 10,000. This demonstrated the viscoplastic nature of both adhesives and particularly the toughened adhesive. The model has good agreement with experiment for 1000 cycles for the toughened adhesive for both cyclic load and recovery. At longer durations (Fig. 5b) and 80% UTS the model overestimated experiment, with a difference of 140 µε at 10,000 cycles (which is within the experimental accuracy of the strain measurement for these durations). For the standard adhesive, the coupons loaded to 80% UTS broke after 1000 cycles, preventing comparison of its recovery response. Comparison of the remaining results from the standard adhesive are, nevertheless, favorable, where the maximum difference was within 4% of the experiment.
The contribution of the viscoelastic and viscoplastic behaviors may be compared from the energy dissipated during cyclic loading. This was done for 1000 cycles at 80% UTS for both adhesives. The viscoelastic energy dissipation per cycle was nearly constant with increasing number of cycles (0.8% and 0.5% of the stored energy for the toughened and standard adhesives, respectively). Although the viscoplastic energy dissipation decreased slightly with increasing cycles, it was more than two orders of magnitude smaller than the viscoelastic energy dissipation for both adhesives. However, the cumulative effects of viscoplasticity are not necessarily small. The cumulative viscoplastic dissipated energy reached 2% of the stored energy over the 1000 cycles considered here.

Conclusion
In this study, we proposed, calibrated, and verified a one-dimensional nonlinear viscoelasticviscoplastic computational model. Schapery's single-integral model with a Prony series was employed for the viscoelastic component, whereas Perzyna's model with the von Mises criterion and nonlinear kinematic hardening was employed for the viscoplastic component. The viscoplastic strain was calculated using the Newton iterative method. A solution algorithm with stress correction for the combined viscoelastic-viscoplastic model was developed. The model showed good agreement with experimental uniaxial creep and cyclic response, where the adhesives exhibited nonlinear viscoelastic and viscoplastic behaviors. Although the adhesive response was dominated by viscoelastic behavior, viscoplastic strain was needed to describe it accurately. The computational model provides a framework for both adhesive characterization and numerical implementation and may be extended to multidimensional analysis.