A computational micromechanics constitutive model for the unloading behavior of paper

International Journal of Solids and Structures 44 (2007) 7615–7632 www.elsevier.com/locate/ijsolstr

A computational micromechanics constitutive model for the unloading behavior of paper M.K. Ramasubramanian *, Yun Wang Department of Mechanical and Aerospace Engineering, North Carolina State University, 3211 Broughton Hall, Raleigh, NC 27695-7910, USA Received 1 September 2006; received in revised form 8 March 2007 Available online 10 May 2007

Abstract In this study, a computational micromechanics material model for the unloading behavior of paper and other nonwoven materials is presented. The asymptotic fiber and bond (AFB) model for paper elastic–plastic behavior [Sinha, S.K., Perkins, R.W., 1995. Micromechanics constitutive model for use in finite element analysis, In: Proceedings of the 1995, Joint ASME Applied Mechanics and Materials Summer Meeting, Los Angeles, CA, USA, Jun 28–30, 1995] has been extended to model the unloading process through a computational algorithm and implemented using the UMAT subroutine in ABAQUS finite element code. For every unloading increment, the material model assumes elastic unloading with a slope equal to the initial elastic modulus. The Jacobian matrix of the constitutive model is updated at every unloading increment by applying the incremental form of AFB model for a planar element with an elastic fiber and bond condition. A uniaxial tensile and a biaxial Mullen burst loading–unloading experiments were carried out for a paperboard sample and simulated using the model. The stress–strain curve and residual strain for the uniaxial loading were in good agreement with experimental results. The finite element model of the burst test with the AFB unloading material model predicted the general shape of the pressure versus deflection curve. However, the model over predicted the residual deflection by more than 50%. The loading portion of the pressure–deflection curve had a significant offset from experimental curves, and the nonlinearity in the unloading curve towards the end was not predicted. The discrepancies with experimental results are attributed to the burst test itself, model parameter estimation inadequacies, boundary conditions used in the FEA, and neglecting time-dependant effects. Nevertheless, the model can be useful in parametric studies relating microstructure to unloading behavior in structural problems.  2007 Elsevier Ltd. All rights reserved. Keywords: Micromechanics; Paper; Nonwovens; Constitutive model; Plasticity; Unloading; Cyclic loading

1. Introduction Paper is an orthotropic material, made up of self-bonding short cellulosic fibers preferentially oriented in the machine direction (MD), and exhibiting a distinct microstructure consisting of fibers, fiber fragments, and *

Corresponding author. Tel.: +1 919 515 5262; fax: +1 919 515 7968. E-mail address: [email protected] (M.K. Ramasubramanian).

0020-7683/$ - see front matter  2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijsolstr.2007.05.002

7616

M.K. Ramasubramanian, Y. Wang / International Journal of Solids and Structures 44 (2007) 7615–7632

