A nonlinear viscoelastic–viscoplastic constitutive model for adhesives under creep

In this study, we consider the nonlinear viscoelastic–viscoplastic behavior of adhesive films in scarf joints. We develop a three-dimensional nonlinear model, which combines a nonlinear viscoelastic model with a viscoplastic model using the von Mises yield criterion and nonlinear kinematic hardening. We implement an iterative scheme for the viscoplastic solution and a numerical algorithm with stress correction for the combined viscoelastic–viscoplastic model into finite element analysis. The viscoelastic component of the model is calibrated using creep-recovery data from adhesive films in scarf joints. The viscoplastic parameters are calibrated from the residual strains of recovered creep tests with varying load durations. A two-dimensional form of the model shows good agreement with the three-dimensional model for the scarf joint considered in this work and is compared with experiment. The numerical results show favorable agreement with the experimental creep and recovery responses of two epoxy adhesive systems. We also discuss the contribution of nonlinear viscoelasticity and viscoplasticity to the stress/strain distribution along the adhesive center lines. Viscoplasticity tends to lower the stress concentration.


Introduction
In recent years, the use of adhesively bonded joints in engineering applications has increased, especially in aerospace industries. Compared with traditional joining methods, adhesively bonded joints lower stress concentrations, reduce weight, and improve resistance to corrosion. Adhesive joining can also complicate disassembly, require careful surface preparation, and exhibit temporal and plastic responses that are not well understood (Ali and Wahab 2015). Understanding the mechanical behavior of the adhesive within a bonded joint is essential to predict joint response and strength. However, constraints from the adherends and the unbonded interfacial region between the adhesive and the adherends limit the ability of adhesives in bulk to describe thin film adhesive response (Hart-Smith 1981;Afendi et al. 2011).
As most adhesives are polymeric materials, they exhibit time dependence and possibly nonlinearity and plasticity under high stress and high temperature. There have been a variety of models incorporated into the finite element method (FEM) to predict the stress/strain states under monotonic and uniaxial loading. Some have used Schapery's single integral model to describe the creep behavior and stress distribution within adhesives Henriksen 1984;Roy and Reddy 1986). Popelar and Liechti modified the free volume approach of Knauss and Emri Knauss (1981), Knauss et al. (1987) by considering the effects of distortion (Popelar et al. 1997) and hydrostatic pressure (Popelar and Liechti 2003) to predict the nonlinear load-displacement response under monotonic shear tests for adhesives in Arcan specimens. The models and experiments were in good agreement. Others have extended rheological models formed with Maxwell and Kevin-Voight elements to nonlinear viscoelasticity to model the creep of adhesives (Majda and Skrodzewicz 2009;Choi and Reda Taha 2013;Houhou et al. 2014;Badulescu et al. 2015).
The foregoing nonlinear viscoelastic work concerns materials with deformation that is fully recoverable. Some polymers exhibit viscous flow, resulting in permanent deformation, which can also be described using viscoelastic principles (Silva et al. 2017;Acha et al. 2007). For other polymers, the unrecoverable deformation occurs above a yield point (which can be time-or rate-dependent Hu et al. 1996), for which several viscoplastic models have been developed. Perzyna's viscoplastic model is often used, given its ability to accommodate a variety of yield criteria (Perzyna 1966). Pandey and Narasimhan (2001) developed a threedimensional elasto-viscoplastic model using Perzyna's model with a pressure-sensitive modified von Mises yield function for the stress distribution in lap joints. Su and Mackie (1993) studied the influence of creep on the thick adherend shear tests in which a two-dimensional elasto-viscoplastic model consisting of rheological elements was developed and calibrated with creep tests on bulk adhesives. Morin et al. (2015) and G'sell and Jonas (1979) developed an elasto-viscoplastic model based on the G'sell plastic model to simulate the stress response to uniaxial compressive loading with a wide range of strain rates for a structural adhesive. The results agreed with experiments.
Some materials exhibit both viscoelastic and viscoplastic behavior, which has received less attention in the literature than the responses separately. Groth (1990) studied the stress distribution in single lap joints using three different rheological viscoelastic-viscoplastic models. Rocha et al. (2019) proposed a similar rheological viscoelastic-viscoplastic-damage model that accurately described rate-dependent plasticity of an epoxy resin. Dufour et al. (2016) presented a fully coupled viscoelastic-viscoplastic damage model in which the generalized Maxwell viscoelastic model was integrated with a viscoplastic model using Raghava's yield criterion. The model was validated on a series of tensile monotonic experiments on Arcan joints at different loading angles. Krairi and Doghri (2014) proposed a thermodynamically based viscoelastic-viscoplastic model coupling a linear Boltzmann's integral, nonlinear hardening, and ductile damage evolution to describe the uniaxial behavior of polymers. Ilioni et al. (2018) proposed a model combining viscoelasticity with viscoplasticity considering a hydrostatic pressure-dependent yield function for the nonlinear behavior of an adhesive in an Arcan joint. The authors found that viscoplasticity had a strong influence on both shear and tensile monotonic behavior. However, there has been less work on combined viscoelastic-viscoplastic behavior involving nonlinear viscoelasticity. Frank and Brockman (2001) used molecular cooperativity with a Schapery model, incorporating the free volume approach, combined with a modified Bodner-Partom viscoplastic model (Bodner andPartom 1975, 1972) to capture the effects of pressure, strain rate and hardening on multiaxial responses for a glassy polymer. The approach found good agreement with experiment.
The foregoing has described the limited nonlinear viscoelastic-viscoplastic models describing adhesive response. This work combines a modified nonlinear viscoelastic Schapery model with Perzyna's viscoplastic model (using the von Mises yield criterion with nonlinear kinematic hardening). The constitutive model is solved numerically by an iterative algorithm with stress corrections. The model is compared to experimental results of creep tests on scarf joints with adhesives of differing degrees of toughness.

