An invariant-based anisotropic material model for short fiber-reinforced thermoplastics: Coupled thermo-plastic formulation

Accepted Manuscript An invariant-based anisotropic material model for short fiber-reinforced thermoplastics: coupled thermo-plastic formulation A. Dean, J. Reinoso, S. Sahraee, R. Rolfes PII: DOI: Reference:

S1359-835X(16)30194-4 http://dx.doi.org/10.1016/j.compositesa.2016.06.015 JCOMA 4341

To appear in:

Composites: Part A

Received Date: Revised Date: Accepted Date:

26 September 2015 31 March 2016 18 June 2016

Please cite this article as: Dean, A., Reinoso, J., Sahraee, S., Rolfes, R., An invariant-based anisotropic material model for short fiber-reinforced thermoplastics: coupled thermo-plastic formulation, Composites: Part A (2016), doi: http://dx.doi.org/10.1016/j.compositesa.2016.06.015

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

An invariant-based anisotropic material model for short fiber-reinforced thermoplastics: coupled thermo-plastic formulation A. Deana , J. Reinosob,∗, S. Sahraeea , R. Rolfesa a Institute b Elasticity

of Structural Analysis. Leibniz Universit¨ at Hannover, Appelstr. 9A, 30167 Hannover, Germany and Strength of Materials Group, School of Engineering. University of Seville, Camino de los Descubrimientos s/n, 41092, Seville, Spain

Abstract This paper presents the development of a new fully-coupled thermo-mechanical invariant-based elastoplastic constitutive model for short fiber reinforced composites (SFRPs). The invariant-based character of the current formulation under quasi-static multiaxial loading conditions allows the incorporation of the anisotropic (transversely isotropic) response of these materials, which results from the employed injection molding process for their production. From the modeling point of view, the thermo-mechanical nonlinear behavior of these materials is performed through the definition of a thermal-dependent anisotropic yield surface and a non-associative plastic potential function. This novel coupled formulation is accordingly derived and implemented into the FE code FEAP by means of the user capability UEL, i.e. the constitutive model is integrated within an user-defined element. The performance of the model is first assessed using a set of benchmark computations, and subsequently validated with available experimental data, showing the potential applicability of the proposed formulation. Keywords: A. Finite element method (FEM); B. Short fiber reinforced thermoplastics; C. Thermo-mechanical coupling; D. Transversely isotropic thermo-plasticity

1. Introduction The achievement of current demands with regard to the design and manufacturing of optimized structural components is gaining a notable relevance in different engineering applications, especially in the automotive and aeronautical industries. One of the most innovative aspects regards the gradual incorporation of composite materials, substituting traditionally employed metals due to their excellent specific strength and stiffness ratios. In the particular case of the automotive industry, where very high production rates are required, this is the case of short fiber reinforced plastics (SFRPs) for the realization of intricate design concepts using injection moulding production techniques. Aiming at fully exploiting the mechanical capabilities in terms of strength and stiffness of SFRPs, the characterization of the mechanical response of such materials is required. This is principally motivated by the fact that the behaviour of these materials strongly depends upon several factors such as the internal microstructural arrangement, i.e. fibre alignments at different locations within the component, the environmental conditions such as temperature, aging agents, moisture, among others. Regarding the mechanical response of glass fiber reinforced poly-amide PA6GF-x or PA66GF-x (where x stands for the fiber content of the specimen), several investigations have analyzed the effect of the fiber orientation through FE-based modeling [22, 2] or using different experimental techniques, see [7, 36] and the references given therein. Experimentally, a great deal of research has been devoted to the characterization ∗ Corresponding

authors Email address: [email protected] (J. Reinoso)

Preprint submitted to ———

June 20, 2016

of the mechanical response of SFRPs subjected to different loading actions [10, 15, 16], fatigue response [20, 25, 4], and fracture behavior [34, 3], among other aspects. At present, many engineering procedures involving short fiber reinforced thermoplastics and their posterior service applications are performed at different temperatures [6]. In line with this, several investigations have characterized the thermal and geometric dependency of the material short glass fiber reinforced components under quasi-static [21, 12, 29], rate [46] and fatigue loading conditions [13] for different fiber orientation angles. Within the context of the production of automotive parts, a novel manufacturing technique is the denominated tempered hybrid clinching procedure, which combines metallic and composite entities to reduce the structural weight of the final part. Basically, clinching is a mechanical interlocking process, where the joining takes place due to the action of an external punch [5], and, in contrast to the conventional joining methods such as welding, brazing or riveting, this technique does not require special surface treatments. This fact makes it an interesting method for highly productive rates for fastening composites and metals. Consequently, the essential challenge of this hybrid clinching technique lies in connecting the materials with different mechanical properties, where the use of temperature could certainly improve the final quality of the resulting joint. Based on the previous arguments, the prospective use of the tempered hybrid metal-composite clinching technique to produce different parts in the automotive industry motivates the development of fully-coupled thermo-mechanical material models, which simultaneously attain the anisotropic character of the composite component (generally PA6GF-x specimens) together with the incorporation of thermal effects in a consistent way. In the related literature, from the numerical standpoint, a significant number of studies have been devoted to the theoretical formulation and numerical treatment of anisotropic plastic models [14, 27, 28, 26, 35, 33, 41] and thermo-mechanical formulations [11, 23, 24, 38, 40]. In this setting, Reese and coauthors developed several few models for the large deformation that accounted for anisotropic inelastic behavior based upon a structural tensorial representation, see [30, 42] and the references therein given. However, the thermo-mechanical anisotropic plasticity has received a limited attention up to now. In this context, most of the numerical investigations with regard to the thermo-mechanical analysis of composites with thermosetting resins have been devoted to the analysis of endless Carbon Fibre Reinforced Polymers (CFRPs) through finite element formulations, see [31, 32], among others. In this paper, a novel invariant-based transversely isotropic thermo-plastic material model, which displays the temperature dependent response of short fiber reinforced thermoplastics is developed. This model relies on the fundamental constitutive hypothesis outlined in [44], but extending the modeling concept to short fiber reinforced composites [43]. In particular, this extension addresses the fully thermo-mechanical coupling in a thermodynamically consistent framework, and also accounts for the effect of plasticity along the fiber direction using an enriched invariant set. The manuscript is organized as follows. Section 2 outlines the thermodynamically consistent formulation of the coupled thermo-mechanically transversely isotropic elasto-plastic model. The description of the algorithmic treatment of the proposed coupled thermo-mechanical elasto-plastic model at material point level within the framework of nonlinear FEM is outlined in Section 3. Section 4 presents the central governing equations and the corresponding FE formulation of the developed formulation. In Section 5, the applicability of the present model is examined through, first, a uniaxial tensile application of a specimen along the fiber direction, and subsequently by means of the experimental-numerical correlation at different temperatures and fiber orientations with regard to the loading direction. Finally, the main conclusions of this investigation are addressed in Section 6. 2. Constitutive formulation This section is devoted to the constitutive formulation of the anisotropic invariant-based model for short fiber-reinforced polymers (SFRPs) for coupled thermo-mechanical analysis. As was previously stated, the anisotropic response of such materials is basically motivated by the mould flow production process, see Figure 1. 2

Through such production process, CT scan analyses have shown that three pseudo-layers with regard to the preferential fiber orientation can be identified [45]: (i) two outer layers which are mainly aligned with the flow direction and (ii) a central layer whose principal fiber orientation is transverse in-plane to the flow direction. These directions represent the preferential material orientations. To simplify this complex internal arrangement [1], an average procedure over the thickness of the specimen can be performed, which approximately leads to an average fiber alignment of around 65% with respect to the mold flow direction and 35% transverse to that direction. At the meso-mechanical level, from a mathematical perspective, fibers orientations can be attained by means of a tensor-based representation, where the three dimensional fibre orientations is defined at material point level using the second order tensor A [8, 37]. This structural operator can be computed by applying the fiber distribution tensor to a group of n fibers in a limited range. Therefore, the structural tensor A can be defined as A := a ⊗ a, (1)

orientation vector , fiber distribution

[-]

where a is the vector that identifies the preferential fibers orientation (the direction with the highest aligned fiber content, which coincides with the molding direction) within the specimen. Based on the representation given in Eq.(1), the material response is invariant with respect to arbitrary rotations around a, as well as to reflections at planes parallel to the fibers and with respect to the planes with normal coincident with a [8, 37]. This set of operations defines the group of symmetries for transversely isotropic material response, which is herewith adopted for SFRPs. Then, the corresponding theoretical formulation, namely the freeenergy function, the yield function and the plastic potential function, which feature the material response, incorporate the structural operator A as argument for its definition. 1 0.9 90

0.8

3

0.7

1

0.6

2

0.5

4

5 0 3

0.4

2

1

0.3 0.2

-90

n ctio inje

0.1

n

ctio dire

at d pl ishe i-fin 6GF30 m e s A of P

e

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

thickness [mm]