voids. Due to the long slender nature of the fibers and the papermaking process of forming the web on a porous fabric from a very dilute fiber suspension, fibers generally conform to the plane of the paper, forming essentially a two-dimensional layered structure. On a larger size scale (millimeters), paper structure exhibits variation in mass density due to the non-uniform nature of dispersion during the papermaking process. These non-uniformities affect the strength and appearance of paper. Paper exhibits a strong nonlinear elastic–plastic behavior under monotonic uniaxial loading (Niskanen and Karenlampi, 1998). Upon unloading from plastic regime, the paper exhibits permanent deformation and time-dependent recovery, and upon reloading, exhibits a hysteresis (Bronkhorst, 2003). 1.1. Elastic behavior The elastic behavior and strength of paper have been studied extensively in the past through a combination of experimental, semi-empirical, and analytical approaches (Niskanen and Karenlampi, 1998). Constitutive models to predict the elastic behavior of paper in terms of its constituents have been developed by several researchers working at different size scales starting from hydrogen bonds to a macroscopic sheet continuum (Cox, 1952; Corte and Kallmes, 1962; Van den Akker, 1970; Page et al., 1979; Perkins, 1986, 1990; Batten and Nissan, 1987; Schulgasser and Page, 1988; Schulgasser and Grunseit, 1990; Suhling, 1990; Lu et al., 1996; Ostoja-Starzewski and Stahl, 2001; Xie et al., 2002). While all the referenced work above are 2-D planar or 2-D laminated composite or mosaic models, a 3-D network model was presented by Heyden and Gustafsson (2001). This model for paper takes a cubic unit cell, with dimensions 1.2–4.8 times the fiber length, and places fibers in them randomly in 3-D. The fiber crossings were considered bonded when the proximity of fibers were within a prescribed tolerance. Fibers and bonds were assumed to be linear elastic. The bonds were assumed to exhibit a stick–slip behavior. The fibers were modeled as linear elastic Bernoulli beams of constant curvature in a finite element algorithm to simulate uniaxial loading in all three directions. The results from such a simulation were the elastic parameters of the network, and strength related parameters such as fracture energy. The model was used to simulate fiber networks containing 173 fibers for out-of-plane density to inplane density ratios of 1.0–0.2. A ratio of 1 will represent fiber fluff and smaller ratios will tend towards two-dimensional densified paper structure. Only elastic moduli alone were calculated for all density ratios. The stress–strain curve was simulated only for ratios of 1.0 and 0.8. Below 0.8, the degrees of freedom increased rapidly, presumably making it computationally expensive. No experimental results were presented to validate the model. Thus, the model is useful for understanding fiber fluff and very low density paper behavior and not applicable for a majority of paper and paperboard microstructures. Using this approach, it would be almost impossible to model realistic macroscopic structural problems involving paper materials. Similar approaches to defining a 3-D unit cell accurately with statistics and geometry have been presented in the area of textiles and nonwovens (Komori and Makishima, 1977; Komori and Itoh, 1991). In all cases, three-dimensional network models very quickly become computationally unmanageable. Hence, these models have only been applied in elastic uniaxial loading cases. 1.2. Inelastic behavior While extensive work has been carried out on the elastic behavior of paper, inelastic behavior of paper has received relatively less attention. Treating paper as a homogeneous continuum, Johnson and Urbanik (1984) formulated a nonlinear elastic plate theory and defined a strain energy density function W. The form of the strain energy density function was not known. The stress–strain relationship for the material was obtained in a hyperbolic tangent functional form through curve fitting uniaxial experimental data. The hyperbolic tangent function was then integrated to obtain the strain energy density function W. Once the strain energy density function is obtained, stress–strain relations in other directions were obtained through the relationship oW rij ¼ oe . Suhling (1990) used a hyperelastic model in a similar framework to predict the nonlinear stress–strain ij curve. In this approach, although the nonlinear stress–strain curve is correctly predicted, the curve is elastic. It cannot predict the unloading path correctly. Castro and Ostoja-Starzewski (2003) extended the hyperbolic tangent model to fit both uniaxial and biaxial test data. A three-dimensional, anisotropic constitutive model has been presented to model the in-plane elastic–plastic deformation of paper and paperboard (Xia et al.,

M.K. Ramasubramanian, Y. Wang / International Journal of Solids and Structures 44 (2007) 7615–7632

7617

2002), which uses a continuum plasticity model with an evolving hardening parameter to predict the anisotropic elastic–plastic behavior of paper. Makela and Ostlund (2003) presented a continuum orthotropic elastic–plastic material model for paper applying the isotropic plastic equivalent (IPE) (Karafillis and Boyce, 1993) concept. By equating the IPE-material and the orthotropic material energy dissipation rates, the incremental plastic strain tensor for the IPE-material and that of the anisotropic material were also related by the IPE transformation tensor. Strain hardening behavior was modeled with a modified form of Ramberg– Osgood relation and the expressions for plastic strain increments were derived. The plane-stress constitutive model was implemented using the user-defined material (UMAT) subroutine in ABAQUS/Standard finite element code. Uniaxial tests were simulated and compared with experiments. Unloading behavior was not considered. A continuum model that combines nonlinear elasticity and bounding surface plasticity, and implemented using the UMAT subroutine in ABAQUS model, for modeling through-thickness elastic–plastic behavior of paper was presented by Stenberg (2003). In this model, elastic plastic unloading and reloading was considered in compression, out-of-plane tension was treated as elastic. Inelastic behavior of paper has been modeled through a micromechanical analysis by Ramasubramanian (1987), Ramasubramanian and Perkins (1988), and experimentally studied by Furukawa et al. (1991). A fiber (inclusion phase) with portions of homogenized crossing fibers (matrix phase) for shear lag load transfer was defined as the representative volume element (RVE). Network geometry, fiber length and orientation distributions, sheet density, and fiber and fiber-to-fiber bond constitutive behavior were taken into account. Fiber behavior can be different in tension and compression thereby the possibility of fiber buckling under compression is included. Monotonic increase in load was simulated. Stress–strain curves and lateral contraction ratios were predicted. Unloading was not considered. The model is computationally intensive in that for every sheet strain value, the fiber axial strain distribution and the fiber-to-fiber bond shear stress distribution in every direction has to be computed and the portion of the fibers and bonds that are elastic and plastic have to be computed. To reduce the number of computations without significant compromise in quality, Sinha (1994) proposed an asymptotic fiber and bond model to handle the analysis of RVE used by Ramasubramanian and Perkins (1988). Instead of determining the lengths of portions of every fiber and bond that is elastic or plastic at every orientation and strain level, the asymptotic model assumes that when the maximum strain in a fiber or bond reach its critical stress value, the corresponding fiber or bond yields completely along its entire length. Thus, the end effects and shear lag mechanism are not addressed. The simplified model reduced considerably the amount of numerical calculation required without compromising the quality of the results (Sinha, 1994). An incremental form of the asymptotic model was implemented in as the material model using the UMAT procedure in ABAQUS finite element code and monotonic uniaxial tensile loading was studied (Sinha and Perkins, 1995). A two-dimensional stochastic computational network model to simulate uniaxial tension tests was presented by Bronkhorst (2003). A two-dimensional stochastically generated network consisting of fibers in a 10 mm · 10 mm area was studied. The fiber-to-fiber bond was assumed rigid and the fibers behavior was modeled with continuum plasticity models. A finite element model of the structure was created in ABAQUS finite element code, modeling fibers as isotropic elastic–plastic beams, and loaded in uniaxial tension. The stress– strain curve generated was compared with experiments. One cyclic tensile test result was presented and compared with experiment. The discrepancy in later stages of unloading was attributed to time-dependent effects not accounted for in the model. While this model could be used to study larger problems, the RVE of the structure is a macroscopic sheet (10 mm · 10 mm, for example). Thus, local effects of fiber-to-fiber bond yielding and fiber failure could not be taken into account. Reducing the size of RVE may be one option, but it might increase the computational requirements as reported in other stochastic network model approaches (Heyden and Gustafsson, 2001). 1.3. Unloading and cyclic loading behavior While there has been good progress in recent years in modeling the inelastic behavior of paper, very little has been done to model the unloading and cyclic loading behaviors (Ramasubramanian and Wang, 1999). As discussed in previous section, Bronkhorst (2003) showed a result from cyclic loading simulation based on

7618

M.K. Ramasubramanian, Y. Wang / International Journal of Solids and Structures 44 (2007) 7615–7632

