Finite Element Implementation of the Exponential Drucker–Prager Plasticity Model for Adhesive Joints

Nomenclature

1  Introduction

Structural adhesives have been widely used in the aerospace, aviation, shipbuilding, and automotive industries due to a variety of benefits, including multi-material bonding, good load distribution across bonded areas, less environmental degradation, and a simple manufacturing process. It is also relatively inexpensive when compared to other joining methods [1,2]. To effectively use adhesives, it is necessary to forecast the behavior of adhesive joints under complex loading combinations. As a polymer, adhesives’ mechanical behavior and failure rely on hydrostatic pressure, strain rate, and temperature [3]. In typical applications, hydrostatic pressure is the most visible factor influencing the mechanical behavior of polymers. It has a significant impact in polymer yield strength, plastic flow, strain hardening, and failure plasticity [4]. The effect of hydrostatic pressure on material yield strength can be explained by the Drucker-Prager yield criterion, which is a modified version of the von Mises criterion as a function of hydrostatic pressure. The model was originally a linear function and was developed to predict the failure of soil [5]. The exponential version of this criterion shows a good prediction of polymer base material behavior [6]. Despite the development of yield functions in previous works [7], the flow function for plastic deformation is limited to a linear function, resulting in a non-associated flow when using the exponential Drucker-Prager yield criterion. With this type of flow, the simulations are prone to diverge when subjected to a serious tri-axial load in tension or negative hydrostatic pressure, which is a critical load case for adhesives [8].

This article presents the numerical implementation of the exponential Drucker-Parger plasticity model in the commercial finite element software ABAQUS through the user subroutine UMAT. The development of subroutines was restricted to elastoplastic behavior, and for yield surface modification, only hydrostatic pressure was considered. The associated flow during plastic deformation was expected to increase the numerical stability for adhesive joint simulations. The experimental findings from specialized testing using the modified Arcan fixture [9], which is specifically intended to produce a distinct state of hydrostatic pressure inside an adhesive layer, were utilized to validate the developed model. Experiments and simulations using the second exponential Drucker-Parger model correlated well.

2  Material and Behavior Model

From the classical plasticity theory, the total strain increment for an elastoplastic model can be decomposed into elastic and plastic strain increments. For 3D problems, they are represented by 2nd tensor forms as follows:

Δε=Δεel+Δεpl(1)

2.1 Elastic Behavior

Based on Hooke’s law, the incrementation of stress-strain relationship in the elastic deformation is written as follows:

Δσ=DΔε(2)

2.2 Flow Rule

A general form of flow rules for plastic behavior is written as:

dεpl=dλ∂g∂σ(3)

The associated flow is achieved when the plastic potential function is the same as the yield function (g = f). For this flow, the plastic flow direction will always be normal to the yield surface. Fig. 1 depicts an associated where g = f and non-associated flow where g ≠ f. The plastic flow direction is perpendicular to the function g in the case of non-associated flow.

Figure 1: Concept of the associated flow rule and non-associated flow rule

The flow rules in commercial finite element software are commonly non-associated flows since they can be used for a variety of purposes. In this case, an implicit return mapping algorithm (backward Euler method) is used for high precision calculations into the program [10]. In some load cases, the difference between a plastic potential function and a yield function causes the calculation to diverge. Eq. (4) is the plastic potential function for the exponential Drucker-Prager model in ABAQUS software, which is also non-associated flow.

g=(ξσytan⁡β)2+(σeff)2−ptan⁡β(4)

The hyperbolic function is the only available option for the plastic potential function in ABAQUS when using a general exponent yield criterion. The 2nd order exponential Drucker-Prager yield criteria (parabolic function) is recommended even though the order of exponential Drucker-Prager yield criteria can be arbitrary [11]. In the next section, this non-associated flow from ABAQUS will be compared to the model that was made with an associated flow.

2.3 Implicit Integration

