LAPPEENRANTA UNIVERSITY OF TECHNOLOGY LUT School of Engineering Science
LUT Computational Engineering
Minhaj Zaheer
FINITE ELEMENT ANALYSIS OF PAPERBOARD PACKAGES UNDER COMPRESSIONAL FORCE
Supervisor: Associate Professor Joonas Sorvari
ABSTRACT
Lappeenranta University of Technology LUT School of Engineering Science LUT Computational Engineering Minhaj Zaheer
FINITE ELEMENT ANALYSIS OF PAPERBOARD PACKAGES UNDER COMPRESSIONAL FORCE
Master’s thesis 2016
45 pages, 35 figures, 5 table
Supervisor: Associate Professor Joonas Sorvari
Keywords: Paperboard, finite element analysis, corrugated board, buckling, single ply paperboard, multiply paperboard, paperboard box.
This master's thesis studies paperboard failure under compressional forces with aid of modeling. Compressional modeling of paperboard box was done in ABAQUS and conditions were applied to test the feasibility of model. Box was modeled by applying crease pattern at the edges of box, stresses and elastic limits were observed and compared with the normal (without creasing) paperboard box. Then, eigenvalue analysis was performed on single panel and paperboard box to determine the critical stresses and critical forces. The simulations results were then compared with the theoretical calculations and it was perceived that the results were close enough to calculated values.
ACKNOWLEDGEMENTS
This project was completed by the funding of Lappeenranta University of technology, LUT Savo Sustainable Technologies, Varkaus Unit. I thank to Tuomo Kauranne, Haario Heikki, and Teemu Leppänen who provided me opportunity to take part in the research at computational engineering department.
I especially thank to Joonas Sorvari, Associate Professor for assistance with ABAQUS modules, and comments that greatly improved my research focus. Under his supervision, I was motivated, my approach was focused and I achieved some good results within limited time.
Minhaj Zaheer
Lappeenranta 1.2.2016
TABLE OF CONTENTS
1 Introduction ... 7
1.1 Background ... 7
1.2 Aim and scope of the thesis ... 7
2 Paperboard ... 8
2.1 History ... 8
2.2 Production ... 8
2.3 Paperboards hygroscopic behavior ... 9
2.4 Application field of paperboard ... 10
3 Governing equations ... 11
4 Buckling analysis ... 16
4.1 Eigenvalue buckling analysis ... 16
4.2 Nonlinear eigenvalue buckling analysis ... 17
4.3 Critical load calculation ... 17
5 Modeling and simulation ... 20
5.1 Finite element method ... 20
5.2 Abaqus ... 21
5.3 Part designing and ply system ... 21
5.4 Elastic material ... 23
5.5 Material orientation ... 26
5.6 Meshing ... 27
5.6.1 Element type ... 27
5.7 Different geometries ... 28
5.7.1 Geometries for different mesh sizes ... 28
5.7.2 Geometries for eigenvalue analysis by varying lengths and heights ... 30
5.8 With and without crease ... 32
5.8.1 Geometry ... 32
5.8.2 Creases ... 33
6 Simulation results ... 35
6.1 Results for different mesh size ... 35
6.2 Results for eigenvalue analysis for different lengths and heights ... 37
6.3 Results for with and without hinges ... 39
7 Conclusion ... 43
8 Future work ... 43
References ... 44
LIST OF TABLES
Table 1. Elastic material model constants . ... 24
Table 2. Elasto-plastic material model constants. ... 25
Table 3. Boundary conditions. ... 29
Table 4. Boundary conditions for box with and without hinges ... 34
Table 5. Critical load values for different geometries. ... 39
1 INTRODUCTION
1.1 Background
Most commonly corrugated board is utilized as material of packaging. Compressive loading is the vital case when it comes to loading, for example, when the boxes are loaded on the top of each other while placing in big containers. During massive, long distance, and long term shipments many forces are acting on the boxes, especially compressive forces at the top of box. This might result a collapse that can either by structural or material failure.
In various applications the demands of corrugated board is increasing day by day. In order to make it compatible with the today’s requirements it is necessary to optimize its material and structural properties. Different material and structural properties can be tested via modeling. The testing of the models can be done with different material-, geometry- and load parameter combinations in finite element analyses (FEA). Based upon the simulation results we can easily forecast that either our model will be suitable for real world or not.
1.2 Aim and scope of the thesis
The aim of this master thesis is to perform a detailed analysis of a paperboard failure in ABAQUS when loaded in plane. Paperboard is firstly modeled according to given parameters, at box edges the crease pattern is applied and the results are compared with the normal (without creasing) box. The aim in this case is to enhance the load bearing capacity of box so it can bear more stresses than normal. In the analysis, both local and global stability is considered, especially the effects of local buckling on the board's total load carrying capacity. Linear eigenvalue analysis is performed by applying compressional load at the top surface of box. Critical stress and load values are determine when buckling will start. For different geometries load and stress values are calculated and compared with the theoretical values. The aim is this case is to predict the critical load value for different geometries of box.
2 PAPERBOARD
2.1 History
In industries paperboard is mostly used because it is recyclable and environment friendly. In 105 A.D. first paperboard was manufactured from cellulose fibers from flax, cotton and vegetable sources by China [6]. In England during 1817, the first paperboard carton was produced. Around the 1860s folding cartons were manufactured. They save space during the shipment, and can easily set up when the customer demands. In 1879 mechanical die cutting and creasing of blanks was developed and in 1915 the gable top milk carton was patented [6].
2.2 Production
In these days manufacturing process of paperboard begins with separating cellulose fibers originated from renewable raw material. The structure of sheet or web consist of an interlaced network of fibers. These fiber are separated by chemical or mechanical methods from naturally occurring source of birch, pine or spruce [6]. There are two basic types of fibers from which paperboard is made, mechanical and chemical processed wood fibers.
Both treatments have different properties of pulp that effects the paperboard properties (see Figure 1) [1].
Figure 1. Pulp properties and influence on paperboard properties [1].
Mostly recyclable material or tree pulp is used to manufacture paperboard. Paperboard can be divided into two categories, white and brown paperboard. Both white and brown paperboard are further classified into two classes.
Figure 2. Paper products chart [8].
2.3 Paperboards hygroscopic behavior
Moisture needs to be addressed while considering paperboards because properties of paperboard will change in the presence of moisture. If there is more humidity in the surrounding environment paperboard can gain the moisture easily. Expansion and contraction of paper is directly linked with changes in the moisture content of surroundings and is typically 0.8% in MD and 1.6% in CD [5]. The moisture content of paperboard are
Tree/Waste paper
Paper Paperboard
Brown (container board) Pulp
Brown (65%
hardwood and 35% softwood)
White (95%
-100%
softwood) White
(paperboard package) Boxboard
SBS (Solid Bleach Sulfate)
dependent on climatic conditions. While for the testing of paper international standards are 23 Celsius and 50% R.H for temperature and humidity [5]. If moisture contents are uneven in the layers of paperboard than paperboard warping can be an issue. Hysteresis is shown by paper, for example, the equilibrium moisture content is slightly different depending on whether it is approached from higher or lower initial moisture content as Figure 3 is demonstrating [5].
Figure 3. Moisture content and relative humidity [5].
2.4 Application field of paperboard
As the paperboard is usually light in weight, it can easily be customized according to the demand. Most of the industries prefer to use paperboard for packaging of products. Wet products can be packed after applying some chemicals on the surface as in Figure 4. The boxes can be recycled after the purpose is done. Apart from that it is also useful in packing of food, electronics and mechanical parts.
Figure 4. Corrugated board packing for fish [7].
3 GOVERNING EQUATIONS
Stress is defined as force per unit area of a body as the area size tends to be zero. It is represented by sign σ and units N/m2. Load bearing structures undergoes stresses, displacement, and strains when external loads are acting on them. Some cases are looking for intensity of stresses and their distribution while others concerns about the deformations affecting the distribution of stresses or the stiffness. These mechanical behaviors follows physical laws like energy conservation, balance of momentum, conservation of mass, and principle governing thermodynamics [9].
There are nine stress components, see Figure 5. Component of stress tensor 𝜎𝑖𝑗 can be given in a matrix form as
𝜎𝑖𝑗=[
σ11 σ12 σ13 σ21 σ22 σ23
σ31 σ32 σ33]
(1)
In the symbol 𝜎𝑖𝑗, 𝑖 and 𝑗 are to be taken as 1, 2 and 3 as we can see from the matrix on right hand side. σ symbol also represents the stress tensor. When the letter in the subscript are repeating as in 𝜎𝑖𝑖 then it represents the sum of three terms, for example
𝜎𝑖𝑖 = 𝜎11+ 𝜎22 + 𝜎33 (2)
Figure 5. Stress components [9].
According to Findley excessive distortion of a solid body under influence of external loading, justify three governing equations: (1) equilibrium; (2) kinematics; and (3) material constitution. Material constitutive equation is the only one that depends on material type. It represent material properties and mechanical behavior [9].
Paperboard follows static equilibrium in which resultant of all forces is zero. Along x- direction the equilibrium can be sum up as
0 )
(
) (
) (
3 2 1 1 1 2 31 1 2 3 3 31 31
3 1 21 3 1 2 2 21 21
3 2 11 3 2 1 1 11 11
dx dx dx F dx dx dx
dx x dx
dx dx dx
dx x dx
dx dx dx
dx x dx
(3)
where 𝜎11 represents stress component along x-direction, 𝜎21 is the stress along x-direction with respect to y component of stress and 𝜎31 is stress along x-direction with respect to z component of stress. Similarly, equation along y- and z-directions can be obtained. 𝐹1, 𝐹2, and 𝐹3 are components of body forces per unit volume in x-, y-, and z-directions. By dividing equations along x-, y-, and z-directions with
dx
1dx
2dx
3, it will give equation on each direction of coordinate axes [9].. 0
, 0
, 0
3 3 33 2
23 1
13
2 3 32 2
22 1
12
1 3 31 2
21 1
11
x F x
x
x F x
x
x F x
x
(4)
Equation (4) can be express in tensor form as
0
j i
ij F
x
) 3 , 2 , 1 ,
(i j (5)
Linear strain is usually defined as change in length to its original length or it is also the measure of displacement between the particles while deformation relative to a reference length. It is a unitless quantity. Engineering shear strain is measure as the change in angle of two lines intersecting at a point where the lines are at right angles [9]. Like stresses strain
has also nominal and shear components. Normally strain is caused by the stress for the materials obeying Hooke's law and that will cause the dilations. Figure 6 is showing how the deformation will occur in two-dimensional rectangular material after stress application.
Figure 6. Deformation after stresses [15].
The normal strain in the x-, y- and z-directions can be defined as 𝜀𝑥= 𝜕𝑢𝑥
𝜕𝑥 , 𝜀𝑦 = 𝜕𝑢𝑦
𝜕𝑦 , 𝜀𝑧 =𝜕𝑢𝑧
𝜕𝑧
where as 𝑢𝑥, 𝑢𝑦, and 𝑢𝑧 represents the deformation along x-, y- and z-directions.
For an orthotropic linearly elastic material the relation between stresses and strains is given by
yz xz xy zz yy xx
=
yz xz
xy z
y yz x
xz
z zy y
x xy
z zx y
yx x
G G
G E
E v E
v
E v E
E v
E v E
v E
0 1 0
0 0
0
1 0 0
0 0
0
0 1 0
0 0
0
0 0
1 0
0 0
1 0
0 0
1 0
yz xz xy zz yy xx
(6)
where 𝜎𝑥𝑥, 𝜎𝑦𝑦, 𝜎𝑧𝑧, 𝜏𝑥𝑦, 𝜏𝑥𝑧 and 𝜏𝑦𝑧 represent components of stress tensor (N/m2) and Ex, Ey, Ez are elastic modulus in x-, y-, and z-directions. Gxy, Gxz, and Gyz are shear modulus expressed in N/m2 and εxx, εyy, εzz, xy,
xz, and yz are components of strain tensor [16].There are six Poisson ratios υxy, υyx, υxz, υzx, υyz and υzx that needs to be determine due to anisotropy. A material is in plane stress state if the stress vector is zero across a particular surface. Plane stress occurs in thin flat plates that are acted upon only by load forces that are parallel to them. Generally for the stress analysis a gently curved thin plate may also be assumed to have plane stress [18]. In plane stress state (6) can be reduced to the x- and y- axes and stress-strain behavior can be given by [10]
s y x
=
s y
x xy
y yx x
E E
E v
E v E
0 1 0
1 0 1 0
s y x
(7)
The paperboard material model was comprised of elastic as well as plastic material properties. Total strain 𝜀 consist of both elastic and plastic parameters during plastic loading.
𝜀 = 𝜀𝑝+ 𝜀𝑒 (8)
𝜀𝑒 = 𝐷−1𝜎 (9) where D is the elastic moduli matrix and 𝜎 is the stress rate.
The plastic behavior of paperboard material model can be calculated with the help of Hill’s criteria of plasticity. According to Hill’s criteria, plastic strain rate 𝑒̇𝑝 is given by [20]
𝑒̇𝑝 = 𝜆𝜕𝑓
𝜕𝜎 (10) where 𝜆 denotes plastic multiplier, 𝜎 represents stress tensor and 𝑓 is yield function which can be given as [20]
𝑓 = 𝐹(𝜎𝑦𝑦− 𝜎𝑧𝑧)2+ 𝐺(𝜎𝑧𝑧− 𝜎𝑥𝑥)2+ 𝐻(𝜎𝑥𝑥 − 𝜎𝑦𝑦)2+ 2𝐿𝜎𝑦𝑧2 + 2𝑀𝜎𝑧𝑥2 + 2𝑁𝜎𝑥𝑦2 − 1 F, G, H, L, M and N are the material constants. These constants can be determined experimentally and they can be related to yield stress as [20]
2𝐹 = 1
(𝜎𝑦𝑦0 )2+ 1
(𝜎𝑧𝑧0)2− 1
(𝜎𝑥𝑥0 )2 2𝐺 = 1
(𝜎𝑧𝑧0)2+ 1
(𝜎𝑥𝑥0 )2− 1
(𝜎𝑦𝑦0 )2 2𝐻 = 1
(𝜎𝑥𝑥0 )2+ 1
(𝜎𝑦𝑦0 )2− 1
(𝜎𝑧𝑧0)2 (11) 2𝐿 = 1
(𝜎𝑦𝑧0 )2 2𝑀 = 1
(𝜎𝑧𝑥0 )2 2𝑁 = 1
(𝜎𝑥𝑦0 )2
where 𝜎𝑥𝑥0 , 𝜎𝑦𝑦0 , 𝜎𝑧𝑧0 are yield stresses in x-, y-, and z-directions.
4 BUCKLING ANALYSIS
Buckling analysis is used when we have to determine the critical load values at which the instability of a structure will occur and it also helps in investigating the different collapse mode shapes. FEM buckling analysis of paperboard requires material constants and geometrical model. There are two types of buckling analysis (1) nonlinear buckling analysis, (2) eigenvalue (or linear) buckling analysis. Both methods yields various different results [3].
4.1 Eigenvalue buckling analysis
In this thesis we will only look for eigenvalue buckling analysis as our structure is elastic one. For most of elastic structures this eigenvalue buckling analysis is used to predict the theoretically buckling strength (the bifurcation point) [3]. This method is a linear perturbation procedure and it is used to determine the structure imperfection sensitivity.
Those models in which there are supporting structures this method cannot be used for investigating the critical loads values. Figure 7 is showing bifurcation point (buckling point) during analysis. This analysis is mainly used for theoretical approaches. Nevertheless, real- world structures are not able to have their theoretical elastic buckling strength due to nonlinearities and limitations [3].
Figure 7. Linear (eigenvalue) buckling curve [3].
4.2 Nonlinear eigenvalue buckling analysis
Nonlinear buckling analysis sometimes give more feasible results as compared to the linear eigenvalue analysis because in this case nonlinear deformation is taken into account. It operates in very simple way: load is increased slowly until it reaches to a level where it is found that the structure is unstable, and after that point with small increment in load very large deflections will occur. This method is used for example for the structures exhibiting the plastic properties. In many cases for elastic material this type of analysis is normally not used.
As shown in Figure 8 the nonlinear instability region is the region in which the equilibrium path goes from one stable point (1) to another new stable point (2). The nonlinear behavior places the critical limit load at point 1 equal to that at point 2, but the load limit corresponds to a new structural shape. The second stable point along the equilibrium path occurs after a large displacement of the structure. The nonlinear buckling nature of shallow trusses includes a post-buckling instability region along the equilibrium path [4].
Figure 8. Nonlinear buckling curve [4].
4.3 Critical load calculation
Many exact solutions for critical loads have been formulated in Timoshenko monographs [10]. Lekhnitskii and Whitney [10] compiled many exact and approximated solutions for anisotropic plates. Moreover these solutions are based on assumptions and far for real condition because load was considered to be uniform at the ends of plates [10]. After that
Yoon and Lee did buckling analysis and considered edges to be stiffer. They both studied short term behavior of orthotropic plates. The objective of our work is to calculate the approximated values for critical stresses and forces for different geometrical shapes of anisotropic plates while in-plane load is applied with selected boundary conditions [10].
Assume a rectangular plate where ‘a’ is length, ‘b’ is width, and ‘h’ is the thickness of orthotropic plate with x and y material directions as shown in Figure 9.
Figure 9. Schematics of an orthotropic plate under in-plane uniaxial compression [10].
Ex, Ey are young’s moduli in x- and y-directions. υx, and υy are Poisson ratio’s in x- and y-directions, where Es is the shear modulus. When the plates are bent the equation for displacement w of plate in the z-direction is given by [10]
𝜕4𝑤
𝜕𝑥4 𝐷11+ 2(𝐷12+ 2𝐷66) 𝜕4𝑤
𝜕𝑥2𝜕𝑦2+ 𝐷22𝜕4𝑤
𝜕𝑦4 = 𝐹(𝑥, 𝑦)
(12) where F(x,y) depends on the type of considered problem. Dij are flexural rigidities of plate
given by 𝐷11= 𝐸𝑥ℎ3
12(1−𝜈𝑥𝜈𝑦)
(13) 𝐷22 = 𝐸𝑦ℎ3
12(1−𝜈𝑥𝜈𝑦) (14) 𝐷12= 𝜈𝑦𝐸𝑥ℎ
3
12(1−𝜈𝑥𝜈𝑦)
(15) 𝐷66 =𝐸𝑠ℎ3
12 (16) For plate under in-plane compression Nx in x-direction is shown in Figure 9. F(x,y) is given by
𝐹(𝑥, 𝑦) = −𝑁𝑥𝜕2𝑤
𝜕𝑥2 (17) From equation (12) we can have
𝜕4𝑤
𝜕𝑥4 + 2Ƞ𝜆1/2 𝜕4𝑤
𝜕𝑥2𝜕𝑦2+ 𝜆𝜕4𝑤
𝜕𝑦4 =𝐹(𝑥,𝑦)
𝐷11
(18)
where 𝜆 =𝐷22
𝐷11
,
Ƞ = 𝐷12+2𝐷66√𝐷11𝐷22
(19a)
For homogenous, orthotropic plate (19a) can written as 𝜆 =𝐸𝑦
𝐸𝑥 , Ƞ =2𝐸𝑠(1−𝜈𝑥𝜈𝑦)
√𝐸𝑥𝐸𝑦 + √𝜈𝑥𝜈𝑦
(19b) The critical buckling stress 𝜎𝑐𝑟 is dependent on geometry and material of plate and buckling mode. Critical buckling is unitless and can be expressed as
𝜎𝑐𝑟
𝜎0 = ∅(𝑅, Ƞ, 𝑚) (20) 𝜎0 = 𝜋2√𝐷11𝐷22
𝑏2ℎ (21)
where 𝜎0 is the reference buckling stress. In Equation (20) R is the aspected ratio and m is the number of half waves in loading direction
𝑅 = 𝜆1/4𝑎/𝑏 (22) In our calculations we considered all edges are simply supported for plate so the boundary conditions are like
𝑤⎹𝑥=0,𝑎 = 0, 𝑤⎹𝑦=0,𝑏= 0
𝜕2𝑤
𝜕𝑥2 ⎸𝑥=0,𝑎 = 0 ,𝜕2𝑤
𝜕𝑦2 ⎸𝑦=0,𝑏 = 0 (23) The exact solution for critical buckling stress 𝜎𝑐𝑟 given by
𝜎𝑐𝑟 𝜎0 =𝑚2
𝑅2 + 𝑅2
𝑚2+ 2Ƞ (24) Plate will buckle in m half wave if this satisfies [10]
√𝑚(𝑚 − 1) <= R <= √𝑚(𝑚 + 1) (25)
5 MODELING AND SIMULATION
5.1 Finite element method
Finite element method is a numerical technique in which a complicated domain can be sub- divided into a series of smaller regions in which the differential equations are approximately solved. By assembling the set of equations for each region, the behavior over the entire problem domain is determined [19]. The FEM is used for finding the approximate solutions of partial differential equations (PDE) as well as integral equations [6]. It divide the structure into number of finite pieces called elements and these elements are reconnected by nodes.
The solution can be calculated with the help of computation tools such as Abaqus.
Figure 10. Flowchart of FEA.
Material properties
Material law (Hooke law &
Hill's citeria)
Finite element model
Simulated displacement Boundary
conditions
Convergence
Optimal solution Yes
No Creating
parts
5.2 Abaqus
Abaqus is finite element analysis software that use mathematical techniques to give solution for real world engineering problems. Nonlinear and linear problems with innovative designs, assemblies and manufacturing can be solved in abaqus. Many engineering companies are using this software to solve complex problems. Abaqus was designed for nonlinear problems and it can utilize for many materials, for example, elastic, plastic, and viscous materials.
Industries like automobile, packaging, aerospace, chemical, and petroleum industries are using this software for their complex engineering problems. There are three stages of finite element analysis
1. Creating input
2. Processing or FE analyzer
3. Post processing (report, image and animation)
Abaqus generates many files in its directory when the simulation is running. They are DAT, INP, JNL, CAE, IPM, PRT, etc. Each file is having some information regarding the simulations, DAT file includes warnings and error information after the simulation is simulated. INP file is the input file generated by Abaqus before simulation and it contains detail information about the steps done during the formation of simulation. We can get exact information about number of elements, number of nodes, number of internal nodes and job summary.
Through Abaqus we can do preprocessing, post processing, and investigation of different complex problems at various time steps. We can also import structures from different softwares like Solid Works, Pro E, Creo, etc. Current research is done by using Abaqus 6.14, no other CAD software was used. Drawing process was performed in Abaqus user interface.
5.3 Part designing and ply system
In our work Abaqus user interface is used to design the complete geometries. This software can provide us different options for creating deformable, analytical rigid and eulerian types of element in geometry. Abaqus can combine different types of surfaces in one. A rigid body usually contain thousands of nodes and these nodes can be controlled by allocating a single
reference node along the body. While running simulation if load is applied to the single reference node it appears to act on all the nodes associated with the reference node.
Discrete part behaves same as the analytical rigid part. Analytical rigid body shape is not based on arbitrary parameters but the geometry is restricted to sketch lines, arcs and parabolas. On the other hand, deformable part can be shaped as axisymmetric, 2-D and 3-D part, required type of geometry can also be selected with basic features and shape. From the other CAD software the geometry can be imported to the Abaqus. This software creates deformable part that can deform under different loading conditions. Load can be mechanical, thermal or electrical depending on the problem.
In this thesis, full assembly is designed and simulated in user interface of Abaqus.
Techniques were derived to reduce the simulation time and enhance precision and accuracy of results. Firstly, a box having single ply was modelled with thickness of 0.4 mm. Thickness usually deals with the how much amount of fibers is put in per square meter of a box. This structure is made by using the deformable material elements and these elements has displacement ability when load is being applied at the top surface of box.
Secondly, multi-ply cardboard boxes are designed and the buckling load is calculated. These boxes has more folding and creasing performance as compared to the single ply boxes. In situations where the same sort of pulp is being utilized as a part of a few layers, every different layer is dealt with and shaped separately to create the highest possible quality [13].
Multi-ply system will increase the stiffness and weight of product. Another benefit of multi- ply system is that the cost of fiber is reduced while manufacturing product and in each layer controlled amount of fibers are added to get the desired characteristics of product. Figure 11 is showing different number of plies used for the manufacturing of corrugated board.
Figure 11. Different number of plies [12].
5.4 Elastic material
Performance and appearance are two distinguishing properties of paperboard. Printability, whiteness, ink absorption and rub resistance are classified as appearance properties whereas the performance properties are related to the physical characteristics of the paperboard.
Material properties calculations should be accurate and precise. Errors in properties calculations lead to several defects in the end product [6].
In some situations paperboard can be assumed to be elasto-plastic material. In these cases the material deformation is directly related to the applied stress. Initially after the stress application there comes a stage where paperboard reaches its elastic limit as shown in Figure 12 and if stress is further applied the material deformation enters into plastic region where the material is permanently deformed and cannot go back to its original position.
Figure 12. Elasto-plastic behavior of paperboard [14].
In our work, paperboard was considered as elastic and elasto-plastic material. The elastic material properties are given in Table 1. The material model was adopted from literature presented by Huang [16]. Elasto-plastic material properties was considered for other part (with and without creasing pattern) as given in Table 2.
Constants Top ply Middle ply Bottom ply
E
xx(MPa) 5920 3202 8760Eyy(MPa) 2670 1233 3030
Ezz(MPa) 228 160 120
Gxy(MPa) 1431 638 1617
G
xz (MPa) 60 30 60Gyz(MPa) 60 30 60
xy 0.45 0.47 0.51
xz 0 0 0yz 0 0 0
Table 1. Elastic material model constants [16].
Table 2. Elasto-plastic material model constants.
Constants Ply
E
xx(MPa) 3203Eyy(MPa) 1233
Ezz(MPa) 160
Gxy(MPa) 638
G
xz(MPa) 30Gyz(MPa) 30
xy 0.47
xz 0yz 0
H
(MPa) 22890
xx(MPa) 400
yy(MPa) 14.2
0
zz(MPa) 14.20
xy(MPa) 28
0
xz(MPa) 4.80
yz(MPa) 3.48
5.5 Material orientation
As paperboard is highly anisotropic material and because of complex fiber orientation it usually exhibits nonlinear mechanical behavior. Material orientation is a critical parameter in order to compute accurate results. Figure 13 shows paperboard material orientation in x-, y-, and z-direction. While modeling in Abaqus, direction 1 (blue) correspond to machine direction (MD), direction 2 (yellow) represents cross direction (CD) and direction 3 (red) as z-direction (see Figure 14).
Figure 13. Material orientation in x-, y- and z-direction [6].
Figure 14. Material orientation in Abaqus.
5.6 Meshing
Mesh tool in Abaqus basically divides the whole assembly into various nodes and small elements as we specify in software. Abaqus also provide mesh control options like tri, quad, and quad-dominated. Stack can also be assign in mesh part and there are many meshing techniques for complex structures. We can also check the mesh quality after assigning the mesh to a geometry by using ‘verifying mesh’ option.
Mesh control provide flexibility in element shape, structure, technique, algorithm and adaptive remeshing rules. After assigning number of seeds and the element size to a geometry we can create the orphan mesh by using create mesh part. Seeding also control the mesh quality at the edges of structures.
Figure 15. A model with biased seeding.
5.6.1 Element type
Mesh element type controls the element shape that helps in convergence of solution. Right element type reduces the analysis time and enhance accuracy. Abaqus provide variety of element types. The main stress/displacement elements available in ABAQUS include the 4- node linear tetrahedron, the 6-node linear triangular prism, the 8-node linear brick, the 10- node quadratic tetrahedron, the 15-node quadratic triangle and the 20-node quadratic brick [17]. Each element is having three degrees of freedom per node. Figure 16 is showing these elements.
Figure 16. Three-dimensional solid elements and shapes [17].
By default Abaqus allocate certain element type to mesh region. The element shape is relate by its characteristics of element type such as shell part after meshing, assigned by triangular and quadrilateral element type as an auto generated tool in Abaqus. Element type of mesh region can be interpreted by equivalent element shape. For example, mesh region of shell can form triangular elements while ignoring the quadrilateral element type.
5.7 Different geometries
5.7.1 Geometries for different mesh sizes
In this thesis, at first the maximum load bearing capacity of a box is determined for a box having dimensions 200 mm × 200 mm × 200 mm. Figure 17 is showing the length and width of box sketched in the Abaqus user interface. 𝐹𝑚𝑎𝑥 (maximum force for each mesh size) value is obtained for each mesh size. Figure 18 shows the different mesh sizes assigned to the box. Mesh seed of 2, 5, 10 and 20 mm was assigned to box. Tri mesh was allocated to them to see the behavior while applying load.
Material orientation and material properties were assigned according to the Figure 14 and Table 2. Elasto-plastic material was considered for this simulation. It will go to its linear elastic limit first and after that when more force will act it will enter into its plastic region.
Boundary conditions are specified in Table 3, the lower part was fixed completely and the displacement of unit 10 mm was applied at the top nodes of box.
Figure 17. Geometry for box (length and width).
Parts Boundary conditions
Top of box U1= 0
U2= - 10 mm U3= 0 UR1=0 UR2=0 UR3=0
Bottom of box U1= 0
U2= 0 U3= 0 UR1=0 UR2=0 UR3=0 Table 3. Boundary conditions.
Figure 18. Different mesh seed size for box.
5.7.2 Geometries for eigenvalue analysis by varying lengths and heights
Geometries of different length and height are modeled in the user interface of Abaqus, Figure 20 is showing those modeled structures by varying heights and lengths. Units for height, length and width are in mm. The thickness of the wall was divided into three plies; top, bottom, and middle ply. Each has its own material properties as mentioned in Table 1.
Middle ply is considered twice thicker than the top and bottom. Material orientation is same shown in Figure 14. In eigenvalue analysis, linear perturbation step is used, number of increments and time step was specified according to the requirement. Unit load is applied at the top of the box as shown in Figure 19 and the lower surface is kept fixed.
mesh seed size 2 mm mesh seed size 5 mm
mesh seed size 10 mm mesh seed size 20 mm
Figure 19. Load at box using constraints.
Varying length Varying height
length=100,height=100,Width=100 length=100,height=100,Width=100
length=200,height=100,Width=100 length=100,height=200,Width=100
length=300,height=100,Width=100 length=100,height=300,Width=100 Figure 20. Different sizes for box.
5.8 With and without crease 5.8.1 Geometry
In this part four shell planar plates of 30 mm × 30 mm were formed in Abaqus user interface and the thickness of wall was kept 0.44 mm (Figure 21). Material orientation and properties were assigned according to Figure 14 and Table 2. Elasto-plastic material model was considered. A plate of 40 mm × 40 mm was placed under the box as Figure 22 is illustrating.
Figure 21. Dimension for shell plate.
Figure 22. Geometry with hinges.
5.8.2 Creases
Crease is usually define as a separation line that distinguish two folded surfaces of paperboard. Crease modeling is bit challenging, hinge joint and the crease line degree of freedom is almost similar. Crease line offers simple three displacement and three rotational degrees of freedom in total six degrees of freedom. However, the five degrees of freedom are fixed between two nodes and they are only allowed to have one degree of freedom (degree of freedom 4) which means the displacement can be only in one specified direction and allow only rotation along crease line [6], see in Figure 23.
Figure 23. Degrees of freedom illustration [6].
In Abaqus, connector tool give option to have hinge joint that is normally connected between the two parallel nodes, but the joint can move in any direction. This is the closest way to correlate paperboard folding in reality with the crease behavior. Performance of crease lines are depend on bending stiffness, angle of creasing, stress properties of paperboard and structure geometry. Figure 24 demonstrate how the crease pattern is applied to edges of plate and Table 4 represents the boundary conditions applied at the top and bottom surface.
Figure 24. Crease pattern on box.
Table 4. Boundary conditions for box with and without hinges.
Parts Boundary conditions
Top of box U1= 0
U2= - 1.5 mm U3= 0 UR1=0 UR2=0 UR3=0
Bottom of box U1= 0
U2= 0 U3= 0 UR1=0 UR2=0 UR3=0
6 SIMULATION RESULTS
6.1 Results for different mesh size
Figure 25 shows the occurrence of deformations after the displacement of 10 mm units during the simulation. When the element size is small forces acting at the boxes were low and most of deformation was at center of the box. As the element size increase more forces start acting on the box especially at the edges. From the simulation result of element size 20 mm it is observed that most forces were acting at the edges the box. Figure 26 is demonstrating that for lower element size the forces acting at the box are lower as compared to the bigger element size. Figure 27 shows maximum force as a function of number of elements. The graph showed inverse relation between 𝐹𝑚𝑎𝑥 and element number. Number of elements increases with decrease in force value. Furthermore, as the element size is increasing the solution converges. So the mesh seed size should be small in order to attain the accurate results.
mesh seed size 2 mm mesh seed size 5 mm
mesh seed size 10 mm mesh seed size 20 mm
Figure 25. Simulation results of different mesh.
mesh seed size 2 mm mesh seed size 5 mm
mesh seed size 10 mm mesh seed size 20 mm
Figure 26. Force vs. time graph for different mesh.
Figure 27. Graph of 𝐹𝑚𝑎𝑥vs. element number.
6.2 Results for eigenvalue analysis for different lengths and heights
This part contains the eigenvalue analysis performed on boxes having different lengths and heights. The purpose of these analysis was to get the value at which the box will first buckle and to check the stresses and strain in the box. Figure 28 shows the deformation in box at it first buckling point and eigenvalue at which the box will buckle. Once we know the eigen- value we can calculate the critical force by multiplying the applied load in simulation with the eigenvalue when it first buckle. In these simulations unit load was applied the top surface of the box so critical load value is the same as eigenvalue. Table 5 illustrates few critical loads value for different geometries.
Critical load values were also theoretically calculated for single plate considering all edges simply supported as mentioned in 4.3. These theoretical values for critical loads were plotted against the ratio of length and width (a/b). The critical load values from simulation were also plotted against the ratio of a and b on the same graph (Figure 29). Height was considered to be fix. It can be seen that trend of theoretical values was almost same as the simulation values. Critical force value decreases with increase in ratio of a and b.
0 100 200 300 400 500 600 700 800
0 20000 40000 60000 80000 100000 120000
Fmax (N)
Element number
Varying length Varying height
length=100,height=100,Width=100 length=100,height=100,Width=100
length=200,height=100,Width=100 length=100,height=200,Width=100
length=300,height=100,Width=100 length=100,height=300,Width=100
Figure 28. Eigenvalue analysis for different geometries.
Table 5. Critical load values for different geometries.
Figure 29. Eigenvalue analysis for different geometries.
6.3 Results for with and without hinges
Simulation of the box having hinges at the edges was tested to check the load bearing capacity of the box, stress on the box and elastic limit of structure. Mesh seed size was 1 mm. These results were then compared with the box having same parameters but without creasing pattern. Top of the box was displaced 1.5 mm which corresponds 5% of the height of the box. From the Figure 30 we can see that the box with hinges have lower maximum stresses at the edges than the box without hinges (Figure 31). Values in the green region are higher for hinges than for the without hinges (Figure 32 & 33).
0 20 40 60 80 100 120 140
0 0.5 1 1.5 2 2.5 3 3.5
Fcr (N)
a/b
Fcr vs a/b
Theoratical Values Simulation values
(mm) (mm) (mm) (N)
100 100 100 51.4
100 100 200 41.7
100 100 300 39
150 100 100 44
200 100 100 43.5
300 100 100 51
Figure 30. Von Mises stress for box with hinges.
Figure 31. Von Mises stress for box without hinges.
Figure 32. In-plane principal strain for box with hinges.
Figure 33. In-plane principal strain for box without hinges.
Graph (Figure 34) illustrates the point at which force acting on creasing pattern box is maximum. After peak the force decrease rapidly and then it remain constant. Figure 35 compares stress vs. strain values for both with and without hinges; for box with hinges a bit higher stress values are obtained as the hinges made the stucture more stiffer than the simple one.
Figure 34. Force vs. time for creasing pattern.
Figure 35. Stress vs. strain for with and without hinges geometry.
0 100 200 300 400 500 600
0 0.013244448 0.048244447 0.08324445 0.118244447 0.153244451 0.188244447 0.223244444 0.258244455 0.293244451 0.328244448 0.363244444 0.398244441 0.433244437 0.468244433 0.50324446 0.538244426 0.573244452 0.608244419 0.643244445 0.678244472 0.713244438 0.748244464 0.783244431 0.818244457 0.853244424 0.88824445 0.923244476 0.958244443 0.993244469
Force (N)
Time (s) Force vs Time
-5 0 5 10 15 20
1 25 49 73 97 121 145 169 193 217 241 265 289 313 337 361 385 409 433 457 481 505 529 553 577 601 625 649 673 697 721 745 769 793
stress (MPa)
Strain (×10-5) Without Hinge With Hinge
7 CONCLUSION
Research was done at first to see behaviour of force acting on the box with varying the mesh size. After the displacement of 10 units in simulation the results were recorded and it was observed that with increase in element size increases the force. Graph between the Fmax and element number also revealed that solution converges with decrease in element size.
Eigenvalue analysis was done on different geometries to get a critical load value at which the box will first buckle and how the deformation will occur. These critical load value were then compared with the theoratical values and the trend was almost same in both cases.
In the last part a research was done to studyhow hinges effect on simulated results. Box with the hinges was compared with box without hinges. The results showed that hinges made the box stiffer and affected almost each studied parameter of box. Simulation and mathematical modeling revealed that material model, anisotropy and boundary conditions play a vital role in testing of paperboard boxes parameters. Modeling was a challenging task, the obstacle like hinges operation, stability of structure, excessive deformation and simulation time were faced in these processes.
8 FUTURE WORK
More research is required on application of creasing pattern to the box so the load bearing capacity of box could be increased. Some other approaches could be made apart from creasing that will be useful in enhancing the load bearing capacity of box when they are placed in shipment. More work can be done in choosing the material parameters of box to broaden its application in different fields.
REFERENCES
[1] "Baseboard physical properties," Iggesund Paperboard, [Online]. Available:
https://www.iggesund.com/globalassets/iggesund-documents/paperboard-
documents/iggesund-anchor-material/reference-manual/reference-manaul-baseboard- physical-properties.pdf.
[2] A. Emblem and H. Emblem, Packaging Technology: Fundamentals, Materials and Processes, Woodhead publishing, 2012, p. 257.
[3] J. C. Amigo, "Stiffness Design of Paperboard Packages using the Finite Element Method," Stockholm, 2012.
[4] "Eco-packaging for the fish industry," Stora Enso, 29 1 2016. [Online]. Available:
http://www.storaenso.com/newsandmedia/eco-packaging-for-the-fish-industry.
[5] A. Urmanbetova, "Terminology On Paper & Pulp:Types of Paper and Containerboard,Containerboard Grades and Tests," 2001.
[6] J. S. Lai, K. Onaran and W. N. Findly, Creep And Relaxation Of Nonlinear Viscoelastic Materials., New York: Dover publications, pp. 22-48.
[7] I. Hwang and . J. S. Lee, "Buckling of Orthotropic Plates under Various Inplane Loads," KSCE Journal of Civil Engineering, Vols. Vol. 10, No. 5, 2006.
[8] "ANSYS Structural Analysis Guide," [Online]. Available:
http://www.ansys.stuba.sk/html/guide_55/g-str/GSTR7.htm.
[9] G. A. Hrinda, "Snap-Through Instability Patterns in Truss Structures," Hampton, Virginia 23831.
[10] "Infinitesimal strain theory," Wikimedia Foundation, Inc., [Online]. Available:
https://en.wikipedia.org/wiki/Infinitesimal_strain_theory.
[11] R. Cook, Finite element modeling for stress analysis, New York: Wiley, 1995, p. 1.
[12] "All About Box - A To Z Janta Packaging, Mumbai, India," A To Z Janta Packaging, [Online]. Available: http://www.corrugatedboxess.com/all-about-box.html.
[13] "Paperboard the Iggesund Way," [Online]. Available:
https://www.holmen.com/globalassets/holmen-documents/iggesund/iam- files/iggesund_ptiw_en2010.pdf.
[14] "Defined physical properties," [Online]. Available:
https://www.iggesund.com/en/knowledge/the-reference-manual/baseboard-physical- properties/Defined-physical-properties/.
[15] "Infinitesimal strain theory," [Online]. Available:
https://en.wikipedia.org/wiki/Infinitesimal_strain_theory..
[16] H. Huang and M. Nygårds, "Numerical Investigation of Paperboard Forming," Nordic Pulp and Paper Research Journal Vol 27 no.2, pp. 211-215, 2012.
[17] S. Johnson, "Comparison of Nonlinear Finite Element," Illinois, 2006.
[18] "Plane stress," Wikimedia Foundation, Inc., [Online]. Available:
https://en.wikipedia.org/wiki/Plane_stress.
[19] "Finite element modeling," Center of Computer Aided Research and Development, [Online]. Available: http://www.zarbakhsh.com/research/computer-simulation/finite- element-modeling/#.V4UnKUbJbM4.
[20] R. D. Broast and P. H. Feenstra, "Studies in anisotropic plasticity with reference to the Hill criterion," International Journal for Numerical Methods in Engineering, vol. 29, pp. 315-336, 1990.