stochastic network model with rigid bonds and elastic plastic fibers. Kang et al. (2006) presented a numerical approach for uniaxial cyclic deformation of discontinuously reinforced (particulate and short fiber) metalmatrix composites. The fibers were treated as an elastic material and the matrix was treated as an elasto-plastic material. Infinitesimal elasto-plasticity with isotropic elasticity, along with kinematic hardening and an evolution rule for back stress was used. The constitutive model was revised in an incremental form and placed in the UMAT subroutine in ABAQUS. Cyclic deformation of metal-matrix composites was simulated through a single fiber and matrix model and verified with experiments. The effects of weak interfacial bonding, fiber orientation, thermal residual stress and other microstructural features, and stress relaxation were not taken into account. Nevertheless, the model offers a good framework to develop a computational constitutional model for unloading and cyclic loading behavior of short fiber composites such as paper. In another study, (Doghri and Tinel, 2006), an incremental micromechanics constitutive model was proposed for rate independent inelastic problems in short fiber reinforced composites. The incremental form of model was implemented through the UMAT procedure in ABAQUS and unloading cases were simulated, and compared with experiments. 2. Objectives In this paper, we present a Computational methodology for simulating the unloading behavior of paper in a micromechanics framework already established for the simulation of inelastic behavior of paper (Ramasubramanian and Perkins, 1988; Sinha, 1994; Sinha and Perkins, 1995). We present the incremental form of the asymptotic fiber and bond model and the unloading algorithm implemented into ABAQUS finite element code using the UMAT subroutine. We report results from simulation of a uniaxial and a biaxial loading–unloading problem, and compare results with tensile tests and biaxial tests on a Mullen burst tester, respectively. 3. Micromechanics model 3.1. Full plasticity model The representative volume element (RVE) or the mesoelement for the model consists of a fiber and a portion of all crossing fibers attached to it as shown in Fig. 1 (Ramasubramanian and Perkins, 1988). The fiberto-fiber bonds are distributed on the fiber and essentially transfers load to the reference fiber in the RVE through a shear lag mechanism. The effect of fiber-to-fiber bonds is smeared or homogenized along the fiber like a surrounding deformable matrix. A bilinear elastic–plastic material model is used to model the constitutive behavior of the fiber and then bond, shown in Figs. 2 and 3 (Ramasubramanian, 1987; Sinha, 1994). In this two-slope model, the elastic modulus of fiber is the same in tension and compression. In order to address localized failure, such as fiber buckling, the compressive strength and tangent modulus (second slope) of the fiber cell wall material could be much smaller than their respective tensile counterparts. The constitutive model of fiber in tension would then be given by

θ Fiber

Machine Direction

RVE Element Boundary Bonds

Fig. 1. Representative Volume Element (RVE) for the micromechanics model.

M.K. Ramasubramanian, Y. Wang / International Journal of Solids and Structures 44 (2007) 7615–7632

7619

σf σ pt

E2t Ef

ε pc

ε pt

E2c

εf

σ pc

Fig. 2. Fiber constitutive model.

τ

τb

G2

Gb γb

γ

Fig. 3. Fiber-to-fiber bond constitutive model.

rf ¼ Ef ef ;

ð1Þ

rf < rpt

rf ¼ rpt þ E2t ðef  ept Þ;

rf P rpt

ð2Þ

and in compression rf ¼ Ef ef ;

jrf j < jrpc j

rf ¼ rpc þ E2c ðef  epc Þ;

ð3Þ jrf j P jrpt j

ð4Þ

where rf = fiber stress, ef = fiber strain, Ef = fiber elastic modulus, E2t = fiber tangent modulus (second slope) in tension, E2c = fiber tangent modulus (second slope) in compression, rpt = fiber yield stress in tension, rpc = fiber yield stress in compression, ept = fiber yield strain in tension, and epc fiber yield stress in compression. The fiber-to-fiber bond, the constitutive behavior can be described by sb ¼ G b cb ;

jsb j < jsp j

sb ¼ sp þ G2 ðc  cb Þ;

ð5Þ jsb j P jsp j

ð6Þ

where, sb = bond shear stress, cb = bond shear strain, Gb = bond shear elastic modulus, G2 = bond tangent modulus (second slope), sp = bond yield stress, and cp = bond yield strain. One half of the RVE is shown in Fig. 4. The length Lf represents the portion of the fiber that is elastic, Lb represents portion of the bond that is plastic, and L is the half length of RVE. The location of these boundaries were computed in a full plasticity model (Ramasubramanian, 1987) for the following three possibilities, namely, elastic fiber and bonds, elastic fiber–plastic bond, and plastic fiber and bond for an RVE. Once, these boundaries were determined, the total work function was numerically

7620

M.K. Ramasubramanian, Y. Wang / International Journal of Solids and Structures 44 (2007) 7615–7632

Lf Lb L Fig. 4. Portion of a RVE showing location of elastic–plastic interfaces from the center of the element.

computed. Determining the location of elastic–plastic boundaries for every macroscopic applied strain and fiber orientation made this computationally intensive and limited the scope of this approach. 3.2. Asymptotic fiber and bond model In the asymptotic model, the average fiber stress was computed once a critical stress was reached in the middle of the fiber, thus homogenizing the fiber stress in the RVE. The average stress is computed by Z 1 L f ¼ r rf dn ð7Þ L 0 Four critical strains values are used to determine the fiber and bond conditions, namely, ecb, ecf, ebf, and efb. Starting with elastic fiber and elastic bond, ecb is associated with the initiation of bond yielding before fiber yielding, and ecf is associated with the initiation of fiber yielding before bond yielding. The strain ebf is associated with the initiation of fiber yielding after the bond has already yielded, and efb is associated with the initiation of bond yielding after the fiber has already yielded. If ebf is greater than efb, then bond would yield first and vice versa. In this paper, we only consider the case of bond yielding while the fiber remains elastic. From equilibrium considerations of an axial section of the RVE, the fiber stress expressions can be derived for each case using the full plasticity model (Ramasubramanian and Perkins, 1988). The expressions can be then simplified using the asymptotic approximation for this work (Sinha, 1994) and the results are shown in Eqs. (8) and (9). Case of bond yielding before fiber (ecf > ecb) Average fiber stress for elastic fiber and elastic bond is given by   tanh aL  f ¼ E f es 1  r ð8Þ aL Average fiber stress for elastic fiber and plastic bond is given by     tanh aL tanh ab L f ¼ Ef ecb 1  r þ Ef ðes  ecb Þ 1  aL ab L

ð9Þ

where es is the sheet strain in the direction of the RVE. In terms of the applied macroscopic strains, and the orientation of the RVE, the sheet strain is given by es ¼ ex cos2 h þ ey sin2 h þ 2exy cos h sin h

ð10Þ

and a2 ¼

Gb we E f Af t b

ð11Þ

a2b ¼

G2 we E f Af t b

ð12Þ

M.K. Ramasubramanian, Y. Wang / International Journal of Solids and Structures 44 (2007) 7615–7632

7621