Figure 1: Micro-Computed tomography of PA6GF30 plate using injection moulding procedure.

2.1. Coupled thermo-mechanical transversely isotropic free energy definition Within the framework of infinitesimal deformation setting, the total strain tensor (ε), which is defined as the symmetric part of the spatial displacement gradient (ε := ∇s u, where u identifies the displacement field), is additively decomposed into elastic (εe ) and plastic (εp ) counterparts as follows: ε = εe + εp .

(2)

An analogous decomposition scheme is adopted taken for the total entropy of the system (η), leading to the following representation: η = ηe + ηp , (3) 3

where η e and η p are the elastic and plastic entropy components, respectively. The additive decomposition given in Eqs.(2) and (3) is based upon the modeling strategy proposed in [38], which provided an alternative version of the second law of the thermodynamics in its Clausius-Plank form (Section2.2). The principal feature that motivated the adoption of such modeling approach relied on its straightforward application to general thermo-mechanical applications and its robustness from a numerical standpoint. In addition to the previous arguments, one appealing aspect regards the decomposition of the internal dissipation into purely mechanical and thermal counterparts. Let Ψ to denote the Helmholtz free energy function. According to the previous decompositions Eqs.(2)(3), the postulation of Ψ is assumed to be a sum of the mechanical (Ψmec (εe , α)), thermal (Ψth (ϑ)), and thermo-mechanical coupling (Ψthm (εe , α, ϑ)) counterparts, which can be expressed as: Ψ(εe , α, ϑ, A) = Ψmec (εe , α, A)+Ψth (ϑ)+Ψthm (εe , α, ϑ) = Ψe (εe , A)+Ψhard (α)+Ψth (ϑ)+Ψthm (εe , α, ϑ), (4) where α and ϑ stand for the set of internal variables that account for hardening and the temperature of the system, respectively; Ψhard (α) is the hardening part of the free-energy function due to plastic effects. It is worth mentioning that in line with Simo and Miehe [38], the free energy function given in Eq.(4) endowed the thermal effects into the corresponding coupling term Ψthm (εe , α, ϑ). This assumption is merely adopted for simplicity reasons. However, alternative free energy functions which include the thermal dependency of the elastic properties can be also defined and incorporated into the current modeling framework. Based on Eq.(4), it can be observed that the mechanical part can be decomposed into a elastic free energy term Ψe (εe , A) and a hardening term Ψhard (α): Ψmec (εe , α, A) = Ψe (εe , A) + Ψhard (α)

(5)

In what follows, the dependency of the free energy function with respect to the structural tensor A is omitted to alleviate the notation. A particular form of Eq.(4) for coupled thermo-plastic material models reads:   1 e ϑ e e e e e Ψ(ε , α, ϑ) = ε : C : ε − Z : C : ε [ϑ − ϑ0 ] + cp (ϑ − ϑ0 ) − ϑ log + Ψhard (α), (6) 2 ϑ0 where ϑ0 is the reference temperature satisfying Ψ(0, 0, ϑ0 ) = 0, with 0 denoting the second-order zero tensor; Ce identifies the elastic constitutive tensor, which is particularized to transversely isotropic material response in the sequel (see Eq.(8)); Z denotes the second order (transversely isotropic) thermal expansion tensor; and cp represents the specific heat capacity. Based on standard arguments [9], the Cauchy stress tensor (σ) is obtained as the first derivative of the free energy function with respect to the elastic strain tensor, while the elastic constitutive operator (Ce ) is defined as the second derivative of the free energy with respect to elastic strain tensor. These relations can be expressed as: σ := ∂εe Ψ; C := ∂εe εe Ψ. (7) For transversely isotropic materials, the elasticity operator adopts the form: Ce := ∂εe εe Ψe = λ1 ⊗ 1 + 2µT I + α(1 ⊗ A + A ⊗ 1) + 2(µL − µT )IA + βA ⊗ A,

(8)

where I stands for the fourth-order identity tensor, whereas IA = Aim Ijmkl + Ajm Imikl ; λ, α, β, µT and µL identify for the elastic constants, whose definitions and relationships with regard to the engineering constants are given in Appendix A. 2.2. The Clausius-Plank inequality The so-called Clausius-Plank inequality for internal dissipation (Dint ) in case of coupled thermo-mechanical response takes the form: Dint = σ : ε˙ − e˙ + ϑη˙ ≥ 0, (9) 4

where e is the internal energy, which can be computed through the Legendre transformation [19] as follows: e = Ψ + η e ϑ,

(10)

Recalling the previous definitions Eqs.(9)-(10) and complying with the standard Truesdell and Noll procedure [39], the following constitutive equations stemming from Eq.(6) can be set up: σ := ∂εe Ψ = Ce : [εe − Z [ϑ − ϑ0 ]] η e := −∂ϑ Ψ = cp log

(11)

ϑ + Z : Ce : εe ϑ0

(12)

Γ := −∂α Ψ = −∂α Ψhard

(13)

where Γ denotes the so-called hardening force, and εϑ := Z [ϑ − ϑ0 ] defines the temperature-derived strain field. The heat flux equation is assumed to obey the Fourier law as follows: q = −K · ∇ϑ

(14)

where K is the second-order transversely isotropic heat conductivity tensor, and ∇ϑ denotes the spatial gradient of the temperature. With the previous definitions at hand, incorporating Eqs.(11)-(13) and making use of Eq.(10), the restriction over the internal dissipation to fulfill the second law of the thermodynamics reads: ˙ + ϑη˙ p ≥ 0 Dint = σ : ε˙ p + Γ ? α

(15)

where the operator ? stands for any arbitrary product. Consequently, the local dissipation can be written mech th as the additive decomposition of the mechanical Dint and the thermal parts Dint as follows: mech th Dint = Dint + Dint ;

mech ˙ with Dint = σ : ε˙ p + Γ ? α;

th Dint = ϑη˙ p

(16)

2.3. Yield function The thermo-elastic domain E, assuming the maximum dissipation principle [38], can be defined as: E = {(α, ε¯p , ϑ) | f (σ, A, ε¯p , ϑ) ≤ 0},

(17)

where ε¯p identifies the equivalent plastic strain. The definition of the equivalent plastic strain in the present formulation is given by r r 2 p 2 p p p ε¯ := kε k = ε :ε . (18) 3 3 Note that the previous definition is similar to the standard equivalent plastic strain in isotropic von Mises (J2) plastic models [38]. Following [43], we postulate the decomposition of the stress tensor σ into plasticity inducing stresses σ pind and reaction stresses σ reac : σ := σ pind + σ reac

with

σ pind :=

1 1 (tr[σ] − aσa)1 − (tr[σ] − 3aσa)A, 2 2

(19)

where tr[•] stands for the trace operator of a second-order tensor. The construction of the transversely isotropic yield surface (f (σ, A, ε¯p , ϑ)), which accounts for the plastic deformation along the fiber direction as well as the thermal dependency, reads f (σ, A, ε¯p , ϑ) = ζ1 I1 + ζ2 I2 + ζ3 I3 + ζ4 I32 + ζ5 I4 + ζ6 I42 − 1 ≤ 0 5

with

nf :=

∂f , ∂σ

(20)

with nf being the normal vector to the yield surface, and Ii (i = 1, 4) identifying the stress invariants, which are defined as: 1 1 3 (tr[σ pind ])2 − tr[A(σ pind )2 ]; I2 := tr[A(σ pind )2 ]; I3 := (tr[σ]) − tr[A(σ)2 ]; I4 := tr[Aσ dev ], 2 2 2 (21) where σ dev denotes the deviatoric part of the Cauchy stress tensor. The six temperature dependent parameters ζi (¯ εp , ϑ), (i = 1, 6) and their corresponding invariants represent different thermo-mechanical loading states. The first parameter ζ1 is related to transverse shear loading, while ζ2 accounts for the material performance under in-plane shear loading. The parameters ζ3 and ζ4 are associated with thermo-mechanical loading states transverse to the fiber direction, and lastly, the parameters ζ5 and ζ6 govern yielding along the longitudinal direction with respect to the fiber orientation. An schematic representation of the yield surface in the invariant space is portrayed in Figure 2, where having defined the material parameters that characterize the thermo-mechanical response, an unique representation is obtained. In this graph, the interaction between different loading states based on the invariant representation is shown. I1 :=

Figure 2: Schematic representation of the yield surface in the invariant space.

A compact representation of the yield function takes the form: f (σ, A, ε¯p , ϑ) =

1 σ : K : σ + L : σ − 1 ≤ 0, 2

(22)

where 9 dev K := ζ1 Pind + (ζ2 − ζ1 )Pind ⊗ Adev ; A + 2ζ4 (1 − A) ⊗ (1 − A) + ζ6 A 2

3 L := ζ3 (1 − A) + ζ5 Adev , (23) 2

where Adev is the deviatoric part of the structural tensor A, and the operators Pind and Pind A are defined as: Pind = I −