Model description for adhesives
Nonlinear viscoelastic and time-dependent plastic strain has been observed during a series of stress-controlled experiments of adhesive bulk tensile coupons Smith 2016, 2018). At low stress levels, these adhesives exhibited only viscoelastic response with full recovery, whereas at high stress levels, both viscoelastic and viscoplastic response was observed. Accordingly, in the following, the total strain of an adhesive was decomposed into fully recoverable viscoelastic and unrecoverable viscoplastic components with the assumption that they are uncoupled: where the superscripts tot, e, and p denote the total, viscoelastic, and viscoplastic components, respectively.

Viscoelastic component
Schapery's single integral model for nonlinear viscoelasticity, which is based on the concepts of irreversible thermodynamics, was employed to describe the recoverable timedependent responses as (Schapery 1969) 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 (temperature effects were not considered here), and D 0 is the instantaneous compliance. The transient compliance D ψt can be expressed in the form of a power-law (Tuttle and Brinson 1986;Shuangyin and Tsai 1994) or Prony series (De Prony 1795) and characterized by creep-recovery experiments (Gamby and Blugeon 1989). In this study, D ψt was expressed as a five-term Prony-series Lai and Bakker (1996) generalized the uniaxial Schapery's model to three dimensions for isotropic materials with the assumption that the hydrostatic and deviatoric components were fully uncoupled. The modified formulation is presented as ε e,t ij = e e,t ij + 1 3 ε e,t kk δ ij where δ ij is the Kronecker delta, e e,t ij , ε e,t kk , and S t ij are the deviatoric strain, volumetric strain, and deviatoric stress at the current time t , respectively.
The instantaneous shear compliance J 0 and bulk compliance B 0 are evaluated from where v is the Poisson ratio assumed to be constant. Similarly, the transient shear and bulk compliance are determined by Substituting Eqs. (6) and (7) into Eq. (5) and employing a hereditary integral proposed by Haj- Ali and Muliana (2004), we obtain the deviatoric strain e e,t ij and the volumetric strain ε e,t kk in the recursive form as where q t− t ij,n and q t− t kk,n are the shear and volumetric hereditary integrals at the previous time increment (t − t ) for each term in the Prony series. They are updated and stored at the end of the current time increment as