where we is the width of the RVE, Af is the cross section area of the fiber, and tb is the thickness of the bond region in the RVE. The critical strain expression for the case of elastic fiber and plastic bond in terms of bond yield stress sp is (Sinha, 1994):   sp we L 1 ecb ¼ ð13Þ Ef Af aL tanh aL The sheet stress rx and ry are then obtained from the integral sum of the appropriate average fiber stress component for fibers in all orientations, Z p=2 q f ðhÞ cos2 hf ðhÞdh rx ¼ s ð14Þ r qf p=2 Z p=2 q f ðhÞ sin2 hf ðhÞdh ry ¼ s ð15Þ r qf p=2 where f(h) is the von Mises function used to describe the fiber orientation distribution, qs is the sheet density and qf is the density of the fiber cell wall material (Ramasubramanian and Perkins, 1988). 3.3. Incremental form of the asymptotic fiber and bond model In the incremental form (Sinha and Perkins, 1995), the incremental sheet strain in the direction of the RVE in terms of applied incremental two-dimensional macroscopic strains field is given by Des ¼ Dex cos2 h þ Dey sin2 h þ 2Dexy cos h sin h The total axial strain in the direction of RVE is given by the sum of the incremental strains n X ðDes Þi es ¼

ð16Þ

ð17Þ

i¼1

The elemental Drij is assumed to be sufficiently small such that it is still linearly related to Deij by the instantaneous slope of the incremental stress–strain curve. The fiber stresses in the incremental form can be derived in a fashion similar to the full asymptotic approach. For the case of elastic fiber and elastic bond (cf. Eq. (8)), the average incremental stress is given by (Sinha, 1994)   tanh aL hD rf i ¼ Ef ðDes Þ 1  ð18Þ aL The incremental sheet stresses are obtained by Z qs p=2 Drx ¼ hD rf ðhÞi cos2 hf ðhÞdh qf p=2 Z qs p=2 Dry ¼ hD rf ðhÞi sin2 hf ðhÞdh qf p=2 The total sheet stress is a sum of all the previous sheet stress increments X ðDrx Þ rx ¼ X ðDry Þ ry ¼

ð19Þ ð20Þ

ð21Þ ð22Þ

4. Computational material model The incremental asymptotic fiber and bond model was incorporated into ABAQUS/Standard for the solution of problems with complex geometries and multi-axial loading, through the UMAT subroutine (ABAQUS, 1998). This subroutine will be called at all material calculation points of elements and the subroutine

7622

M.K. Ramasubramanian, Y. Wang / International Journal of Solids and Structures 44 (2007) 7615–7632

must update the stresses and solution-dependent state variables to their values at the end of the increment. The user-defined material   model is incorporated by providing complete expressions for the terms in the Jacobian matrix oDrij =oDeij . The Jacobian matrix represents the incremental change in material properties for an element and is determined by applying the incremental micromechanics model for a planar element in our case. The constitutive relationships based on the asymptotic fiber and bond model in its incremental form can be written in a matrix form as shown in Eq. (23). 9 2 9 8 38 k 11 k 12 k 13 > > = = < Drx > < Dex > 6 7 Dry ¼ 4 k 21 k 22 k 23 5 Dey ð23Þ > > > > ; ; : : Drxy Dexy k 31 k 32 k 33 where k 11 k 12

q ¼ s qf q ¼ s qf

k 13 ¼ k 22 ¼ k 23 ¼

qs qf qs qf qs qf

Z

1

1

1

1

0

Z

F sin2 h cos2 hf ðhÞf ðkÞdh dk

ð25Þ

F sin h cos3 hf ðhÞf ðkÞdh dk

ð26Þ

F sin4 hf ðhÞf ðkÞdh dk

ð27Þ

F sin3 h cos hf ðhÞf ðkÞdh dk

ð28Þ

p=2

Z

p=2

Z

p=2

p=2

0

Z

ð24Þ

p=2

0

Z

F cos4 hf ðhÞf ðkÞdh dk

p=2

0

Z

p=2

p=2

0

Z

Z

1

Z

p=2

p=2

k 33 ¼ k 12

ð29Þ

The F functions are derived from micromechanics which are essentially the fiber stresses in the RVE for different loading cases (Sinha, 1994). In our study, we have considered the fiber to be elastic and bond to be elastic–plastic. The two relevant cases are given below. Fiber and bond elastic (es 6 ecb)   tanh aL F ¼ Ef 1  ð30Þ aL Fiber elastic and bond plastic (es > ecb)   tanh ab L F ¼ Ef 1  ab L

ð31Þ

4.1. General unloading algorithm When a fiber is unloading, es would have an opposite sense to that of Des. Both fiber and bond would behave elastically. The F function that corresponds to elastic fiber and bond condition (Eq. (30)) is used in the computation of the Jacobian matrix and the corresponding change in rij. The algorithm for unloading calculations is shown in Fig. 5. For every time step, each iteration of the Finite Element equilibrium conditions, and at each integration (Gauss) point of the mesh, ABAQUS calls UMAT as shown in Fig. 5. Consider now an undeformed sheet subjected to tensile loading in the machine direction. The values of the model parameters are defined in the input file for use in the UMAT subroutine. In this work, ecb is greater than ecf in tension and compression. Therefore, bond would yield first when the critical strain ecb is reached. The fiber orientation distribution, f(h), is discretized into 30 increments varying from p/2 to p/2. For a given loading condition on the sheet, the algorithm begins with h = p/2 and calculates the corresponding es (h). The critical strains are calculated. Suppose es(h) is loaded in tension (sgn(es) = sgn(Des) and Des > 0), beyond the yield point ecb.

M.K. Ramasubramanian, Y. Wang / International Journal of Solids and Structures 44 (2007) 7615–7632

Call UMAT Subroutine

Calculate εs(θ) and Δεs(θ ) Orientation dependent model parameters Ef, E2, G b, G2, a, ab, Critical strain ε cb

Read Input Parameters ε x , ε y , ε xy , σ x , σ y , and σ xy for the kth time step; Δε x , Δε y , and Δε xy for the (k+1)th time step

Sign of σ f ( θ) = Sign of Δε f(θ) ?

Define model parameters E f and E 2 . Gb and G2 τ p ,ε pt ,ε pc , ρ s , ρ f , κ

I=1

7623

YES Loading process

NO

NO

Unloading Process: Calculate F function using the expression for elastic fiber and bond

I= I+1= 31?

θ ( I ) = -π / 2 + ( I -1) (π / 30) YES Determine kij of Jacobian Matrix using Simpson’s 1/3 rule.