The numerical integration method is required in order to solve the elastoplastic differential equations. The common methods used in computational processes are the forward Euler method, the backward Euler method, and the midpoint method. While the forward Euler method is the simplest computational process, it has a limit on computational stability and high error responses. The backward Euler method is more complex than the forward Euler method but highly accurate for the increase of stress that affects the yield surface expansion. The method is also very stable in terms of computation [12]. For this reason, the backward Euler method is very popular for solving elastoplastic. Its equation is written as in Eq. (5) and Fig. 2 illustrates its approximation step.

Figure 2: Numerical integration of the elastoplastic equations

σ|k+1=tσ+Δσ−DΔεpl|k+1(5)

With 0σ=tσ+Δσ being fixed over iterations, the Eq. (5) can be written as the Eq. (6).

σ|k+1=0σ−DΔεpl|k+1(6)

The residual vector is then defined as follows:

R|k+1=σ|k+1−0σ+DΔεpl|k+1≅ [ 0 ]6×1(7)

According to the associated flow rule, the plastic strain increment can be defined as follows:

Δεpl|k+1=Δλk+1∂f∂σ|k+1=Δλk+1a|k+1(8)

where ∆λ is the plastic multiplier and a is the flow vector. Substituting the flow vector into Eq. (7).

R|k+1=σ|k+1−0σ+Δλk+1Da|k+1(9)

Then, applying the Newton-Raphson method in the residual vector (Eq. (9)).

R|k+1=R|k+∂R∂σT|kδσ|k+∂R∂Δλ|kδλk(10)

Thus

∂R∂σT=∂σ∂σT−∂0σ∂σT+ΔλD∂a∂σT=[I]6×6+ΔλDA=Q(11)

where ∂σ∂σT=[I]6×6, ∂0σ∂σT=[ 0 ]6×6 and ∂a∂σT=∂2f∂σT∂σ=A

A is a second order derivativities tensor of yield function, thus

∂R∂Δλ=∂(ΔλDa)∂Δλ=Da(12)

Substituting Eqs. (11) and (12) into Eq. (10), then the increased stress can be finally determined by rearranging Eqs. (13) to (14).

R|k+1=R|k+Q|kδσ|k+Da|kδλk≈0(13)

δσ|k=−Q−1|kR|k−Q−1|kDa|kδλk(14)

2.4 Yielding Function

The exponential Drucker-Prager models are expressed as follows:

f(σ|k+1,Δλk+1)=a(σeffk+1)b−pk+1−ptk+1(15)

where σeff is von Mises stress, p is hydrostatic pressure, pt is hardening constant, a and b are material parameters that are independent to plastic deformation.

σeff=12((σ11−σ22)2+(σ22−σ33)2+(σ33−σ11)2+6(σ122+σ132+σ232))(16)

p=−13(σ11+σ22+σ33)(17)

The von Mises stress and hydrostatic pressure are defined as Eqs. (16) and (17). The strain hardening for the Drucker–Prager model is defined by the multilinear method in plastic shear hardening where hydrostatic pressure is zero [13]. Then, Eq. (15) can be rewritten as Eq. (18).

pt=aσeffb(18)

The plastic shear stress equivalent plastic strain relationship in the hardening rule is defined as follows:

ptk+1=pt0+Hpt(εeqpl)k+1=pt0+Hpttεeqpl+HptΔλk+1(19)

The Newton-Raphson method is then applied to yielding function in Eq. (15).

f(σ|k+1,Δλk+1)=f(σ|k,Δλk)+f(σ,Δλ)∂σT|kδσ|k+∂f(σ,Δλ)∂Δλ|kδλk(20)

Defining f(σ,Δλ)∂σT|k=aT|k and ∂f(σ,Δλ)∂Δλ|k=−Hpt, substituting into Eq. (20).

f(σ|k,Δλk)+aT|k(−Q−1|kR|k−Q−1|kDa|kδλk)−Hδλk=0(21)

By rearranging Eqs. (21) to (22), the plastic multiplier increment can be finally determined.

δλk=f(σ|k,Δλk)−aT|kQ−1|kR|k(aT|kQ−1|kDa|k+Hpt)(22)

Eqs. (14) and (22) will be used in the following section for finite element implementation.