1 1 3 (1 ⊗ 1) + (A ⊗ 1 + 1⊗A) − (A ⊗ A) ; 2 2 2

ind ind ind Pind A := PAijkl = Aim Pmjkl + Amj Pimkl .

(24)

2.4. Plastic potential function In line with [43, 44], the present model introduces the definition of a non-associative flow rule. Therefore, the evolution of the plastic strains is not governed by the gradient of the yield surface. This leads to the 6

definition of the temperature dependent plastic potential function g = g(σ, A, ϑ). Analogously to the definition of the yield function but omitting the linear terms associated with I3 and I4 , the plastic flow function is assumed to take the form: ∂g g(σ, A, ϑ) = ι1 I1 + ι2 I2 + ι3 I32 + ι4 I42 − 1 ≤ 0 with ng := , (25) ∂σ where ιi (i = 1, 4) are the set of temperature dependent plastic potential parameters that are given as follows: ιi = ιi(0) (ϑ0 ) [1 + wι (ϑ − ϑ0 )] , (26) where ιi(0) (ϑ0 ) denote the plastic potential parameters at reference temperature and wι identify a set of weight functions. The plastic flow potential can be expressed in condensed format as: g(σ, A, ϑ) = M : σ − 1 ≤ 0

(27)

with

9 dev M := ι1 Pind + (ι2 − ι1 ) Pind ⊗ Adev . (28) A + 2ι3 (1 − A) ⊗ (1 − A) + ι4 A 2 Figure 3 shows a representation of the plastic potential in the invariant space, where again the interaction between the different loading states can be observed.

Figure 3: Schematic representation of the plastic potential function in the invariant space.

2.5. Internal variables evolution equations The evolution equations of the internal variables, namely the plastic strains (εp ), the hardening variables (α), and the plastic entropy (η p ) are given by the following rate expressions: ε˙ p = γ

∂g(σ, A, ε¯p ) = γng = γM : σ with ∂σ ∂g(σ, A, ε¯p ) ˙ =γ α ∂Γ

ng = M : σ

(29) (30)

4

η˙ p = γ

∂g(σ, A, ε¯p ) ∂g(σ, A, ε¯p ) X ∂M = γgϑ = γ ∂ϑ ∂M ∂ϑ i=1

(31)

where γ represents so-called plastic multiplier. The partial derivative of the operator M with respect to the temperature in Eq.(31) reads: 4

  ∂ι1 ∂M ∂M X ∂ιi ∂ι2 ∂ι3 9 ∂ι4 = = Pind − Pind + Pind + 2(1 − A) ⊗ (1 − A) + Adev ⊗ Adev . A A ∂ϑ ∂ιi i=1 ∂ϑ ∂ϑ ∂ϑ ∂ϑ 2 ∂ϑ

(32)

Finally, the customary Kuhn-Tucker loading/unloading conditions are given by γ ≥ 0;

f (σ, A, ε¯p ) ≤ 0;

γf (σ, A, ε¯p ) = 0

(33)

whereas the consistency condition renders γ f˙(σ, A, ε¯p ) = 0 7

(34)

2.6. Parameters identification The characterization of the current thermo-plastic model requires the determination of material parameters, which define the thermal dependent yield surface Eq.(20) and plastic potential function Eq.(25). In this section, the procedure followed in this study is briefly explained. The reader is referred to Vogler et al. [44] for a more comprehensive description. 2.6.1. Yield surface parameters To determine the six thermal dependent yield surface parameters ζi (i = 1, 6), for each of the triggering points of the yield locus shown in Figure 2, the yield stress vs. the corresponding plastic strain at different temperatures ϑ has to be obtained. According to the yield function definition Eq.(20), the necessary yielding stresses to obtain the complete identification of the yield surface are identified as follows: (i) uniaxial longitudinal tensile yield stress Yutl (εputl , ϑ), (ii) uniaxial longitudinal compressive yield stress Yucl (εpucl , ϑ), (iii) uniaxial transverse tensile yield stress Yutt (εputt , ϑ), (iv) uniaxial transverse compressive yield stress Yuct (εpuct , ϑ), (v) in-plane shear yield stress Yis (εpis , ϑ), and (vi) transverse shear yield stress Yts (εpts , ϑ). For each of these (six) stress states, the plastic strain is related to the equivalent plastic strain through certain correction factors, εpi = χi · ε¯p . These correction factors χi have a direct relation with the choice of plastic potential parameters ιi (i = 1, 4) [44]. Subsequently, with the previous information at hand, six hardening curves featuring the yield stress vs. equivalent plastic strain response can be prepared, giving the definition of each of the yield stresses Yi = Yi (¯ εp , ϑ) mentioned above. Based on this, after having these six material characterization tests at different temperature ϑ, the explicit determination of the yield function parameters ζi (i = 1, 6) for any thermo-mechanical loading state (represented by ε¯p and ϑ) is simply given by solving a 6 × 6 algebraic system of equations. The solution is obtained by substituting each of these material tests in the yield function in Eq.(20) for each ζi (i = 1, 6)). It is worth noting that for practical objectives and due to the fact that it is difficult to accomplish all of these material characterization tests at any possible temperature ϑ, in the algorithmic treatment with regard the corresponding FE analysis, a linear interpolation between the available data could be invoked. 2.6.2. Plastic potential parameters The reason behind using non-associative flow rule is the determination of a realistic plastic deformation behavior, specifically in terms of the plastic contractlity/dilitancy and the identification of the so-called plastic Poisson coefficients. Therefore, the choice of the plastic potential parameters is a notable issue. To determine the appropriate values of the plastic potential parameters ιi (i = 1, 4), a trial and error optimization process using the so called potential stresses Y s , Y ul and Y ut is used [44], which can be expressed as: 1 1 1 ι1 1 ι4 ι1 = 2 , ι2 = 2 , ι4 = 2 , ι3 = 2 − − . (35) 4 4 Y Y Y Y s

s

ul

ut

It should be noted that every single potential stress is a yield stress-like quantity, which is used to construct the plastic potential surface but does not have a direct physical meaning. Thus, similarly to the characterization procedure outlined for the characterization yield surface, after determining the set of the plastic potential parameters, the correction factors χi have to be calculated for each of the six simple characterization tests previously mentioned. The specific procedure to set up the parameters ιi associated with the potential function replicates that outlined in [44] based on a optimization procedure, and consequently is omitted here for the sake of the conciseness. However, it is worth noticing the thermal dependency of these parameters, so that the procedure has to be accordingly repeated for each service temperature taken into consideration. 3. Numerical integration of the transversely isotropic thermo-plastic model: return mapping algorithm In this section, the numerical treatment regarding the integration of the transversely isotropic thermoplastic model introduced in Section 2 is outlined according to the implicit backward Euler time integration 8

algorithm. The algorithmic treatment of the present model is performed within the context of nonlinear FEM using the standard incremental-iterative Newton-Raphson solution procedure, and is formulated at material point level (that is denoted by x ∈ Rndim , with ndim denoting the Euclidean dimension space). (i) Let us to define the prescribed motion of an arbitrary body B (Figure 4) within time interval [tn , tn+1 ], (i) with t ∈ R+ , where tn identifies the previous converged time step and tn+1 stands for the current prospective time step at the global Newton-Raphson iteration i (see Section 4, where the Initial Boundary Value Problem (i) (IBVP) for thermo-mechanical applications is set up). The time increment is denoted as: ∆t = tn+1 − tn . The superscript i is omitted in the following to simplify the notation at material point level. The strain and temperature rates within the given time step are computed as follows: ε˙ =

εn+1 − εn ; ∆t

ϑn+1 − ϑn ϑ˙ = . ∆t

(36)

Recalling basic ingredients regarding nonlinear FEM analysis, we also assume that the internal variables εpn , ε¯pn , αn and ηnp and the prospective total strains εn+1 and temperature ϑn+1 are given data. The generic thermo-plastic constitutive value problem at material point level is drafted in Algorithm (1). Given: εpn , ε¯pn , αn and ηnp , εn+1 and temperature ϑn+1 ; p Find: εpn+1 , ε¯pn+1 , αn+1 and ηn+1 , at the end of the time interval [tn , tn+1 ]; such that; r e

ε˙ = ε˙ − γn ˙ g;

p

η˙ = γg ˙ ϑ;

ε¯˙p = γ˙

2 kng k 3

(37)

and; γ˙ ≥ 0;

f (σ, A, ε¯p ) ≤ 0;

γf ˙ (σ, A, ε¯p ) = 0

(38)

Algorithm 1: Thermo-plastic constitutive value problem at material point level. Note that in Eq.(37)3 , the equivalent plastic strain serves as hardening parameter. The point of departure for the algorithmic treatment of the model is the adoption of the implicit backward (fully implicit) Euler method to obtain the discrete version of the rate expressions given in Eq.(37) within the interval [tn , tn+1 ]. These discrete expressions read: r 2 p p e e p p εn+1 = εn + ∆ε − γn+1 ng,n+1 ; ηn+1 = ηn + γn+1 gϑ,n+1 ; ε¯n+1 = ε¯n + γn+1 kng,n+1 k. (39) 3 with γn+1 ≥ 0;