Calculate Δσ ij

σ ijk +1 = σ ijk + Δσ ij Return to ABAQUS Fig. 5. Incremental form of asymptotic model UMAT subroutine implementation flow chart.

Bond would yield and the total strain of es could be broken into two components, i.e. the elastic and plastic strain component. The es,el is recovered when the sheet is unloaded. The kij in the Jacobian matrix are obtained by evaluating the integral of the F functions numerically using the Simpson’s 1/3 rule at every material calculation point. The incremental stress is then calculated and the total stress value is returned to ABAQUS. The algorithm is capable of addressing cyclic loading. During cyclic loading, the critical strains have to be modified accordingly so that the ensuing plastic deformation would occur at the correct stress limit. Moisture and temperature effects can be easily added to the incremental constitutive relationships (Sinha, 1994). 4.2. Cyclic loading Consider a fiber with an orientation h, to the MD and loaded in tension until the bond in the RVE has yielded. When it is unloaded, the fiber strain would have a nonzero plastic component, es,pl. This nonzero value of the fiber plastic strain component is added to the critical strains in compression. After complete unloading, if the fiber is subjected to compression, the updated critical strains allow the RVE to yield at the correct stress limit in compression.

7624

M.K. Ramasubramanian, Y. Wang / International Journal of Solids and Structures 44 (2007) 7615–7632

The elastic and plastic strain components of es are also stored for decision making in the algorithm and updating the critical strains for subsequent loading cycles. In this study, only the bond plastic critical strain ecb for the RVE is considered. It is updated as follows. ecb ¼ ecb;initial þ es;pl

ð32Þ

For a fiber subjected to multiple cycles of loading and unloading process, the critical strains are updated by the same equation, except the ecb,initial. is replaced by the previous value of ecb. When subjected to compression after unloading, the original critical strains are adjusted by adding the cumulative absolute value of the plastic strains, eoffset (cf. Eq. (33)). In tension, the offset value is subtracted from the original critical strains.  n  X i  ð33Þ eoffset ¼ es;pl  i¼1

At any instance of unloading, the F function for elastic fiber and elastic bond is used to calculate the material contribution to the Jacobian matrix. 5. Results 5.1. Uniaxial loading–unloading The nonlinearity of the sheet stress–strain behavior is due to the nonlinear behavior of the fibers and fiberto-fiber bonds, and could be modeled by either allowing the fiber or the bond to yield first. Since, in this study the later case is used to demonstrate the method, the model parameters are reduced to j, Ef, Gb, sp, and G2. The model parameters were determined by a curve fitting procedure using uniaxial stress–strain data for the material, in this case, a bleached softwood kraft board paper cup stock material (density = 720 kg/m3, thickness = 0.48 mm). Tests were conducted on an InstronTM tensile tester with rectangular samples, 25.4 mm · 177.8 mm, at a displacement rate of 25.4 mm/min. In order to determine the orientation parameter j, and the lateral contraction ratio m, the samples were first relieved from drying induced stresses by exposing to high humidity (85% RH, 23 C) for 24 h and later conditioning the samples in a TAPPI standard testing conditions (50% RH, 23 ± 1.0 C) for 24 h. The ratio of EMD/ECD of these samples measured through uniaxial tests is then a function of j and the Poisson’s ratio m, and the relationship is given by (Wang, 2000)  R p=2  2 cos h  mMD sin2 h cos2 hf ðhÞdh EMD p=2 ¼ R p=2  2  ECD sin h  mCD cos2 h sin2 hf ðhÞdh

ð34Þ

p=2

The orientation distribution f(h) is a function of the orientation parameter j Perkins, 1986). A family of EMD/ ECD versus j (0.0 6 j 6 1.0) curves were generated for a range of discrete values for m 0.06m61.0) using Eq. (34). For the experimentally measured value of EMD/ECD (2.06), the range of possible values for j were determined from the plots and the average value was chosen, in this case, 0.6. Once the value of j is determined, the value of lateral contraction ratio is determined by interpolation (Wang, 2000). In the elastic region, the elastic modulus is a function of fiber and bond elastic properties. In a sensitivity study by Sinha (1994), Ef was found to have a strong influence the sheet elastic modulus, while Gb had negligible effect. Meanwhile, the sheet yield stress was dependent on sp and the slope of the stress–strain curve in the plastic region was influenced by G2. The values for these parameters were determined using the procedure discussed by Wang (2000) through curve fitting uniaxial tensile stress–strain curves, once the elastic parameters were determined. The hardening parameters (due to drying restraint effects HEf0, HGb0,Hsp0, HG20) could be obtained in similar manner using the experimental stress–strain data from samples that without relieving drying stresses. Table 1 shows the list of parameters and their values chosen through this procedure for the simulation. Experiments were carried out by first loading the sample into the plastic regime and unloading subsequently at a rate of 25.4 mm/min. The displacements at zero force are difficult to obtain in a tensile tester because the samples tend to buckle when nearing zero force. Therefore, the displacement at zero force was approximated by using the displacement when the applied load is close to zero.

M.K. Ramasubramanian, Y. Wang / International Journal of Solids and Structures 44 (2007) 7615–7632

7625

Table 1 Experimentally determined and estimated material parameters for UMAT input Parameter

Values

j l axM. ayD exD, evD Ef0 HEf0 Gbo HGbo sp0 Hsp0 Gb2 HG20

0.6 0 (MD) 0.045, 0 0,0 18 GPa 0 15 MPa 0 6.5 MPa 2.0 10 kPa 70

5.2. Asymptotic fiber and bond model calculations The uniaxial loading–unloading problem can be simulated directly using the asymptotic fiber and bond model. Since we assume that unloading is elastic, with a slope equal to the initial elastic modulus, the computation of residual stress is quite straight forward and can be computed as follows. For a uniaxially strained sheet in the x-direction, Eqs. (8), (10) and (14) can be combined to determine the sheet stress as: rx ¼

Z

qs E f ex qf

p=2



cos2 h  m sin2 h

p=2



  tanhðaLÞ 1 cos2 hf ðhÞdh aL

ð35Þ

where m is Poisson’s ratio and the sheet elastic modulus is given by q Ex ¼ s Ef qf