3  Mechanical Properties of Adhesive

The structural adhesive used in this study is an epoxy-based adhesive called SikaPower-497, manufactured by Sika. The mechanical properties of adhesive were identified from a modified-Arcan test using a universal testing machine (Instron 5567). The specimen deformation was followed by an image correlation system (ARAMIS, manufactured by GOM) as shown in Fig. 3.

Figure 3: Modified-Arcan test [7]

The modified Arcan fixture was designed to induce a different state of stress within an adhesive layer by varying its angle. A direction of 0° and 90° represents a triaxial and pure shear mode. A direction of 30° and 60° induces a mixed state of stress, which is a tensile-shear mode in an adhesive joint. A direction of 120° represents a compressive-shear mode. The experimental yield and failure surfaces are plotted on the von Mises stress-hydrostatic pressure axis as shown in Fig. 4. They both depend on hydrostatic pressure with the non-linearity that corresponds to the exponential Drucker-Prager model. The elastic properties and yield surface parameters are summarized in Tables 1 and 2, respectively.

Figure 4: The Drucker-Prager yield and failure surface [7]

The parameters of the nonlinear hardening equation (Eq. (23)) were identified using simulations of the entire modified Arcan fixture to consider the deflection of the fixture. They are summarized in Table 3.

σeff=σy+q(1−e−cεeqp)+Hεeqp(23)

4  Finite Element Implementations

The exponential Drucker-Prager model with an associated flow rule was implemented in ABAQUS software via a user subroutine, UMAT. With this developed model, the high-order exponential Drucker-Prager is also possible. It also provides higher simulation stability at high-order exponential and requires fewer model parameters since the yield and plastic potential function are the same (associated flow). In subroutine UMAT, elastic and plastic behavior have to be established from input data (material properties). For elastic behavior, it is the same as in ABAQUS. Concerning plastic behavior, an iterative calculation approach is needed to find the solution of a stress tensor increment for plastic behavior (δσ) and the plastic multiplier increment (δλ). At the end of each increment, the subroutine UMAT returns the updated values of stress tensor (σ) and equivalent plastic strain (εeqpl) to the program for the next increment calculation. The flow chart of user subroutine UMAT can be summarized in Fig. 5.

Figure 5: Flow chart of user subroutine UMAT

4.1 Validation of Model

The finite element simulations of the modified Arcan test have been carried out as a validation of the implemented Drucker-Prager model. The four cases were performed: (1) an ABAQUS material model with 2nd order exponential Drucker-Prager yield function and non-associated flow rule (Eq. (4)); (2) a developed user subroutine UMAT with 2nd order exponential Drucker-Prager yield function and associated flow rule; (3) an ABAQUS material model with 4th order exponential Drucker-Prager yield function and non-associated flow rule (Eq. (4)), and (4) a developed user subroutine UMAT with 4th order exponential Drucker-Prager yield function and associated flow rule. The 4th order in the cases was selected because this high order fits relatively well the experimental data of the initial yield surface (Fig. 6). All mechanical properties for the simulations are summarized in Tables 1 to 5.

Figure 6: The plastic potential function

The modified Arcan fixture and adherends are made of steel (Young’s modulus = 200 GPa, Poisson’s ratio = 0.3 and yield stress = 570 MPa). The detailed geometries and boundary conditions are presented in Fig. 7. The adhesive layer has a thickness of 0.3 mm with a total cross-sectional area of 70 × 10 mm2. Because the modified-Arcan fixture is symmetrical in the x-y plane, the finite element model was only constructed in half, and the fixed boundary condition in the z-direction (U3 = 0) across the entire surface is applied. The tie condition is used for the constraint between the fixture and the adhesive layer. The model has 97,860 elements in total, with 350 elements for the adhesive layer. The element volume ranges from 0.3 to 15.5 mm3 and is kept constant at 0.3 mm3 for the adhesive layer. All fixture components, including an adhesive layer, were assigned a 3D 8-node solid element (C3D8) as shown in Fig. 8.

Figure 7: The detailed geometries and boundary conditions