fn+1 (σ n+1 , A, ε¯pn+1 ) ≤ 0;

γn+1 fn+1 (σ n+1 , A, ε¯pn+1 ) = 0

(40)

where ∆ε = εn+1 − εn . In case of plastic flow (γn+1 > 0) within the given time interval, the constraints outlined in Eq.(40)2 and Eq.(40)3 leads to the following condition: !

fn+1 (σ n+1 , A, ε¯pn+1 ) = 0.

(41)

Following classical arguments within the context of Computational Plasticity, a two step predictorcorrector algorithm is performed. The first step concerns with the computation of the elastic trial step (denoted with the superscript tr in the sequel) according to the following scheme: e εe,tr n+1 = εn + ∆ε

and 9

ε¯p,tr ¯pn n+1 = ε

(42)

e,trial σ tr n+1 = Ce : εn+1 − Z : Ce [ϑn+1 − ϑ0 ] .

(43)

The corresponding trial yield function can be expressed as f (σ tr ¯pn , ϑn+1 ) = n+1 , A, ε

1 tr tr tr σ : Ktr : σ tr n+1 + L : σ n+1 − 1 ≤ 0. 2 n+1

(44)

Note that in Eq.(44) the operators Ktr and Ltr are function of the trial equivalent plastic strain ε¯p,tr n+1 . Therefore, if the elastic trial state fulfill the condition expressed in Eq.(41), i.e. lies within the elastic domain, the trial elastic state is accepted as the actual solution. Otherwise, the solution has to be obtained via the so-called plastic corrector step. Based on this, the computation of the plastic multiplier γn+1 follows the procedure outlined in Algorithm 2. The basic operators that are required for the computation of Algorithm 2 are given in Appendix B. e 1. Compute εe,tr n+1 = εn+1 + γn+1 ng,n+1 . 2. Add the term Z : [ϑn+1 − ϑ0 ] at both sides of the previous expression and substitute e ng,n+1 = M : σ n+1 → εe,tr n+1 − Z : [ϑn+1 − ϑ0 ] = εn+1 − Z : [ϑn+1 − ϑ0 ] + γn+1 M : σ n+1 .

e e e 3. Compute Ce : εe,tr n+1 − Z : [ϑn+1 − ϑ0 ] = C : εn+1 − Z : [ϑn+1 − ϑ0 ] + γn+1 C : M : σ n+1 . tr e 4. Identify: σ n+1 = σ n+1 + γn+1 C : M : σ n+1 . −1

tr tr 5. Compute σ n+1 = [I + γn+1 Ce : M] : σn+1 = H : σn+1 . 6. Solve the equation to determine the consistency parameter γn+1 via local iterative process (local Newton-Raphson index denoted by the superscript k).

(a) (b) (c) (d) (e)

(k=0)

p,(k=0)

Set k = 0 and initial values (σ n+1 = σ tr ¯n+1 = ε¯pn , γ (k=0) = 0) n+1 , ε p Compute f (k) (σ n+1 , A, ε¯n+1 , ϑn+1 ). IF f (k) (σ n+1 , A, ε¯pn+1 , ϑn+1 ) ≤ TOL GOTO 7, ELSE (k) Set residual for local Newton-Rapshon iteration: Rn+1 = f (k) (σ n+1 , A, ε¯pn+1 , ϑn+1 ). (k) Perform linearization of Rn+1 :  (k)

∂fn+1