Z

p=2



2

2

cos h  m sin h

p=2





 tanhðaLÞ 1 cos2 hf ðhÞdh aL

ð36Þ

where Ex is the elastic modulus of the sheet. The asymptotic fiber and bond model described in Section 3.2 is used to generate the uniaxial tensile loading portion of the stress–strain curve. Upon unloading the sheet from the nonlinear portion of the stress–strain curve, the sheet unloads elastically (elastic fibers and bonds). Then, the elastic strain recovered by the sheet at the end of unloading would be given by the sheet stress at the point of unloading divided by the sheet elastic modulus, given by R p=2

ex;el ¼

f ðhÞ cos2 hf ðhÞdh r

 R p=2  Ef p=2 cos2 h  m sin2 h 1  tanhðaLÞ cos2 hf ðhÞdh aL p=2

ð37Þ

The integrals in Eqs. (35)–(37) are evaluated numerically using multiple applications of Simpson’s 1/3 rule. Once the elastic recovered strain is obtained, the non-recoverable plastic strain, epl, is given by epl ¼ etotal  eel

ð38Þ

MD and CD tensile elastic–plastic unloading experiments described above for paper cup stock material were simulated using the asymptotic fiber and bond model. The results are shown with the experimental curves in Figs. 6 and 7. The deviation between the predicted plastic residual strain and the experimental values varies

7626

M.K. Ramasubramanian, Y. Wang / International Journal of Solids and Structures 44 (2007) 7615–7632 3.5E+07 model experimental 3.0E+07

MD Stress, Pa

2.5E+07

2.0E+07

1.5E+07

1.0E+07

5.0E+06 0.005

0.01

0.015

MD Strain Fig. 6. MD unloading curve and residual strain from model and experiment. The unloading was done at a target strain value of 0.01726 m/m.

from 5% to 13% for the MD curves and between 1% and 2% for the CD curves. The measured plastic strains at zero stress are all higher than the model predictions as the experiments did not take the sample down to zero stress. The asymptotic model predicts uniaxial unloading effects satisfactorily. 5.3. Finite element simulation of uniaxial unloading The unloading algorithm implemented within UMAT in ABAQUSTM was used to simulate the uniaxial tensile test. The sample was modeled using a 4-node bilinear plane stress element (CPS4). Each node of this element has two degrees of freedom, i.e. ux and uy. The nodes on the left boundary were constrained in the two degrees of freedom while the right boundary was displaced (cf. Fig. 8) by the amount of the average MD displacement obtained in experiments, i.e. 3.120 · 103 m. For unloading, the right boundary was displaced in the opposite direction until the average experimental displacement at unloading was reached, i.e. 1.752 · 103 m. The results from the center element on the right boundary were used to obtain the stress–strain curves. A plot of the 100-node model solution and comparison with the closed form asymptotic fiber bond model solution are shown in Fig. 9. Results show excellent agreement with the asymptotic fiber and bond calculation results, validating the finite element implementation of the incremental asymptotic fiber bond model. 5.4. Biaxial loading–unloading In order to subject the paper sample to a state of biaxial loading, we used the standard Mullen burst tester. In this test, a paper sample is gripped between two circular disks (like two flat washers) and a hydraulic piston moves a fluid against a neoprene diaphragm contacting the circular unsupported sample area in the disk thereby increasing the pressure until the sample fails. The maximum pressure is recorded as the burst strength. The standard Mullen tester, shown in Fig. 10, was instrumented to measure the sample central deflection as a function of applied pressure. A schematic of our test setup is shown in Fig. 11. The measurement of central deflection of the sample presents a difficulty due to limited space available between the sample and the clamp. A 12.7 mm · 12.7 mm metal foil target was placed in the center on the sample with a double sided adhesive

M.K. Ramasubramanian, Y. Wang / International Journal of Solids and Structures 44 (2007) 7615–7632

7627

2E+07 1.8E+07

model experimental

1.6E+07

CD Stress, Pa

1.4E+07 1.2E+07 1E+07 8E+06 6E+06 4E+06 2E+06 0.01

0.02

0.03

CD Strain Fig. 7. CD unloading curve and residual strain from model and experiment. The unloading was done at a target strain value of 2.5 · 102 m/m.

Fig. 8. The tensile test sample was modeled with 100 4-node bilinear plane stress elements in ABAQUSTM.

tape. A capacitance sensor (HPB-500 Button Probe, by Capacitec), in the form of a circular disk, measuring 12.7 mm in diameter, was placed above the target. The sensor measures change in capacitance caused by a change in distance between the sensor and the metal target. The output signal is converted to voltage by the amplifier and measured by a data acquisition system (CIO-DAS1400 DAQ board from Computer Boards, Inc., LabTech NotebookTM data acquisition software, and a PC). The hydraulic pressure is measured with using an in-line pressure transducer in the hydraulic circuit and interfaced to the same data acquisition system. The output voltage of the pressure transducer varies linearly with the applied pressure. In the experimental work, the sample was clamped down in the burst tester. The 3.0-cm diameter sample area was subjected to the applied pressure and the displacement data was collected. 5.5. Finite element simulation of mullen burst loading–unloading test The test area was modeled with 96 shear deformable shell elements in ABAQUS. The innermost ring was made up of 6-node triangular thin shell elements (STRI65), while the rest were 8-node doubly curved thin shell elements (S8R5). The boundary of the model was constrained in the displacement degrees of freedom only, i.e. ux,uy, and uz. Although the sample was clamped down during the test, the rotational degrees of freedom were

7628

M.K. Ramasubramanian, Y. Wang / International Journal of Solids and Structures 44 (2007) 7615–7632

Fig. 9. Asymptotic fiber bond model and ABAQUS finite element implementation comparison, MD orientation.

Fig. 10. Standard mullen burst tester and view of the clamp area.

not restrained in the finite element analysis to simulate correctly the typical central failure location in a burst test. The shell elements used in ABAQUS, STRI65 and S8R5, require shear elastic moduli values to calculate the transverse shear stiffness. Since the asymptotic fiber and bond model does not predict out-of-plane shear properties, experimental values obtained for the samples using ultrasonic measurements (Baum et al., 1981). For the coated paper cup stock used in this work, we determined G23 = 0.109 GPa and G13 = 0.119 GPa. The corresponding transverse stiffness are given by 5 k ts11 ¼ G13 t; 6