Figure 8: The element size of the modified Arcan test in direction 0°

4.2 Results and Discussion

The results of stress field on a modified Arcan fixture loaded in directions of 0°, 60°, 90°, and 120° are shown in Fig. 9. The deformations of the fixture are also comparable to the experimental test and cannot be neglected. This confirms the necessity of the entire fixture simulation.

Figure 9: The modified Arcan simulation in direction: (a) 0°, (b) 60°, (c) 90°, and (d) 120°

The comparison between experiments and simulations is shown in Fig. 10 for directions of 0°, 60°, 90°, and 120°, respectively. They are the plots of average stress and average strain as the stress-strain fields from a modified Arcan test are not uniform across the section. The simulation results for directions of 0°, 60°, and 90° from cases 1 and 2 are in good agreement with experimental results, while in cases 3 and 4, the discrepancy is more obvious. However, none of the cases can reproduce the results in the direction of 120°. Considering the comparisons with ABAQUS, these confirm the validity of the developed user subroutine UMAT. They also show that the epoxy-based adhesive SikaPower-497 behavior is compatible with the 2nd order exponential Drucker-Prager model.

Figure 10: The modified Arcan simulation result in direction: (a) 0°, (b) 60°, (c) 90°, and (d) 120°

From the simulations, the ABAQUS material models with 2nd and 4th order exponential Drucker-Prager yield function and non-associated flow rule give almost identical results as the developed UMAT with 2nd and 4th order exponential Drucker-Prager yield function and associated flow rule. For further investigation, higher order exponential model simulations were performed for 120° loading. At the 9th order of exponential function (Fig. 11), the result does not converge for the ABAQUS material model. However, with the developed user subroutine UMAT, a converged result can be obtained despite a long computational time. This implies that the associated flow provides more numerical stability than the non-associated flow, especially for a very high order of exponential function and at highly positive hydrostatic pressure loads. In common application, there is no significant difference between the associated and non-associated. However, the associated flow still has the advantage of having fewer material parameters since it uses the same function to describe the yield surface and plastic potential function.

Figure 11: The simulation result of 9th order exponential model

5  Conclusions

The finite element implementation of the exponential Drucker-Prager model for the elastoplastic material has been carried out and validated experimentally with epoxy-based adhesive joints. The numerical implementation of the ABAQUS software was accomplished using a user subroutine, UMAT. This subroutine defines the plastic behavior of the adhesive joint using a yield and plastic potential function. With an associated flow, these functions are identical, providing the advantages of fewer plastic flow parameters and better numerical stability at very high orders of exponential function compared to when a non-associated flow is used. The developed UMAT is also capable of simulating high order exponential function of the Drucker-Prager model which makes it more flexible to apply to other elastoplastic materials.

For the adhesives used in this study (SikaPower-497), its hydrostatic pressure dependent yield surface and plastic behavior have been identified using the modified Arcan tests. The experimental results from different directions of the modified Arcan fixture are used as a validation for the developed UMAT. The simulation result from the 2nd order exponential Drucker-Prager model shows a good agreement with experiments. For the 4th order exponential Drucker-Prager model, it predicts an initial yield surface better than the 2nd order model but poorly interprets the plastic deformation of the adhesive. Concerning the results at direction 120°, neither the ABAQUS material model nor the developed user subroutine UMAT can reproduce the experimental results. This divergence is also expected since the experimental results at direction 120° show the mix of cohesive and adhesive failure modes. In this direction, the compressive stress induced by the fixture strengthens the adhesive joint, so its strength is close to the interface strength leading to the mixed failure mode. The interface element option in ABAQUS [14–16] was not used in this work since the developed UMAT focuses only on material behavior. This limitation opens the way to further development of this UMAT. Another important factor that affects the simulation results is the hardening rule [17–19]. Its parameters can also be further investigated based on hydrostatic pressure and effective plastic strain to provide better simulation results.

Funding Statement: This research was funded by King Mongkut’s University of Technology North Bangkok. Contract No. KMUTNB-PHD-62-07.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.