ˆ (k) ] = R(k) + ∆γ (k) L[R n+1 n+1

(k) ∂σ n+1

(k)

:

(k)

∂σ n+1

+

(k) ∂γn+1

(k)

∂fn+1

(k) ∂Kn+1

····

∂Kn+1 (k) ∂γn+1

(k)

+

(k)

∂fn+1

:

(k) ∂Ln+1

∂Ln+1 (k) ∂γn+1

= 0.

(k)

(f) Compute ∆γ (k) =

−Rn+1 (k)

∂fn+1 (k)

∂σ n+1 (k+1) γn+1

(k) γn+1

(g) Set = + ∆γ (h) k ← k + 1 GOTO (b)

(k)

∂σ n+1

:

(k)

(k)

∂γn+1

(k)

+

∂fn+1 (k)

∂Kn+1

(k)

····

∂Kn+1 (k)

∂γn+1

(k)

+

∂fn+1 (k)

∂Ln+1

(k)

:

∂Ln+1

.

(k)

∂γn+1

. (k)

p,(k)

p,(k)

p,(k)

p 7. Update the internal variables: σ n+1 = σ n+1 , ε¯pn+1 = ε¯n+1 , ηn+1 = ηn+1 , εpn+1 = εn+1 . 8. Compute algorithmic tangent operator, see Appendix C.2.

Algorithm 2: Plastic corrector step: algorithmic update of the plastic multiplier and the internal variables.

4. Thermo-mechanical balance principles and finite element formulation In this section, the fundamentals of the coupled thermo-mechanical problems are introduced. Section 4.1 briefly outlines the continuum balance principles of linear momentum and energy, whereas the corresponding finite element discretization is described in Section 4.2.

10

4.1. Governing equations of coupled thermo-mechanical problem Within the small deformation setting, let B ∈ Rndim to identify the configuration occupied for by an arbitrary body in the Euclidean space of dimensions ndim , see Figure 4. The position of the material points is given by the vector x ∈ Rndim . The boundary of the solid is denoted by ∂B ∈ Rndim −1 , which is subdivided into the disjointed parts ∂Bt ⊂ ∂B and ∂Bu ⊂ ∂B, with ∂Bt ∪ ∂Bu = ∂B and ∂Bt ∩ ∂Bu = ∅, where ∂Bt and ∂Bu stand for the portions of the body boundary where the Neumann or the Dirichlet boundary conditions are prescribed, respectively. The volume forces are denoted by Fv : B → Rndim , whereas the prescribed tractions and displacements conditions are identified with t : ∂Bt → Rndim −1 and u : ∂Bu → Rndim , respectively. In the sequel, we particularize ndim = 3 for three-dimensional analysis.

Figure 4: Schematic representation of the coupled thermo-mechanical problem of an arbitrary body B.

The local form of the balance of linear momentum reads: div [σ] + Fv = 0 in B,

(45)

where div [•] is the divergence operator. The second central equation is the balance of internal energy, which accounting for the Legendre transformation, Eq.(10), takes the form: mech ϑη˙ e = R − div [q] + Dint ,

(46)

where R identifies a given heat source. Substituting in Eq.(46) the constitutive equation for the elastic entropy, Eq.(12), we end up with the evolution equation of the temperature, which renders: mech cp ϑ˙ + ϑZ : Ce : ε˙ e − R + div [q] − Dint =0

(47)

The weak form of the balance principles, Eqs.(45) and (47), can be constructed through the standard Galerkin procedure as follows: Z Z Z t · δu d∂Ω − Fv · δu dΩ = 0 (48) Gu (u, δu) = σ : δεe dΩ − B

ϑ

Z

G (ϑ, δϑ) = B

B

∂B

Z Z Z Z e e ˙ cp ϑδϑ dΩ+ ϑZ : C : ε˙ δϑ dΩ− Rδϑ dΩ− q·∇δϑ dΩ+ B

B

B

∂B

Z hδϑ d∂Ω− B

mech Dint δϑ dΩ = 0

where n represents the outer normal vector, and with h = q · n representing the surface thermal flux. 11

(49)

4.2. Finite element formulation Let the arbitrary body B to be discretized into ne non-overlapping elements, such that B ≈

n Se

B (e) .

e=1

The kinematic and thermal fields (u, ϑ), their respective variations (δu, δϑ) and their increments (∆u, ∆ϑ) are interpolated by means of the standard finite element approach within the element domain as follows: h

u≈u =

n X

A

h

N dA = Nd; δu ≈ δu =

A=1

ϑ≈

n X

A

h

N δdA = Nδd; ∆u ≈ ∆u =

A=1 n X

ˆ δϑ ≈ ˆ ϑ, N A ϑˆA = N

A=1

n X

n X

N A ∆dA = N∆d

(50)

A=1

ˆ ∆ϑ ≈ ˆ ϑ, N A δ ϑˆA = Nδ

A=1

n X

ˆ ˆ ϑ, N A ∆ϑˆA = N∆

(51)

A=1

where n denotes the node number at the element level; N A (ξ) is the interpolation shape functions with ξ := {ξ 1 , ξ 2 , ξ 3 } ∈ A identifying the coordinates of the isoparamateric domain A := {ξ ∈ R3 | − 1 ≤ ξ i ≤ +1; i = 1, 2, 3}; dA and ϑˆA stand for the nodal displacements and temperature, which are respectively ˆ N represents the matrix containing the standard element shape functions collected into the vectors d and ϑ; ˆ is the matrix for the temperature interpolation within the element. related to the displacements, and N The variation and the increment of the strain field are approximated by means of the displacement-strain operator B as follows: δε ≈ Bδd, ∆ε ≈ B∆d, (52) whereas the spatial gradients of the thermal field are given by ˆ = Bϑ ϑ; ˆ ˆϑ ∇ϑ ≈ ∇ξ N

ˆ = Bϑ δ ϑ; ˆ ˆ ϑ δ (∇ϑ) ≈ ∇ξ Nδ

ˆ = Bϑ ∆ϑ. ˆ ˆ ϑ ∆ (∇ϑ) ≈ ∇ξ N∆

The discrete versions of Eqs.(48)-(49) substituting the Fourier law Eq.(14), can be cast into: Z  Z Z ˆ = δdT ˆ u (d, δd, ϑ) G BT σ dΩ − NT t d∂Ω − NT Fv dΩ = 0 B

B

(54)

B

∂B

Z    e Z e  ˆ δ ϑ) ˆ = δϑ ˆT ˆ T cp ϑn+1 − ϑn dΩ + ˆ T ϑZCe εn+1 − εn dΩ− ˆ ϑ (d, ϑ, G N N ∆t ∆t B B  Z Z Z Z
T ˆ dΩ + ˆ T R dΩ + ˆ T h d∂Ω − ˆ T Dmech dΩ = 0 N Bϑ KBϑ ϑ N N B

(53)

∂B

B

(55)

int

The consistent linearization procedure of Eqs.(54)-(55) is performed through the directional derivative concept [9], adopting the following representation: ˆ ∆ϑ) ˆ =G ˆ + ∆d G ˆ ˆ ∆ϑ)∆ ˆ ϑ ˆ ˆG ˆ u ](d, δd, ∆d, ϑ, ˆ u (d, δd, ϑ) ˆ u (d, δd, ∆d, ϑ)∆d ˆ u (d, δd, ϑ, L[ + ∆ϑˆ G

(56)

ˆ δ ϑ, ˆ ∆ϑ) ˆ =G ˆ δ ϑ) ˆ + ∆d G ˆ δ ϑ)∆d ˆ ˆ δ ϑ, ˆ ∆ϑ)∆ ˆ ϑ ˆ ˆG ˆ ϑ ](d, ∆d, ϑ, ˆ ϑ (d, ϑ, ˆ ϑ (d, ∆d, ϑ, ˆ ϑ (d, ϑ, L[ + ∆ϑˆ G

(57)

where ∆∗ [•] denotes the directional derivative operator with respect to the field ∗. The residual mechanical ˆ and Rϑ = G ˆ δ ϑ). ˆ The stiffness matrices ˆ u (d, δd, ϑ) ˆ ϑ (d, ϑ, and thermal vectors are identified with Ru = G u ˆ ˆ ∆ϑ), ˆ Kϑu = ˆ ˆ u (d, δd, ϑ, of the system are denoted as follows: Kuu = ∆d G (d, δd, ∆d, ϑ), Kuϑ = ∆ϑˆ G ˆ δ ϑ) ˆ , and Kϑϑ = ∆ ˆ G ˆ δ ϑ, ˆ ∆ϑ). ˆ ˆ ϑ (d, ∆d, ϑ, ˆ ϑ (d, ϑ, ∆d G ϑ The fully coupled thermo-mechanical system takes the form:      Kuu Kuϑ ∆d Ru (58) ˆ = Rϑ Kϑu Kϑϑ ∆ϑ where Kuu and Kϑϑ are the mechanical and thermal element stiffness matrices, respectively, whereas the coupling matrices between both fields are denoted by Kuϑ and Kϑu . Details of the operators stemming from the system outlined in Eq.(58) are given in Appendix C.1. This nonlinear set of equations has been implemented into the FE code FEAP through the user-defined capability UEL using a 20-node element to minimize prospective locking effects. 12

5. Applications The objective of this section is to examine the predictive capabilities of the coupled transversely thermoplastic model herewith proposed. First, the performance of the model is verified through a benchmark problem consisting of a uniaxial tensile loading test, where the preferential fiber orientation is aligned with the external solicitation (Section 5.1). The results obtained from the simulations at room temperature are compared with the corresponding experimental data, while the model performance is further investigated numerically at different temperatures. Subsequently, the material model is validated through experimentalnumerical correlation with regard to experimental data taken from literature concerning monotonic tensile tests at different temperatures and with different material orientations (Section 5.2). It is worth mentioning that in the present simulations, an average material orientation vector a over the specimen thickness is assumed [12]. 5.1. Model verification: monotonic tensile test In the first application, the response of a bone-shaped specimen subjected to monotonic tensile loading is investigated. The specimen was manufactured using the short-fiber reinforced material PA6GF-30, whose elastic material properties are listed in Table 1 (within the proposed range addressed in the corresponding ISO norm), while the thermal properties are given in Table 2. Due to the lack of the required material tests with regard to different loading conditions, reasonable assumptions were made to obtain rational simulation results. The plastic potential parameters that characterize the yield function, ζi , at room temperature are listed in Table 3, whilst the plastic potential function parameters, ιi , are detailed in Table 4 (Section 2). The material identification of the anisotropic yield surface was carried out based on the data at room temperature obtained by IFUM (Institute of Forming Technology and Machines, Leibniz Universit¨at Hannover, Germany) within the consortium of the project SPP1640 ”Joining by plastic deformation”, funded by the German Research Foundation (DFG). Then, before performing the simulations and incorporating the material hardening data, a precheck with regard to the consistency and convexity of the yield surface at different thermo-mechanical loadings was performed, see Figure 5. In this graph, it can be seen the adequate convex form of the yield surface, which was obtained during the characterization phase at room temperature (since no material data were available at different temperatures for this specific material).

50

plastic potential yield surface

40 30 20 10

0

10

20

30

40

50

80 70

60

60

50

50

40

40

30 30

20

20

10

10

-80

-60

-40

-20

0

20

40

60

-150

-100

-50

0

50

Figure 5: Characterization of the yield surface corresponding to the material PA6GF-30.

13

100

The geometric characteristics of the specimen under study were: axial length L = 60 mm, width w = 20 mm, distance between shoulders of d = 35.4 mm, and thickness t = 1 mm (Figure 6.a). The finite element discretization consisted of 615 elements subjected to mechanical and thermal loadings. The user-defined second-order fully integrated element topology (3D 20-node brick element) here developed was used for the spatial discretization. The boundary conditions imposed onto the numerical model to reproduce the experimental supporting conditions were: (i) fully constrained displacements at the clamped edge (Figure 6.a), and (ii) restrained displacements at the mobile edge (Figure 6.a) except the imposed longitudinal displacement. The geometric description of the specimen and finite element discretization are depicted in Figure 6.a. The angle (θ) was identified with loading direction angle with respect to the preferential fiber orientation.

]

a

b

[

70

60

50

40

30

20

10

0 0 0.5 1

[%]

1.5 2 3 3.5

o, RT TestData θ=0RT Test Simulation RT Sim θ=0o,RT Simulation Sim θ=0o,T=303 303K K Simulation Sim θ=0o,T=313 313K K

2.5 4

c Lengths in mm

Figure 6: Application of SFRP PA6GF-30 under thermo-mechanical uniaxial loading conditions. (a) Specimen definition. (b) Experimental–numerical correlation at room temperature and numerical results at different temperatures. (c) Longitudinal axial stress map. Note that θ identifies the loading direction angle with respect to the preferential fiber orientation.

E11 (MPa) 3970.20

E22 (MPa) 3580.30

G12 (MPa) 2221.40

ν12 0.32

ν23 0.34

Table 1: Mechanical properties for PA6GF-30.

Figure 6.b shows the stress-strain evolution curves corresponding to the simulations with regard to uniaxial tensile loading at different temperatures ranging from room temperature (TRT ≈ 20◦ C) to T = 80◦ C. The depicted data correspond to one single element in the middle of the specimen (parallel range) 14

cp (J/kg K) 1380

Z11 (m/mK) 27 · 10−6

Z22 (m/mK) 28 · 10−6

Z33 (m/mK) 28 · 10−6

K11 (W/mK) 0.23

K22 (W/mK) 0.23

K33 (W/mK) 0.23

Table 2: Thermal properties for PA6GF-30.

ζ1 4.521 × 10−4

ζ2 4.521 × 10−4

ζ3 6.624 × 10−3

ζ4 9.429 × 10−5

ζ5 3.899 × 10−3

ζ6 2.196 × 10−4

Table 3: Yield parameters ζi at room temperature for PA6GF-30.

ι1 1.073 × 10−4

ι2 1.073 × 10−4

ι3 6.782 × 10−6

ι4 2.638 × 10−5

Table 4: Plastic potential parameters ιi at room temperature for PA6GF-30.

to avoid the triaxiality effects within the vicinity of the clamps and to have an approximate uniaxial stress state. The simulation results at room temperature are compared with the experimental data obtained by IFUM for the SFRP material mentioned above, observing a very satisfactory level of agreement along the whole evolution. Subsequent simulations were also carried out for higher temperatures (Figure 6.b). In those computations it was observed that a more ductile response was predicted for higher temperatures. This is consistent with the experimental trend associated with thermal effects on the material response for a similar material as was reported in [12], so that a reasonable qualitative agreement was found. Finally, Figure 6.c shows a uniaxial longitudinal stress map of the simulation at room temperature, where a quasi-uniform stress distribution at the central area of the specimen can be identified. 5.2. Model validation: experimental-numerical correlation of monotonic thermo-mechanical tensile tests The second application here considered deals with a series of monotonic tensile loading of PA66GF35 bone-shaped specimens with different preferential material orientations (θ) with respect the loading direction, see Figure 7.a for the geometrical description. Identical finite element discretization schemes for all the computations with respect to that described in Section 5.1 were employed. These specimens were tested by De Monte et al. [12] under quasi-static conditions at room temperature (RT ) and T = 130◦ C, where the thermo-mechanical behavior of such material and the thickness effects on the experimental response were investigated. In the sequel, we focus our attention on the specimens with 1 mm in thickness, which allow the adoption of the so-called Halpin-Tsai-Nielsen [17] equations with a satisfactory level of accuracy. However, in this investigation, the material properties listed in Tables 5 and 6 for the mechanical and thermal properties, respectively, were adopted. These material properties lie within the specified ranges ruled by the corresponding ISO norms. The fiber distribution was considered to be homogeneous over the cross section of the specimen through the averaging procedure aforementioned. In addition to the previous data, the yielding and parameters, ζi , here used at room temperature and T = 130◦ C are given in Tables 7 and 8, respectively. Finally the plastic potential parameters ιi are reported in Table 9. E11 (MPa) 9919.13

E22 (MPa) 5270.13

G12 (MPa) 2415.49

ν12 0.4

ν23 0.43

Table 5: Mechanical properties for PA66GF-35.

cp (J/kg K) 1380

Z11 (m/mK) 25 · 10−6

Z22 (m/mK) 28 · 10−6

Z33 (m/mK) 28 · 10−6

K11 (W/mK) 0.67

Table 6: Thermal properties for PA66GF-35.

15

K22 (W/mK) 0.67

K33 (W/mK) 0.67

ζ1 2.648 × 10−4

ζ2 2.648 × 10−4

ζ3 3.272 × 10−3

ζ4 2.523 × 10−5

ζ5 1.338 × 10−3

ζ6 2.588 × 10−4

Table 7: Yield parameters ζi at room temperature for PA66GF-35.

ζ1 1.409 × 10−3

ζ2 1.409 × 10−3

ζ3 7.355 × 10−3

ζ4 1.224 × 10−4

ζ5 2.886 × 10−3

ζ6 1.203 × 10−4

Table 8: Yield parameters ζi at T = 130◦ C for PA66GF-35.

ι1 1.00 × 10−4

ι2 1.00× 10−4

ι3 -2.03 × 10−6

ι4 2.50 × 10−5

Table 9: Plastic potential parameters ιi for PA66GF-35.

]

a

b

[

180

160

140

120

100

80

60

40

20

0 0 0.5 1 1.5

[%]

2

TestRT θ = 0o, RT Test , RT Sim θ = 0oRT Simulation TestRT θ =90o, RT Test Sim θ = 90o, RT Simulation TestRT θ = 45o, RT Test o, RT Sim θ = 45RT Simulation

2.5 3 3.5

c

[

] 80

70

60

50

40

30

20

10

0 0 2 4

[%] 6 10

Lengths in mm

8

Test RT Test θ = 0o, 403K Simulation 403K Sim θ = 0o,RT Test Test RT θ = 90o, 403K Simulation Sim θ = 90o, 403K Test Test RT θ = 45o, 403K o, 403K Simulation Sim θ = 45RT

12

Figure 7: Application of SFRP PA66GF-35 under thermo-mechanical uniaxial loading conditions. (a) Specimen definition. (b) Experimental–numerical correlation at room temperature for different preferential fiber orientations. (c) Experimental numerical correlation at T = 403 K for different preferential fiber orientations. Note that θ identifies the loading direction angle with respect to the preferential fiber orientation.

In line with the previous application, a preliminary material parameter characterization stage prior to running the simulations was required. In this case, this procedure was carried out at different temperatures due to availability of existing material data [12]. The variation of the yield locus with temperature is depicted in Figure 8, where, again, a convex surface definition at different temperatures was obtained. 16

403

10

K

RT

-40

-30

-20

-10

5

0

10

30

20

initial yield surface at RT initial yield surface at 403 K initial yield surfaces at different temperature ranges between RT and 403K

15

3K

40

RT

5

-60

-40

-30

-20

-10

0

10

20

Figure 8: Characterization of the yield surface corresponding to the material PA66GF-35.

The experimental-numerical correlation of the monotonic tensile case are shown in Figures 7.b and 7.c for RT and T = 130◦ C, respectively, where the longitudinal stress vs. strain evolutions are depicted. In these graphs, an overall satisfactory level of accuracy between the numerical and the experimental data can be observed, evidencing the practicability and representativeness of the developed formulation at different temperatures. Note however that mid deviations with respect to the reference data were obtained for the off-axis angle tests for both temperatures here investigated. These differences were mainly attributed to the high sensitiveness of the parameter identification, which was appreciated during the preparation of the hardening data according to the procedure described in Section 2.6. In addition to this, it is worth noting that, as was reported in [12], the proposed model enabled capturing the material behavior at RT which is characterized for a more brittle and stiffer response in comparison to the corresponding performance at higher temperatures. 6. Concluding remarks In this paper, the development of a new invariant-based transversely isotropic thermo-plastic material model was presented. This model displayed the temperature dependent response of SFRPs for their use in engineering applications experiencing thermal changes. The model formulation was characterized for the incorporation of plastic deformation in fiber-direction by means of an enhanced invariant set. 17

Regarding the nonlinear material modelling, the plastic algorithm relied on a non-associative flow rule that allowed the estimation of the so-called plastic Poisson’s ratios. The current formulation preserved the thermodynamic consistency and was successfully implemented into the FE code FEAP using the user-defined routine UEL, following a monolithic thermo-mechanical numerical procedure. The applicability of the formulation developed was examined through several applications undergoing thermo-mechanical uniaxial loadings at different temperatures. The comparison of the simulations results with the experimental data for two types of SFRPs showed a satisfactory level of accuracy along the loading evolutions. Further research activities on this area would regard the application of the model to more complex structural conceptions as well as loading conditions. An additional potential area for subsequent developments would regard the consideration of dynamic and viscous effects in the current modelling strategy. In the particular case of thermo-viscoplastic constitutive formulation, a plausible rheological approach would be the Perzyna model, as proposed in [45], though different alternatives are intended to be explored. Acknowledgments The authors would like to acknowledge to Dr.-Ing. Matthias Vogler for many helpful comments and discussions. RR, SS and AD acknowledge the German Research Foundation (DFG) for the financial support through the priority program 1640 joining by plastic deformation with contract No. RO 706/6-1. JR gratefully acknowledges the Spanish Ministry of Economy and Competitiveness/FEDER (DPI2012-37187) and the Andalusian Government (Projects of Excellence No. TEP-7093 and P12-TEP-1050). Appendix A. Conversion of engineering constants for transversely isotropic materials In this appendix, the relationships between the engineering properties and the material constants introduced in Eq.(8) for transversely isotropic materials are presented. These relations read: λ = E22 (ν23 + ν31 ν13 ) /D

(A.1)

α = E22 [ν31 (1 + ν32 − ν13 ) − ν32 ] /D

(A.2)

β = E11 (1 − ν32 ν23 ) /D − E22 [1 − ν21 (ν12 + 2(1 + ν23 ))] /D − 4G12

(A.3)

µL = G12 ;

µT = G23

(A.4)

with 2 D = 1 − ν32 − 2ν13 ν31 − 2ν32 ν13 ν31

(A.5)

Appendix B. Derivatives for the local return mapping algorithm of the invariant-based coupled anisotropic thermo-plastic model In this section, we present some terms that are required for the computation of the local Newton-Raphson iteration for the algorithmic treatment of the invariant-based coupled anisotropic thermo-plastic model. We start the derivation with the linearization of the incremental yield function at the local iteration k (k) ¯pn , ϑn+1 ). The term given in Algorithm 2, where the local residual is expressed as: Rn+1 = f (k) (σ tr n+1 , A, ε (k)

∂fn+1 (k)

∂σ n+1

stemming from the numerical procedure takes the form: (k)

∂fn+1

(k)

(k) ∂σ n+1

(k)

(k)

= Kn+1 : σ n+1 + Ln+1 ,

(B.1)

where (k)

(k)

σ n+1 = Hn+1 : σ trial n+1 , 18

(B.2)

and p,(k) ¯n+1 ε

¯p,(k) ε n

=

+

(k) γn+1

r

2 (k) kMn+1 : σ n+1 k 3

(B.3)

(k)

The second term here expanded is the corresponded to

∂σ n+1

that reads

(k) ∂γn+1

(k) h i ∂σ n+1 (k) (k) = −Hn+1 : (Ce : Mn+1 ) : σ n+1 . ∂γn+1

(B.4)

The next terms within the denominator of the linearization adopt the following representation: (k)

∂fn+1 (k)

=

∂Kn+1

(k)

∂fn+1

i 1 h (k) (k) σ n+1 ⊗ σ n+1 ; 2

(k)

∂Ln+1

(k)

∂Kn+1

(k)

= σ n+1 ;

(k)

=

∂γn+1

(k)

p,(k)

p,(k)

(k)

∂Kn+1 ∂εn+1

.

(B.5)

∂¯ εn+1 ∂γn+1

(k)

The term (k)

∂Kn+1 p,(k)

∂Kn+1 p,(k) ∂¯ εn+1

given in Eq.(B.5)3 adopts the form (k)

=

∂¯ εn+1 where

(k) (k) (k)  ind  ∂ζ1(k) ∂ζ4 9 dev ind ind ∂ζ2 dev ∂ζ6 = P − P +P +2 (1 − A)⊗(1 − A) + A ⊗A A A (k) p,(k) p,(k) p,(k) p,(k) p,(k) ∂¯ εn+1 ∂¯ εn+1 ∂¯ εn+1 ∂¯ εn+1 2 ∂¯ εn+1 i=1,2,4,6 ∂ζi (B.6) # r r " p,(k) (k) (k) ∂¯ εn+1 ∂σ n+1 2 2 Mn+1 : σ n+1 : Mn+1 (k) (k) = kM : σ k + γ : . (B.7) n+1 n+1 n+1 (k) (k) (k) 3 3 ∂γ kMn+1 : σ k ∂γ

X

∂Kn+1 ∂ζi(k)

n+1

n+1

n+1

Finally (k) ∂Ln+1 (k) ∂γn+1

where

(k)

∂Ln+1 p,(k)

p,(k) ∂¯ εn+1

(k) X ∂L(k) n+1 ∂ζi

=

∂¯ εn+1

=

  (k) (k) p,(k) ∂Ln+1 ζi=3,5 ∂¯ ε

(k)

i=3,5

∂ζi

p,(k)

n+1 (k) ∂γn+1 (k)

= (1 − A)

∂¯ εn+1

.

(B.8)

(k)

∂ζ3

p,(k)

∂¯ εn+1

3 ∂ζ5 + Adev p,(k) . 2 ∂¯ ε

(B.9)

n+1

Appendix C. Element stiffness matrices and computation of the material tangent operators In this appendix, the stiffness matrices and the consistent tangent material operators corresponding to the coupled thermo-mechanical problem introduced in Eq.(58) are outlined. Appendix C.1. Element stiffness matrices Z Kuu = Z Kϑu =

Z Kϑϑ =

B

B

T

B B



∂σ ∂ε



Z B dΩ;

Kuϑ =

T



B B

  Z ˆT ˆT N ∂ [Z : Ce : εe ] N ˆ ˆ Nϑ B dΩ − ∆t ∂ε B ∆t

  Z ˆT ˆT N ∂ [Z : Ce : εe ] ˆ N ˆ ˆ Nϑ N dΩ+ ∆t ∂ϑ B ∆t

∂σ ∂ϑ



ˆ dΩ N

mech,∗ ∂Dint ∂ε

"

(C.1)

!

mech,∗ ∂Dint cp + [Z : Ce : εe ] − ∂ϑ

B dΩ

#!

(C.2)

Z ˆ N dΩ+

Bϑ


T

KBϑ dΩ

B

(C.3) mech,∗ mech with Dint = ∆tDint .

19

Appendix C.2. Consistent material tangent operators This section presents the derivation of the consistent tangent operators that are required for the consistent linearization of the coupled thermo-mechanical IBVP. From the consistency condition, one can express
∂fn+1 ∂fn+1 p ∂fn+1 ¯pn +1 = : dσ n+1 + dϑn+1 + p d¯ ε = 0. dfn+1 σ n+1 , ϑn+1 , ε ∂σ n+1 ∂ϑn+1 ∂¯ εn+1 n+1

(C.4)

Let to define A?n+1 = [Z : Ce + γn+1 [[Ce : Mϑ,n+1 ] : σ n+1 ]] ; where

An+1 = [[Ce : Mn+1 ] : σ n+1 ] ,

(C.5)

∂Mn+1 ∂ϑn+1

(C.6)

fϑ,n+1 =

∂fn+1 ∂Kn+1 ∂fn+1 ∂Ln+1 ∂fn+1 = ···· + : , ∂ϑn+1 ∂Kn+1 ∂ϑn+1 ∂Ln+1 ∂ϑn+1

(C.7)

fε¯pn+1 =

∂fn+1 ∂fn+1 ∂Kn+1 ∂fn+1 ∂Ln+1 = ···· p + : . ∂¯ εpn+1 ∂Kn+1 ∂¯ εn+1 ∂Ln+1 ∂¯ εpn+1

(C.8)

Mϑ,n+1 = We expand the following terms:

and

The following definitions of auxiliary operators are also introduced: Bn+1 = [Mn+1 : σ n+1 ] ;

M?n+1 = Mϑ,n+1 : σ n+1 : Bn+1 ;

  ˆ n+1 = M? − D? : Hn+1 : A? M n+1 n+1 n+1 ;

D?n+1 = Bn+1 : Mn+1 .

ˆ n+1 = −D? : [Hn+1 : An+1 ] ; D n+1

(C.9)

U?n+1 = D?n+1 : [Hn+1 : Ce ] . (C.10)

The derivatives in Eq.(C.4) can be recast into ∗ ∗ ∗ dγn+1 = 0, dϑn+1 + fγ,n+1 : dεn+1 + fϑ,n+1 dfn+1 = fε,n+1

where

(C.11)

r

2 γn+1 U? 3 kBn+1 k n+1 r   2 γn+1 ˆ ∗ ? p fϑ,n+1 = −fσ,n+1 : Hn+1 : An+1 + fϑ,n+1 + fε¯n+1 Mn+1 3 kBn+1 k "r # r 2 2 γn+1 ˆ ∗ fγ,n+1 = −fσ,n+1 : [Hn+1 : An+1 ] + fεpn+1 kBn+1 k + Dn+1 . 3 3 kBn+1 k ∗ fε,n+1

e

= fσ,n+1 : [Hn+1 : C ] + fε¯pn+1

(C.12) (C.13) (C.14)

Therefore dγn+1 is given by dγn+1 =

∗ ∗ dϑn+1 : dεn+1 + fϑ,n+1 fε,n+1 , m ˆ

where

∗ . m ˆ = −fγ,n+1

(C.15)

Finally, one can express:   ∗ fε,n+1 ∂σ n+1 e = = Hn+1 : C − An+1 ⊗ ∂εn+1 m ˆ  ∗   fϑ,n+1 ∂σ n+1 ∗ = . = Hn+1 : −An+1 − An+1 ∂ϑn+1 m ˆ

Ctep u,n+1 Ctep ϑ,n+1

20

(C.16)

(C.17)

From the terms stemming from the linearization of the internal dissipation, we introduce the operators: ∗   fε,n+1 + γn+1 Mn+1 : Ctep u,n+1 m  ∗  h i fϑ,n+1 = Bn+1 + γn+1 [Mϑ,n+1 : σ n+1 ] + γn+1 Mn+1 : Ctep ϑ,n+1 m

Xn+1 = Bn+1 ⊗

Yn+1

(C.18) (C.19)

Accordingly, one can express

and

mech,∗   ∂Dint = σ n+1 : Xn+1 + εpn+1 − εpn : Ctep u,n+1 , ∂εn+1

(C.20)

mech,∗   ∂Dint = σ n+1 : Yn+1 + εpn+1 − pn : Ctep ϑ,n+1 . ∂ϑn+1

(C.21)

Finally,   ∂ Z : Ce : εen+1 = Z : Ce − [Z : Ce ] : Xn+1 , ∂εn+1   ∂ Z : Ce : εen+1 = − [Z : Ce ] : Yn+1 . ∂ϑn+1

(C.22) (C.23)

References [1] Advani S.G., Tucker C.L. (1987) The use of tensors to describe and predict fiber orientation in short fiber composites. J.Rheol., 31:751-784 [2] Advani, S.G, Talreja, R. (1993) A continuum approach to determination of elastic properties of short fibre composites. Mech. Of Compos. Mater. 29(2):171–183. [3] Arif, M.F., Meraghni, F., Chemisky, Y., Despringre, N., Robert, G. (2014) In situ damage mechanisms investigation of PA66 GF30 composite Effect of relative humidity. Comp. Part B 58:487–495. [4] Arif, M.F., Saintier, N., Meraghni, F., Fitoussi, J., Chemisky, Y., Robert, G. (2014) Multiscale fatigue damage characterization in short glass fiber reinforced polyamide-66. Comp. Part B 61:55–65. [5] Behrens, B.-A., Bouguecha, A., Eckold, C.-P., Peshekhodov, I. (2012) A new clinching process especially for thin metal sheets and foils. The 15th international ESAFORM Conference on Material Forming, University of Erlangen. [6] Bellenger, V., Tcharkhtchi, A., Castaing, P. (2006) Thermal and mechanical fatigue of a PA66/glass fibers composite material. Int. J. of Fatigue 28(10):1348-1352. [7] Beyer, U. (2012) Multi-Material-F¨ ugen mittels Flach-Clinch-Technologie, Dissertation, Band 6, Technische Universit¨ at Chemnitz. [8] Boehler, J.P. (1987) Applications of Tensor Functions in Solid Mechanics. CISM Courses and Lectures, vol. 292, Springer. [9] Bonet J., Wood R.D. (1997) Nonlinear continuum mechanics for finite element analysis Cambridge University Press. [10] Botelho E.C., Figiel L., Rezende M.C., Lauke B. (2003) Mechanical behavior of carbon fiber reinforced polyamide composites. Comp. Sci. Tech. 63:18431855. [11] Celigoj, C. C. (1998) Finite deformation coupled thermomechanical problems and generalized standard materials. Int. J. Numer. Meth. Eng. 42:1025-1043. [12] De Monte, M., Moosbugger, E, Quaresimin, M. (2010) Inuence of temperature and thickness on the off-axis behaviour of short glass bre reinforced polyamide 6.6Quasi-static loading. Compos: Part A 41:859-871. [13] De Monte, M., Moosbugger, E, Quaresimin, M. (2010) Inuence of temperature and thickness on the off-axis behaviour of short glass bre reinforced polyamide 6.6Cyclic loading. Compos: Part A 41:1368-79. [14] Eidel, B., Gruttmann, F. (2003) Elastoplastic orthotropy at finite strains: multiplicative formulation and numerical implementation. Comput. Mater. Sci. 28:732-742. [15] Fu S.Y., Lauke B. (1998) The elastic modulus of misaligned short-fiber-reinforced polymers. Compos. Sci. Tech. 58(3):389– 400. [16] Gusev A.A., Hine, P.J. Ward I.M. (2000) Fiber packing and elasticity properties of a transversely random unidirectional glass/epoxy composite. Composites Sc. and Tech. 60:535–545. [17] Halpin J.C., Kardos J.L. (1976) The Halpin–Tsai equations: a review. Polym. Eng. Sci. 16:344352. [18] Hine, P.J., Lusti, H.R., Gusev, A.A. (2002) Numerical simulation of the effects of volume fraction, aspect ratio and fibre length distribution on the elastic and thermoelastic properties of short fibre composites. Comp. Sc. Tech. 62(10-11):14451453. [19] Holzapfel G. A. (2000) Nonlinear solid mechanics. John Wiley & Sons Ltd. ISBN: 978-0-471-82319-3.

21

[20] Horst, J. J., Spoormaker, J. L. (1996) Mechanisms of fatigue in short glass ber reinforced polyamide 6. Polym. Eng. Sci. 36:2718-2726. [21] Karsli, N.G. Aytac, A. (2013) Tensile and thermomechanical properties of short carbon fiber reinforced polyamide 6 composites. Compos. Part B 51:270-275. [22] Lazzarin P., Molina G., Molinari L., Quaresimin M.(1998) Numerical Simulation of SMC Component Moulding. Trans Technol Publ, Key Eng. Mater. 144:191200. [23] Imbrahimbegovic, A., Chorfi, L., Gharzeddine. F. (2001) Thermomechanical coupling at finite elastic strain: Covariant formulation and numerical implementation. Commun. Numer. Meth. Engng. 17:275-289. [24] Imbrahimbegovic, A., Chorfi, I. (2002) Covariant principal axis formulation of associated coupled thermoplasticity at finite strains and its numerical implementation. Int J. Solids Struct. 39:499–528. [25] Launay, A., Maitournam, M. H., Marmo, Y., Raoult, I., Szmytka, F. (2011) Cyclic behaviour of short glass bre reinforced polyamide: Experimental study and constitutive equations. Int. J. Plast. 27:1267–1293. [26] Lu, J., Papadopoulos, P. (2004) A covariant formulation of anisotropic finite plasticity: theoretical developments. Comput. Methods Appl. Mech. Engrg.193:5339-5358. [27] Menzel, A., Steinmann, P. (2003) On the spatial formulation of anisotropic multiplicative elasto-plasticity. Comput. Methods Appl. Mech. Engrg. 192:3431-3470. [28] Miehe, C., Apel, N., Lambrecht, M. (2002) Anisotropic additive plasticity in the logarithmic strain space: modular kinematic formulation and implementation based on incremental minimization principles for standard materials. Comput. Methods Appl. Mech. Engrg. 191:5383-5425. [29] Mlekusch, B. (1999) Thermoelastic properties of short-fibre-reinforced thermoplastics. Comp. Sc. and Tech. 59(6)911-923. [30] Reese, S. (2003) Meso-macro modelling of fibre-reinforced rubber-like composites exhibiting large elastoplastic deformation. International Journal of Solids and Structures 40 (4), 951–980. [31] Rohwer K., Rolfes R., H. Sparr (2001) Higher order theories for thermal stresses in Layered Plates. Int. J. Solids Struct.38:3673–3687. [32] Rolfes R., Noack J., Taeschner M. (1999) High performance 3Danalysis of thermo-mechanically loaded composite structures. Comp. Struct. 46:367–379. [33] Sansour C., Karsaj I., Soric J. (2007) On a formulation for anisotropic elastoplasticity at finite strains invariant with respect to the intermediate configuration. J. Mech. Phys. of Sol. 55:2406-2426. [34] Sato N., Kurauchi T., Sato S. Kamigaito, O. (1991) Microfailure behaviour of randomly dispersed short fbre reinforced thermoplastic composites obtained by direct SEM observation. J. Mater. Sci. 26:3891-3898. [35] Schr¨ oder, J., Gruttmann, F., L¨ oblein, J. (2002) A simple orthotropic finite elastoplasticity model based on generalized stressstrain measures. Comput. Mech. 30:48–64. [36] Shen H., Nutt S., Hull D. (2004) Direct observation and measurement of fiber architecture in short fiber-polymer composite foam through micro-CT imaging. Compos. Sci. Technol. 64:2113–2120. [37] Spencer, A.J.M. (1971) Theory of invariants, in: A.C. Eringen (Ed.), Continuum Physics, vol. 1, Academic Press, New York, 239-353. [38] Simo, J.C., Miehe, C. (1992) Associative coupled thermoplasticity at finite strains: Formulation numerical analysis and implementation. Comput. Methods Appl.Mech. Engrg. 98:41-104. [39] Truesdell, C., Noll, W. (1965) The nonlinear field theories of mechanics, in: S. Flgge (Ed.), Encyclopedia of Physics, vol. III/3, Springer, Berlin. [40] Ulz, M. (2009) A GreenNaghdi approach to finite anisotropic rate-independent and rate-dependent thermo-plasticity in logarithmic Lagrangean strainentropy space. Comput. Methods Appl. Mech. Engrg. 198:3262–3277. [41] Vladimirov, I. N., Pietryga, M. P., Reese, S. (2010) Anisotropic finite elastoplasticity with nonlinear kinematic and isotropic hardening and application to sheet metal forming. International Journal of Plasticity; 2010; 0749-6419; Vol. 26, 659-687 [42] Vladimirov, I.N., Reese, S. (2012) Production simulation by means of a structure tensorbased framework of anisotropic plasticity. PAMM 12 (1), 823–826. [43] Vogler, M., Andrade, F., Sch¨ opfer, J., Kolling, S. Rolfes, R. (2011) A novel transversely-isotropic 3D elastic-viscoplastic constitutive law for modeling fiber matrix composite [Conference Proceedings, Download: www.dynalook.com/]. In 8th European LS-DYNA Users Conference. Strasbourg, France, 23–24 May. [44] Vogler, M., Rolfes, R. and Camanho, P.P. (2013) Modeling the inelastic deformation and fracture of polymer composites part I: plasticity model. Mech.of Mater. 59:50–64. [45] Vogler, M. (2014) Anisotropic Material Models for Fiber Reinforced Polymers. PhD Thesis, Institute of Structural Analysis, Leiniz Universit¨ at Hannover, Hannover, Germany. [46] Wang Z., Zhou Y., Mallick P.K. (2002) Effects of temperature and strain rate on the tensile behavior of short fiber reinforced polyamide-6. Polym Compos 23(5):858-871. [47] Zienkiewicz O.C., Taylor, R.L. (1989) The Finite Element Method. Volume 1 and 2. 4th Edition, McGraw-Hill, London.

22