5 k ts22 ¼ G23 t; 6

k ts12 ¼ 0:0

ð39Þ

where t is the average sample thickness (0.48 mm). Loading unloading experiments and simulations were carried out at four different maximum pressure values, namely, 0.49, 0.51, 0.54, and 0.72 MPa. Results are shown in Fig. 12.

M.K. Ramasubramanian, Y. Wang / International Journal of Solids and Structures 44 (2007) 7615–7632

7629

Fig. 11. Test setup for measuring pressure–deflection curve with a Mullen tester.

Fig. 12. ABAQUS FEA model with incremental form of asymptotic fiber bond material model compared with mullen burst test.

7630

M.K. Ramasubramanian, Y. Wang / International Journal of Solids and Structures 44 (2007) 7615–7632

From Fig. 12, it can bee seen that the agreement between FEA simulation and the experimental results are not very good, although a general shape of a nonlinearly increasing pressure versus deflection behavior and a residual deflection upon unloading are evident. However, several important observations can be made to explain the discrepancies. Historically, Mullen burst test is not the most reproducible test. It was selected because it offers a quick and convenient way to conduct a reasonable approximation to a standard biaxial test, and it simulates some end use applications for paperboard. The hydraulic loading as indicated by the pressure of the fluid is only an indirect way to measure the force applied to the sample. The differences in measured pressure of 0.49, 0.51, and 0.55 MPa are not significant in a burst tester, and should be considered essentially the same test. While the first two experimental curves (0.49 and 0.51 MPa) are sufficiently close in terms of peak central deflection (about 2.3 mm), and residual deflection after unloading (about 1 mm), the third test in the same group (Peak pressure = 0.55 MPa), exhibited an anomalous behavior; the peak deflection was 2.6 mm and the residual deflection was 1.8 mm, exceeding that of the higher pressure test at 0.72 MPa. The highest pressure test at 0.72 MPa measured a peak deflection of 3 mm and residual deflection of 1.75 mm. The model results, on the other hand, were consistent with one another; the residual deflections were 1.65, 1.7, 1.8, and 2.4 mm, for pressures 0.49, 0.51, 0.54, and 0.72 MPa peak pressures, respectively. Thus, the burst test itself contributes significantly to the discrepancy. The model parameters were determined by curve fitting uniaxial tensile test data. Inaccuracies in estimating parameters, particularly, the bond shear modulus and bond hardening slope that the inelastic portion is sensitive to (Sinha, 1994), will contribute to the discrepancy. The boundary condition used in simulating the Burst test (rotation unconstrained) to simulate correctly the failure location (high stress) near the center of the sample and not at the boundary, makes the sample more compliant in the model and produce a pressure deflection curve that lags the experimental curve as observed. The unloading part of the curve is highly nonlinear towards the end. This feature is not captured in the model due to not including time-dependent effects. The experimentally obtained pressure–displacement curves show a step increase before 0.1 MPa. The specific cause of this step increase is unclear. In spite of these experimental artifacts, and model parameter estimation inaccuracies, the incremental computational approach has shown promising results in the one-dimensional and twodimensional unloading applications. The model can be used in parametric studies relating microstructural parameters to elastic–plastic response in unloading in applications such as creasing, bending and scoring of paperboard. 6. Summary and conclusions The asymptotic fiber and bond model for the elastic–plastic behavior of paper has been extended to include the unloading behavior and cyclic loading. The incremental form of the asymptotic fiber and bond model and the unloading algorithm have been implemented in ABAQUS finite element code through the UMAT subroutine. Uniaxial loading–unloading tests were carried out with a paper cup stock material (density 720 kg/m3, t = 0.48 mm.). The asymptotic fiber–bond model was used to predict the unloading behavior, specifically, the stress–strain curve and residual plastic strain after unloading. The finite element code with the incremental form of the asymptotic model was also used to predict the unloading curve and the residual strain under uniaxial loading. The discrepancy between the experimental measurements and asymptotic model prediction was between 5% and 13% in MD tests and between 1% and 2% in CD tests. The finite element solution agrees well with the numerical solution of the asymptotic model. Mullen burst tester was used to obtain the pressure versus central deflection in a biaxial loading– unloading test for the cup stock paper material. The finite element analysis with the incremental form of the asymptotic fiber and bond material model shows some promise with the methodology, but also shows significant discrepancies with experimental results. The discrepancies are attributed to poor reproducibility of the Mullen burst test, inaccuracies in the model parameter estimation, and boundary conditions in FEA simulation, and not accounting for time-dependent effects in the model. However, the model results themselves were consistent among each other when tested at different maximum pressures. In spite of this, the finite element implementation of the incremental form of the asymptotic fiber and bond model, together with the unloading algorithm, can be used for modeling complex geometries subjected to loading and unloading.

M.K. Ramasubramanian, Y. Wang / International Journal of Solids and Structures 44 (2007) 7615–7632

7631