Viscoplastic component
The evolution of the viscoplastic strain rate was given by Perzyna's model (Perzyna 1966) with the associated flow ruleε where η is a viscosity parameter, (f ) is the overstress function expressed with the yield function f , and · is the McCauley bracket such that where σ 0 y is the initial yield stress, and N is a rate sensitive constant. From Eq. (12) we can see that the viscoplastic component is activated when the yield function is larger than zero.
The von Mises yield criterion with nonlinear kinematic hardening was employed for the yield function f since it has been shown to agree with the plastic response of the adhesives considered in this work (Mohapatra 2018(Mohapatra , 2014. The function f is written as where α t ij is the deviatoric backstress tensor at the current time. The evolution law of the backstress is defined from the Armstrong-Frederick model (Frederick and Armstrong 2007) where c is the initial kinematic hardening modulus, κ determines the rate at which the kinematic hardening modulus decreases with the increasing plastic deformation, andε p,t e is the rate of the effective viscoplastic strain, 3 Numerical algorithm

Discretization
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 . Accordingly, the total strain from Eq. (1) can be presented as After some algebraic manipulation, the viscoelastic deviatoric and volumetric strain increments from Eqs. (8) and (9) can be written as Assuming that the time increment t is small, we obtain the following relationship between the nonlinear viscoelastic parameters: The approximate stress increment at the beginning of each time increment may be obtained by substituting Eq. (21) into Eqs. (19) and (20): Using the backward Euler integration method, the viscoplastic strain from Eq. (12) can be represented with the current increment as Similarly, the current backstress tensor is transferred into an incremental form

Viscoplastic formulation
From Eq. (16) we can see that to calculate the effective viscoplastic strain ε p,t e , the value of ∂f ∂σ t ij is needed. The derivative is presented as The yield function f depends on not only the current stress σ t ij but also on the backstress α t ij . Simultaneously, α t ij is a function of σ t ij and ε p,t e , whereas ε p,t e depends on f . The interactions between each variable do not allow ∂f ∂σ t ij and ε p,t e to be found directly. Newton's iteration was employed to calculate the viscoplastic strain according to the flowchart shown in Fig. 1. At the beginning of the iteration, we initialize by where γ is a constant. The increment of the effective viscoplastic strain ε p,t e is expressed with γ as Accordingly, the derivative of the effective strain increment with respect to the stress is given by The derivative of the deviatoric backstress with respect to the stress is obtained from Eq. (25) whereas the derivative of the deviatoric stress with respect to the stress is expressed as Substituting Eqs. (29)-(31) into Eq. (26), we obtain ∂f ∂σ t ij as a function of the current stress and backstress tensors as with the initial value of γ , the difference is defined as To minimize R, allowing the iteration converge, γ is corrected at the end of every iteration according to The superscripts K and K + 1 present the Kth and (K + 1)th iterations, respectively. It should be noted that the number 4 in Eq. (34) was used to slow down the rate of convergence, preventing γ from being negative without changing the iteration result. Substituting γ K+1 into Eq. (28), the (K + 1)th loop begins. Once the iteration converges, namely R < 10 −7 , the viscoplastic strain increment ε p,t ij can be calculated from Eq. (24).

Viscoelastic-viscoplastic formulation
The viscoelastic-viscoplastic model was implemented into FE using a UMAT in ABAQUS 2016. The stress and the Jacobian matrix are required for each time increment in the UMAT. Figure 2 shows the iterative algorithm for the combined viscoelastic-viscoplastic system. At the beginning of the nth time increment t n , we assume the adhesive behaves viscoelastically. A trial stress increment accounting only for viscoelasticity, σ tr , was calculated from Eqs. (22) and (23) with the nonlinear parameters and variables from the previous time history t − t n . Based on the trial stress, those parameters and the increment of the viscoelastic strain ε e, t ij were updated. To implement viscoplasticity, the backstress α t ij and the yield function f were obtained by the trial stress. The iteration of viscoplastic strain (Fig. 1) was activated when f was larger than zero.

Fig. 2 Flowchart of the viscoelastic-viscoplastic system
To decide if the calculation has converged, a strain comparison was introduced by where the total strain increment ε tot,t ij is obtained as a combination of both components, and ε ext,t ij is calculated from the stiffness matrix and external force and then passed into the UMAT. The stiffness matrix is updated for each increment from the Jacobian in the UMAT.
The root of the squared differences is defined by When ξ * < 10 −7 , the iteration is assumed to be converged. If the convergence is not reached, then the stress is corrected by where the superscripts K and K + 1 represent the Kth and (K + 1)th iteration, respectively. In Eq. (38), The derivative of the viscoelastic strain increment to the stress, ∂ ε e,t ij ∂σ t kl , is calculated according to the condition of the effective stress exceeding the linear viscoelastic range. The effective stress σ t is given as For nonlinear viscoelasticity, the derivative can be obtained from Eqs. (19) and (20) as In this study the nonlinear parameters are defined as linear functions of the current effective stress σ t : where a i and b i are constants. For linear viscoelasticity, Therefore Eq. (41) becomes The derivative of the viscoplatic strain increment is obtained from Eq. (24) as In Eq. (45) the high-order derivative is given by where (54) Fig. 3 Geometry of the studied scarf joint. All dimensions are in mm

Plane strain environment
The 3D constitutive model was used to simulate the creep-recovery behavior of a scarf joint as shown in Fig. 3. Given the high computational cost of the 3D model, the suitability of a 2D model was considered. For a shear dominated stress state, as occurs in the adhesive of the scarf joint considered here, plane stress and strain formulations provide similar solutions. Nevertheless, the relatively stiff aluminum adherends create a nearly plane strain condition in the adhesive and were used here. The stress and strain components were reduced according to

Mesh and boundary conditions
The scarf joint consisted of two adherends and an adhesive film with a thickness of 0.21 mm. The finite element meshes for the 2D and 3D models were created using the preprocessor in ABAQUS 2016 as shown in Fig. 4. In the 2D model ( Fig. 4-a) the adherends and adhesive were mostly meshed with 4-node plane strain bilinear reduced integration with hourglass control elements (CPE4R). At the near surface-transition zone from the adherend to adhesive 3-node linear plane strain triangle elements were used for the wedge shapes. To capture the deformation within the adhesive layer, the adhesive was divided into four equal layers through the thickness. Whereas a convergence study showed minimal improvement from two layers in the adhesive, the refinement in the relatively small adhesive layer had a negligible computational cost. Similarly, the 3D model ( Fig. 4-b) was meshed with 8-node linear brick, reduced integration with hourglass control elements (C3D8R), with 6-node triangular prism elements (C3D6) for the transition zone. The surface on one end of the scarf joint was constrained in all directions and rotations. The surface on the other end was allowed to deform in the longitudinal direction of the joint (x-axis). Load was applied as a pressure over the surface of the free end. Accordingly, the average shear stress in the adhesive was obtained by τ adh 12 = sin 10 • cos 10 • σ joi xx , where τ adh 12 is the shear stress in the adhesive, and σ joi xx is the applied pressure.

Parameter calibration for viscoelasticity
This work considered the response of a "standard" (Cytec FM300-2) and "toughened" (Hysol EA9696) film adhesive in a scarf joint of which the adherends were made of 2024 aluminum. The adherends were modeled elastically with a modulus of 73 GPa and Poisson's ratio of 0.33. Both adhesives were thermoset and epoxy-based. A Poisson's ratio of 0.38 and 0.43 (Mohapatra 2018) was used for the standard and toughened adhesive, respectively. The model was calibrated from 10000 s creep tests of scarf joints at 20% and 50% of the ultimate shear strength (USS). Considering the linear viscoelastic response with a constant shear stress σ adh 12 = 20% USS, Eq. (5) reduces to from which D o was found at t ≈ 0, whereas the Prony series parameters were determined at t 0.
For nonlinear creep, Eq. (5) reduces to ε e 12 = 1 2 g 0 J 0 σ adh 12 + g 1 g 2 σ adh whereas the associated recovery is The nonlinear parameters g 0 , g 1 , g 2 , and a were fit to linear functions of the effective stress by Eq. (42). To provide continuity between the linear and nonlinear viscoelastic response, g 1 for both adhesives was represented as where a 1 , b 1 , and c 1 are constant. Accordingly, the instantaneous strain at 50% USS was used to obtain g 0 from Eq. (59). The recovery strain was used to find g 2 and a at 50% of USS from Eq. (60). Finally, the creep response was used to obtain g 1 using Eq. (59). It should be noted that the permanent strain was subtracted from the above experimental responses at 50% USS before fitting the viscoelastic parameters, which are listed in Table 1.

Time-dependent plasticity
The viscoplastic parameters c and k were determined from monotonic scarf joint tests at 4448 N/min. Assuming that the adhesive was under pure shear, Eqs. (14) and (15) were combined and rewritten as It should be noted that c and k describe the development of the yield surface, whereas η and N in Eq. (12) influence the time-dependent plasticity (viscoplasticity). To calibrate these two parameters, the effect of viscoplasticity in the adhesive films was considered.
A series of creep tests at 50% USS with different durations were conducted on the scarf joints. The applied load was ramped to 50% USS over 50 s for the toughened and standard adhesives and held constant until the load was removed at double the loading rate, after which the scarf joints were allowed to recover. The average adhesive strain was measured    Table 2.

Numerical analysis
The calibration described above considered only one load duration for viscoelasticity (10,000 s) and only the permanent deformation for viscoplasticity. Using the calibrated parameters above, the predicted total shear response of the adhesive films to a load duration of 10000 s are compared with experiment for creep-recovery tests of scarf joints in Fig. 6. The experiments at 20% and 50% USS were repeated three times for the toughened and standard adhesive. As shown in Fig. 6, the nonlinear viscoelastic-viscoplastic (NEP) model showed good agreement with experiments at 20% and 50% USS for both adhesives. To show the contribution of nonlinear viscoelasticity, dotted lines in Fig. 6 describe linear-viscoelastic (LEP) response. Neglecting nonlinear viscoelasticity led to underestimating total strain for the toughened by 21% and overestimating strain for the standard adhesive by 4%. The standard adhesive experienced nonlinear stiffening, resulting in a decreasing instantaneous compliance (g 0 < 1 in Eq. (2)) in the nonlinear viscoelastic-viscoplastic model. Thus the linear viscoelastic-viscoplastic model, with a constant instantaneous compliance, overestimated the strain response. The 2D and 3D models, compared in Fig. 6, were within 1% of each other. Given the similarity of the results, and the computational efficiency of the 2D model, it was used through-  Figure 7 compares the predicted creep-recovery responses of varying load durations at 50% USS with experiment. The predicted strain response was within the experimental repeatability for creep and recovery. The predicted recovery followed the experimental response for each duration. The residual strain corresponding to the viscoplastic deformation developed nonlinearly and reached 7.2% and 3% of the total strain up to the load duration of 10000 s for the toughened and standard adhesives, respectively. Figure 8 compares the distributions of shear stress and strain along the center line of the adhesives at 50% USS for three cases: NEP, LEP, and nonlinear viscoelastic (NE). The NEP model shows that time reduces the shear stress concentration but magnifies the strain concentration, as observed elsewhere (Su and Mackie 1993). For both adhesives ( Fig. 8(a, c)), the stress distribution predicted by the NE model changed little from 1000 to 10000 s, showing that viscoplasticity accounted for the reduced stress concentration at longer durations. In Fig. 8(b, d) the NE model produced a smaller strain concentration than the models incorporating viscoplasticity; this shows that viscoplastic response magnifies the strain concentrations in these adhesives. Note that the higher strain of the standard adhesive from NE (in comparison to the other models) is a result of nonlinear hardening for this adhesive (Fig. 8d).

Conclusion
In this study, we proposed a three-dimensional nonlinear viscoelastic-viscoplastic model combining the modified Schaper integral model and Perzyna viscoplastic model with the von Mises yield criterion and nonlinear kinematic hardening, an iterative scheme assuming only viscoelasticity initially with a correction of stress in both viscoelasticity and viscoplasticity to solve the combined system numerically. The nonlinear viscoelastic and viscoplastic parameters were calibrated from creep-recovery and monotonic tests on scarf joints with a toughened and a standard adhesive film. The existence and development of viscoplasticity was discussed by comparing the viscoplastic component of the model with residual strain from creep tests on the scarf joints. The viscoplastic strain increased nonlinearly with increasing load duration. A two-dimensional model considering plane strain was developed and compared with the 3D model in simulating the creep-recovery of the adhesives in scarf joints. A good agreement between the prediction and the experimental results was obtained. The finite element analysis showed that viscoplasticity reduced the stress concentration and magnified the strain concentration.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.