Acknowledgements This research was supported in part by unrestricted a gift from James River Corporation (now GeorgiaPacific Corporation). The authors thank Professor Richard Perkins of Syracuse University for insightful discussions and suggestions throughout the course of this study, and his valued friendship. References ABAQUSTM, 1998. User-defined mechanical material behavior, Section 12.8.1, 25.2.30, Version 5.8, Hibbitt, Karlsson & Sorensen, Pawtucket, RI, USA. Batten, G.L., Nissan, A.H., 1987. Unified theory of the mechanical properties of paper and other H-bond-dominated solids 3. Tappi Journal 70 (11), 137–140. Baum, G.A., Brenna, D.C., Habeger, C.C., 1981. Orthotropic elastic constants of paper. Tappi 64 (8), 97–101. Bronkhorst, C.A., 2003. Modelling paper as a two-dimensional elastic–plastic stochastic network. Int. J. Solids Struct. 40, 5441–5454. Castro, J., Ostoja-Starzewski, M., 2003. Elasto-plasticity of paper. Int. J. Plast. 19, 2083–2098. Corte, H. Kallmes, O. J. 1962. Statistical geometry of a fibrous network. In: Bolam, F. (Ed.), Formation and structure of Paper, vol. 1, Tech. Sect. Brit. Paper and Board Makers Association, London, pp. 13. Cox, H.L., 1952. The elasticity and strength of paper and other orthotropic fibrous materials. Brit. J. Appl. Phys. 3, 72–79. Doghri, I., Tinel, L., 2006. Micromechanics of inelastic composites with misaligned inclusions: Numerical treatment of orientation. Comput. Methods Appl. Mech. Eng. 195, 1387–1406. Furukawa, I., Mark, R.E., Crosby, C.M., Perkins, R.W., 1991. Inelastic behavior of machine-made paper related to its structural changes. Japan Tappi J. 45 (5), 582–590. Heyden, S., Gustafsson, P.J., 2001. Stress–strain performance of paper and fluff by network modeling, In: Proc. 12th Fundamental Research Symposium, Oxford, September 2001, pp. 1385–1401. Johnson, M.W., Urbanik, T.J., 1984. A nonlinear theory for plastic plates with application to characterizing paper properties. J. Appl. Mech. 51 (3), 146–152. Kang, G., Gou, S., Deng, C., 2006. Numerical simulation for uniaxial cyclic deformation of discontinuously reinforced metal matrix composites. Mater. Sci. Eng. A 426, 66–76. Karafillis, A.P., Boyce, M.C., 1993. A general anisotropic yield criterion using bounds and a transformation weighting tensor. J. Mech. Phys. Solids 41 (12), 1859–1882. Komori, T., Itoh, M., 1991. Theory of the general deformation of fiber assemblies. Text. Res. J 61 (10), 588–594. Komori, T., Makishima, K., 1977. Numbers of fiber-to-fiber contacts in general fiber assemblies. Text. Res. J. 47 (1), 13–17. Lu, Wentao, Carlsson, Leif A., de Ruvo, Alf, 1996. Micro-model of paper. Tappi J. 79 (2), 197–205. Makela, P., Ostlund, S., 2003. Orthotropic elastic–plastic material model for paper materials. Int. J. Solids Struct. 40, 5599–5620. Niskanen, K., Karenlampi, P., 1998. In-plane tensile properties, In: Paper Physics, Book 16, Papermaking Science and Technology, A series of 19 Books, Niskanen, K. (Ed.), published by Fapet Oy, PO Box 146, FIN-00171, Helsinki, Finland, 1998, pp. 138–191. Page, D.H., Seth, R.S., De Grace, J., 1979. Elastic modulus of paper (1). controlling mechanisms. Tappi J. 62 (9), 99–102. Perkins, R.W., 1986. Fiber Networks (Models for Predicting Mechanical Behavior), Encyclopedia of Materials Science and Engineering, Bever, M.B. (Ed.)-in-Chief, Pergamon Press, 1986, pp. 1712–1719. Perkins, R.W., 1990. Micromechanics Models for Predicting the Elastic and Strength Behavior of Paper Materials, Materials Interactions Relevant to Pulp, Paper, & Wood Ind. (Caulfield, Passaretti, & Sobczynski, ed.)/Materials Res. Soc. Symp. (San Francisco) Proc. vol. 197, April 18–20, 1990, pp. 99–118. Ramasubramanian, M.K. 1987. Computer simulation of the uniaxial stress–strain behavior of ribbonlike fiber Nonwovens, Syracuse Univ. Ph.D. Thesis, 1987, 203p Ramasubramanian, M.K., Perkins, R.W., 1988. Computer simulation of the uniaxial elastic–plastic behavior of paper. J. Eng. Mater. Technol. 110 (2), 117–123. Ramasubramanian, M. K., Wang, Y.Y., 1999. Constitutive models for paper and other ribbon-like nonwovens – a literature review. In: Perkins, R. (Ed.) Mechanics of Cellulosic Materials, AMD-vol. 231, MD-vol. 85, Am. Soc. Mech. Eng., New York, pp. 31–42. Schulgasser, K., Page, D.H., 1988. The influence of transverse us properties on the in-plane elastic behavior of paper. Composites Sci. Technol. 32, 279–292. Schulgasser, K., Grunseit, Z., 1990. The influence of micro-formation on the inplane elastic behavior of paper. J. Mater. Sci. 25 (5), 2433– 2440. Sinha, S.K., 1994. Viscoplastic Micromechanics and Continuum Models of Fibrous Cellulosic Materials. Syracuse University, 322 p. Sinha, S.K., Perkins, R.W., 1995. Micromechanics constitutive model for use in finite element analysis, In: Proc. 1995, Joint ASME Appl. Mech. Mater. Summer Meeting, Los Angeles, CA, USA, Jun 28–30, 1995. Ostoja-Starzewski, M., Stahl, D.C., 2001. Random fiber networks and special elastic orthotropy of paper. J. Elasticity 60, 131–149, 2000. Stenberg, N., 2003. A model for the through-thickness elastic–plastic behaviour of paper. Int. J. Solids Struct. 40, 7483–7498. Suhling, J.C., 1990. Continuum Models for the Mechanical Response of Paper and Paper Composites: Past, Present, and Future, Materials Interactions Relevant to Pulp, Paper, & Wood Ind. (Caulfield, Passaretti, & Sobczynski, ed.)/Materials Res. Soc. Symp. (San Francisco), Proc. vol. 197, April 18–20, 1990, pp. 245–255.

7632

M.K. Ramasubramanian, Y. Wang / International Journal of Solids and Structures 44 (2007) 7615–7632

Van den Akker, J.A., 1970. Structure and textile characteristics of papers. Tappi 53 (3), 388–400. Wang, Y.Y. 2000. Constitutive Modeling of the Unloading Behavior of Paper Material Using the Asymptotic Fiber and Bond Model, Ph.D. Dissertation, Mechanical and Aerospace Engineering, NC State University, NC, USA. Xia, Q.S., Boyce, M.C., Parks, D.M., 2002. A constitutive model for the anisotropic elastic–plastic deformation of paper and paperboard. Int. J. Solids Struct. 39 (15), 4053–4071. Xie, Z., Gilliksson, M., Hagglund, R., 2002. Determining the elastic constants of paper with optimization methods. Inverse Problems Eng. 10 (5), 393–411.