n....:;¥:.i{: ::

 

 

 

 

 

 

 

 

 

 

 

Wm

.1. “1r .

 

 

 

 

 

 

 

1 1.? 2
int, .) a

4 72:...

m.
E mieuw

, , , .3 ..
2. 3E”... , . . a.

 

i. .. L» .53....

5.... . w Umfimflzvuw. ”Mafia

1!. 9.. é

Mews

||||||||||||lllllllllllIllllIllllllllllllUlllllllllll

3 1293 01055 957

This is to certify that the

thesis entitled

FINITE DIFFERENCE AND FINITE ELEMENT

MODELING OF CLOSED CELL CUSHIONS WITH

BLOCK AND RIBBED GEOMETRIES BASED ON
STRAIN DEPENDENT TERMS

presented by

ANDREW W . CHEN

has been accepted towards fulfillment
of the requirements for

MASTER degree in PACKAGING

 

 

 

1’
M or professor

Date OCTOBER 21; 1993

0-7639 MS U is an Affirmative Action/Equal Opportunity Institution

 

 

MBR‘ARY
Michigan State
University

 

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU Is An Atfinnetive Action/Equal Opportunity Institution
czwiddmedmpma-p.‘

 

FINITE DIFFERENCE AND FINITE ELEMENT
MODELING OF CLOSED CELL CU SHION S WITH
BLOCK AND RIBBED GEOMETRIES BASED ON
STRAIN DEPENDENT TERMS

By

Andrew W. Chen

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF PACKAGING SCIENCE
School of Packaging

1993

ABSTRACT

FINITE DIFFERENCE AND FINITE ELEMENT
MODELING OF CLOSED CELL CU SHION S WITH
BLOCK AND RIBBED GEOMETRIES BASED ON
STRAH\I DEPENDENT TERMS

By

Andrew W. Chen

Closed cell cushions behave in a manner such that idealized
models cannot be used. Hooke's Law is one such example. The
compression of a closed cell cushion is non—linear, visco-elastic and
time dependent. Because of this dependency, the model to replicate
cushion curve data must account for thermodynamics and heat
transfer. The derivation of this model is a based on three laws:
Newton‘s second law of motion, the Gas Law, and the First Law of
Thermodynamics.

Solving three differential equations derived from the three
laws based on finite difference yield results that are comparable to
existing cushion curve data. The problem is the derivation is tedious
and cannot be applied to non—block cushions. Hence, the polytropic
model is developed to enable the use of a more powerful numerical
method, the finite element method (FEM). The polytropic model is
put into a constitutive equation form where the variables are all
strain dependent. The result of the FEM compares well with

published cushion curve data and experimental results.

To Christina, Charles and Shang for all of their support

iii

 

ACKNOWLEDGMENTS

This thesis would not have been possible without the assistance of
Dr. Larry Segerlind, Dr. Ronald Averill and Eric Leaman who helped
in the development of the finite element program and Dr. Singh who
was always open for questions. Finally, much appreciation goes to
Dr. Gary Burgess. His patience and knowledge of packaging dynamics
and solid mechanics has lead to the successful completion of this
thesis.

iv

TABLE OF CONTENTS

LIST OF TABLES
LIST OF FIGURES -
CHAPTER 1 INTRODUCTION AND LITERATURE REVIEW
CHAPTER 2 MATERIALS AND METHODS

2.1 Derivation of the Governing Differential Equations

2.2 The Yield Pressure 6

2.3 The Total Cell Wall Surface Area S

2.4 The Convection Heat Transfer Coefficient H

2.5 The Total Air Mass m

2.6 Solving the Governing Differential Equations With

the Forward Finite Difference Method
2.7 The Polytropic Process
2.7.1 Method 1
2.7.2 Method 2

2.8 Force vs Compression Results

2.9 Finite Element Method

2.10 Finite Element Model

2.11 Finite Element Program and Input
CHAPTER 3 RESULTS

3.1 Simulation of Block and Ribbed Cushions

3.2 Block Cushion

3.3 Ribbed Cushion
CONCLUSION
APPENDIX A: FIN ITE ELEMENT PROGRAM
APPENDIX B: ARCEL 512 CUSHION CURVES
LIST OF REFERENCES

PAGE
vi
vii

12
14
14
14

18
27
28
29
32
34
35
4O
46
46
46
54
63
66
79
82

Table
Table
Table
Table

Table

Table

Table
Table

Table

Table

Table

Table

Table

Table

LIST OF TABLES

1 Gas Related Constants
2 Initial conditions
3 Solution to system of differential equations
4 Percent error between finite difference results
and Arcel 512 cushion curve'data for a
12 inch drop onto a 2 inch thick block cushion
5 Percent error between finite difference results
and Arcel 512 cushion curve data for a
36 inch drop onto a 2 inch thick block cushion

6 Variability in pxk at selected times from

Table 3 (k=1.072)
7 Input data for 2 inch and 4 inch thick block cushion
8 Force vs. compression results for the

10 in x 10 in x 2 in block thick cushion

9 Force vs. compression results for the 8 in x 8 in x 4 in
block thick cushion

10 Comparison of finite element deceleration results to
Arcel data

11 Geometry and initial modulus for the rib cushion

12 Force vs. compression results for a 4 inch thick
rib cushion with cross section dimensions as
specified in Table 11

13 Comparison between finite element results and the
equivalent volume results

14 Comparison between the finite difference,
finite element, and actual results at
extreme loading conditions

vi

PAGE
16
20
2 O

25

26

31
48

50

51

53
58

59

62

64

LIST OF FIGURES

Figure 1 Real cushion behavior compared to Hookean cushion

Figure
Figure
Figure

Figure
Figure

Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure

2 Schematic of weight in a free fall onto a cushion
3 Free body diagram of cushion and weight
4 Foam section with a hypothetical array
of square cells
5 Cushion curve for a 12 inch drop onto a
2 inch block cushion
6 Cushion curve for a 36 inch drop onto a 2
inch block cushion
7 The effect of k on the force vs. compression curve
8 Schematic of 1 dimensional linear element
9 General variables used for a block cushion
10 Flow chart of finite element solution
11 Force vs. compression curve for a cushion
12 Mesh variables for a block cushion
13 Three dimensional view of the rib cushion
14 Two dimensional view of the rib cushion
15 Finite element mesh for the rib cushion
16 Comparison between ribbed and block
cushions with the same overall dimensions

vii

PAGE
4
6
8

13

22

23
33
36
38
42
45
47
55
56
57

60

1

CHAPTER 1
INTRODUCTION AND LITERATURE REVIEW

Cellular polymers are multi-phase material systems that
consist of a polymer matrix and a fluid phase. The fluid is generally
a gas and is either trapped or continuous in the polymer matrix [9].
This has lead to the development of closed-cell and open-cell
cushions. Closed-cell cushions are foams in which the air is trapped
within individual unbroken cells. This type of foam is widely used in
the packaging industry for cushioning fragile items. Open cell foams
in which the air is able to circulate freely within the foam structure
behave differently because the air flows rather than compresses
[2,22]. In this study, only closed cell foams were studied.

The behavior of a closed cell cushion is generally described by
cushion curves [2]. Cushion curves relate the shock transmitted to a
product to the product weight, cushion thickness and bearing area of
the cushioning material. Because of the vast collection of cushion
curve information, the engineer and designer of protective packaging
can be overwhelmed. Naturally one would therefore try to reduce
the size of cushion curve data. A reasonable approach would be to
model the spectrum of cushion curve data in the form of descriptive
equations or some type of cushion model [3].

The straight forward approach toward modeling cushion
behavior would be to relate the cushion compression to the force the
cushion exerts upon a product. The solution of this solid mechanics
problem must consider: [5]

1. Equations of motion

2. Geometry

2

3. Material constitutive laws (stress-strain)
If the relationship between stress and strain is linear, the
constitutive law is known as Hooke's Law [15]. Hooke's Law is an
idealized model which some physical systems obey. The problem
with using Hooke's Law to describe the compression of a cushion,
however, is that cushion behavior is known to be non-linear, visco-
elastic, time dependent and strongly related to temperature. Despite
these dependencies, there are models that try to "fit" foam cushions
into the Hookean model [9,21,22].

The compression of foam cannot be accurately modeled by the
Hookean equations that assume a linear relationship between stress
and strain. To show this, it is sufficient to look at the prediction for
peak G-level to a mass dropped onto a cushion using Hooke's Law
and to compare it to published cushion data [2]. The basic steps are:

a. stress = modulus x strain (Hooke's Law)
O'=E8 eq 1
where o 2 stress (psi), 8 2 strain and
E 2 material modulus (psi)
b. energy absorbed by cushion =

(energy density) x (cushion volume)

0.2
U =—— Ab
2E( ) eq 2
where A is the cushion bearing area and b is the
thickness [26].
c. energy absorbed 2 potential energy [23]
U = Wh eq 3

where W is the weight dropped from a height h

3

d. Combining a, b, c and solving for the peak stress and

using o=F/A where F is force,

_ 2hWEA
F— ‘7.— eq 4

e. Force = mass x acceleration 2 weight x G
F = WG eq 5

f. Solving for peak deceleration,
2 E ‘
b(s) eq

where static loading is defined as s=W/A.

According to equation 6, G should continue to decrease as the
static loading increases. This is shown in Figure 1 under linear
cushion behavior. This agrees with known cushion performance, but
only up to a point. It fails to predict the increase in deceleration as
static loading increases. Since the only assumption involved in the
analysis is Hooke's Law, it must be cbncluded that equation 1 is in
error [22].

To model the cushion accurately, one needs to develop more
realistic constitutive equations. In order to develop the constitutive
equations one must first consider the nature of the compression

process. This'inevitably leads to considering heat transfer.

 

    
   

real cushion

Deceleration (G's)

linear cushion

 

 

static loading (psi)

Figure 1. Real cushion behavior compared to a Hookean cushion

As a closed—cell cushion compresses, the air volume decreases
and the pressure and temperature increase. Since air temperature is
now higher than the cell wall temperature, heat is transferred from
the air to the walls. This raises the wall temperature only slightly
because of its large heat capacity so that walls can be considered as a
nearly constant temperature heat sink. So during the entire
compression and expansion process, energy is extracted from air and
therefore not given back to the weight (see Figure (2)). As a
consequence, the stress depends not only on strain, but strain rate as
well [2]. Unfortunately, there has been very little work done on

models with strain rate dependent terms.

 

W = weight

 

 

 

   
 
 
 

A = bearing area
p = pressure
T = air temperature

h = drop height

x o = cushion thickness

 

rigid surface

Figure 2 Schematic of weight in a free fall onto a cushion

 

7

CHAPTER2
MATERIALS AND METHODS
2.1 Derivation of the Governing Differential Equations

In modeling a closed cell foam, three coupled differential
equations must be solved. The three equations are derived from
Newton's second law of motion, the Ideal Gas Law, and the First Law
of Thermodynamics [2]. Three equations are required since there are
three time dependent variables involved when the cushion is
compressed. The variables are the thickness of the cushion denoted
by x, the air pressure p within the cells, and the air temperature, T.
Figure 2 shows a weight in a free fall from a height h onto a block of
closed—cell cushion of thickness x0 with a bearing area A. Newton's
law can now be applied and the first governing equation can be
derived.

In order to solve for the deceleration of the weight, a force
balance is performed on the weight. The free body diagram in Figure
3 shows the weight and a cut away portion of the cushion. Newton's

second law states that:

2
zeld—g— eq 7
g dt
Wdzx
pA+0A—p0A_W:_g—-E§- eq 8

where o is the yield pressure required to compress the foam
structure (less the air), p0 is normal atmospheric pressure, p is the
cell air pressure during compression (assumed to the same for each

cell) and g is the acceleration due to gravitational pull.

 

lioliii
i.

      
 

V=AX

pA+oA

Figure 3 Free body diagram of cushion and weight

Solving for acceleration and rearranging equation 7. gives
dzx 0' + p —— po
zeal—S ‘1] eq 9
Equation 9 cannot be solved because the variable p is a function of
time.
The second governing equation which relates the air pressure
to the other variables is the ideal gas law, [18]
pV=nRT eq 10
where V is the volume of the air in the cushion, n is the fixed
number of moles of air trapped in the cushion, R is the universal gas
constant, and T is the absolute temperature. Given that the quantity

of air in a closed cell foam is constant, equation 10 states that pV/T

is constant at every instant. Since V=Ax, this leads to

_ z: _T_
P—Po x T. eq 11

Equation 11 relates the pressure to the position, but the new
variable T has appeared. Consequently, another relation must be
derived. The third equation required to solve the system of
equations is the energy equation. The energy equation is applied to
the fixed quantity of air within the foam using the first law of
thermodynamics applied to the closed system, [25]

work done + heat added = change in internal energy
During the compression and expansion of the cells surrounding the
air, there is work done on the air. This is known as simple boundary
work [25]. Heat is added to the air as a result of compression or

expansion and takes place through convection heat transfer to the

 

10

cell walls. The change in the internal energy of the air is the direct
result of a change in temperature [10,25]. Mathematically, the first
law of thermodynamics can be written as,
Sq + SW = du eq 12
For a closed system over a small time period, the work done, the heat
added and the change in internal energy can be expressed as,
' SW = -pdV eq 12a
Sq = -HS(T- T0)dt eq. 12b
where H is the convection heat transfer coefficient, S is the total
surface area of the cells in the foam matrix, T is the instantaneous air
temperature and To is the ambient air temperature. Also,
du = vadT eq 120
where m is the total mass of the air, and CV is the constant volume
specific heat for air. Substitution into equation 12 yields,
—pdV + -HS(T- To)dt = mCVdT eq 13
Rearranging with dV=Adx and dividing by dt, the third governing

equation becomes,

 

 

 

dT ( pA )dx ( HS )
—=— —— T—To
dz? mc. dt mc. ( ) “I 14
In summary, the three governing equations are:
dzx 0+ p— p0
Newton's Law dtz =8[_;———1] eq 9
x. T
The Ideal Gas Law I): ”(TEXT—j eq 11

 

 

dT pA )dx ( HS )
— = — — — T — To
Energy Balance dt [va dt m C. ( ) eq 14

l 1
The three coupled equations may be reduced to two equations
since equation 11 is not a differential equation. Substituting p in

equation 11 into equations 9 and 14 yields

 

 

 

dzx 0—pa ) (gpaxo)T

dt2 g( s STo x eq 15
dT A ano T dx HS HS
— = — — —— + T. — T 16
dt MCv To x dt MCV "ICV eq

Since the terms in parenthesis are constants for a particular cushion,

these two equations can be written as:

2
d—x = K1+ 10(5)
X

 

 

 

 

(1,2 eq 15a
a7" de
—=—K3 —— K4—K5T
dt (x dt)+ ( ) eq 168
where
0— 0
K1: [—p‘—:l eq 17a
s
K2- (gpoxa)
7 ST. eq 17b
+ A ana
K3: va To eq 17c
(11sz
K4: mc, " eq 17d
(”5)
K5: va eq17e

In order to examine the effect of cushion geometry and cell
dimensions on the constants K1 through K5, certain terms in the
above expressions need to be examined. Specifically, they are 6, s, m,

and H.

12

2.2 The Yield Pressure 6
' The yield pressure, 0, is the stress required to compress the
cell matrix alone and contributes significantly to decelerating the
product. The polymer matrix can be considered a thin wall
membrane which supports no bending moments [20,21,22]. The cell
can be assumed to under go both buckling and yielding [9]. Modeling
of the cell matrix is difficult due to the variability in cell size and the
complexity of geometry. Many geometry's have been proposed such
as the spherical cell connected by a cylindrical strut model, the
pentagonal dodecahedron cell and the tetrakidecahedronal cell [9,22].
The geometry will vary depending on the material and the foaming
process. Matonis [14] has proposed an equivalent model to simplify
the analysis of the cell structure. The model proposed is a cubical
plate model which consists of layered cubical arrays as in Figure 4.
Figure 4 shows a hypothetical section from a closed cell
cushion. The yield pressure 0, is defined as the vertical surface
pressure which causes the cell walls to yield. A simple force balance
yields,
ozoyp(Aw/A) eq 18
where oyp is the material yield point stress, Aw is the total cross
sectional area of only the walls and A is the apparent bearing area
(1*W). The total amount of wall support area in Figure 4 is,
Aw = [1(w/a) + w(l/a)]t = 2(lwt)/a eq 19
Rearranging and substituting results in

G = (2t/a)0'yp CC] 20

13

cell
+|H A k AJ—V

 

 

 

 

 

 

 

 

 

 

 

Figure 4 Foam section with a hypothetical array of square cells

14
Typical values for the cell dimensions for Arcel 512 [16] are, a=45
mils, t=0.3 mils [21] and Gyp=1500 psi [9,20]. Applying equation 20
to Arcel 512 gives a yield pressure of approximately 20 psi.
2.3 The Total Cell Wall Surface Area S
The total cell wall surface area can be calculated by considering
the average surface area of one cell and multiplying by the number

of cells in a volume of cushioning material.

52(602) Lfix—o :6lwx0 =65 eq 21
aa a a a

 

where V0 is the original volume of the cushioning material.
2.4 The Convection Heat Transfer Coefficient H

The heat transfer rate is related to the overall temperature
difference between the gas in the encapsulated cell and the cell wall
by the convection heat transfer coefficient, H [9]. This term has been
found to be 4 BTU/(hr-ft2-0F) for cells with the dimensions given in
2.2 [2].

2.5 The Total Air Mass m

The mass of the air is estimated to be the volume of the
cushion times the density of the air. The foam itself contributes very
little to the volume of the cushion. This can be realized by taking an
open cell cushion and compressing the cushion. The open cell
cushion can be made to flatten out almost completely. The mass is,

m=p(l*w*x0) eq 22

where the air density is p=.0000428 lb/in3 =.074 lb/ft3 [17].
Substituting the results for o, s, m, and H into the expressions for K1

through K5 gives,

 

 

 

1 20 ‘-t
K1: — ’P — o—
gs[ a p s] eq 23
K2- gpoxo 1
‘- To S eq 24
p0
K =
3 (mm) eq 25
6HT. 1
K4— —— ——
( Cvp )(a) 6“ 26
6H 1
K5: —
(010an 6“ 27

The terms on the right hand sides of equations 23 through 27 have
been grouped using parenthesis in order to illustrate the effects of
material properties, cell dimensions, cushion dimensions and
environmental terms on the constants. Without this simplification,
one cannot deduce the effects of modifying the cell or cushion
structure. The constant K3 for example is dependent on air _
properties alone where as K4 and K5 depend on cell size as well. The
constants K1 and K2 depend on air properties, cell dimensions,
cushion dimensions, and the weight of the falling mass.

Since K1-K5 depends on certain properties of air, substituting
the physical constants in Table 1 into the expressions for Kl—KS

further simplifies them.

16

'_1"_able 1. Gas Related Constants
H = 4 Btu/hr/ft2/0F (heat transfer coefficient)
p0 = 14.7 psi (initial air pressure of cells)
To = 530 OR (initial cell air temperature)
CV = .1715 Btu/lboF (constant volume specific heat of air)
p = .0000428 lb/in3 = .074 lb/ft3 (density of air at 530 0R)

17

Substituting g=386.4 in/sec2 and the values in Table 1 into
equations 23 to 27 gives the constants in terms of quantities which

relate to cell sizes and cushion size:

K1=386.4é—)—[2-(0—’;1]—14:7-(s):l eq 28
K2 = 1072696) eq 29
K3 =.04 e q 30
K4 = 335060661) eq 31
K 5 = 6322(%) eq 32

where the static loading, 3 and yield pressure, oyp are in psi and the
cell dimensions, a and t, are in mils.
Now for a particular cushion like Arcel 512 [16] where a = 45

mils, t=0.3 mil, oyp =1500 psi, the expressions are further simplified

 

 

 

to
2048 in
K1=( S —386.4)[Sccz] eq 33
1 lb2
K2=10.72o — ——-—-
(x°)(s)[secz-°R-in2] eq 34
K3=.04 eq 35
”R
K4=74458[scc:| eq 36
. 1
K5=140[Q] eq 37

18

At this. point, it is evident that only K1 and K2 depend on the
types of variables normally associated with cushion design: the static
loading s, and the cushion thickness x0. An additional parameter not
contained in K1 and K2 but important to cushion design, is the drop
height information. The drop height parameter will enter the
analysis through the governing equations as an initial condition on

velocity. Specifically, the initial or impact velocity is [4],
vi =s/2og'h eq 38

Once the weight, drop height, and the original cushion thickness are
known, the constants K1-K5 are known and the differential equations
can be solved numerically for the pressure, p, temperature, T, and
thickness, x. The procedure for doing this is outlined in the next
section. A computer program based on this procedure was written to
accept drop data (W, w, 1, x0, and h) and generate the shock pulse.

2.6 Solving the Governing Differential Equations With the Forward

Finite Difference Method
From the three governing equations, a series of values for

thickness, temperature and pressure can be generated. In order to
generate these values, the forward finite difference method was
used. The technique involves initializing the variables and stepping
the equations forward in time with an increment of At. This process
will continue until the experimenter decides to stop the march
forward. Referring back to the summary of the differential
equations, notice that only pressure can be found from the three
equations (equations 9, 11, 14). The thickness and temperature

variables are in differential form. In order to calculate the thickness

 

19

and temperature, equations that relate position to acceleration and
temperature to temperature rate are required.

The equations that relate position and velocity to acceleration
are equations derived from kinematics in one dimension [12,13,23].
These equations are only valid for uniformly accelerated motion. For
some differential time dt, these equations hold true and for some

finite time At, they are only approximate.

 

_ dx ldzx 2
x—XO+ZAt+§dtozAt Cq 39
£—_dl+d_2x.At
dt dz, cit,2 eq 40

In a similar manner, the new temperature will depend on the rate of

temperature change.

T = To +£At
dt

0

eq 41

The three differential equations (9,11,14) coupled with the three
kinematic equations (39,40,41) are the six equations necessary to
solve the finite difference problem. The equations were solved'with
the Macintosh version of the Excel 4.0 spread sheet. An example of
the calculation with the spread sheet is shown in Table 3. The drop
data, material constants, environmental constants and gas related

constants used for this example are shown in Table 2.

 

20

386.

4

0.
1500

2
0.001
-166.

K1
K2
K3
K4
K5

Table 2 Initial conditions

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

At T (R) X (in) P (psi) dX/dt G dT/dt d2x/dt2

0 0 530.00 2.000 14.700 -166.80 19.00 17871 7342
1 0.001 547.87 1.837 16.545 -159.45 20.85 16718 8055
2 0.002 564.59 1.681 18.626 -151.40 22.93 15694 8859
3 0.003 580.28 1.534 20.977 -142.54 25.28 14730 9767
4 0.004 595.01 1.397 23.630 ~132.77 27.93 13734 10792
5 0.005 608.75 1.269 26.601 -121.98 30.90 12587 11940
6 0.006 621.33 1.153 29.882 -110.04 34.18 11135 13208
7 0.007 632.47 1.050 33.414 -96.83 37.71 9187 14573
8 0.008 641.66 0.960 37.060 -82.26 41.36 6534 15981
9 0.009 648.19 0.886 40.575 -66.28 44.87 2997 17340
10 0.01 651.19 0.829 43.596 -48.94 47.90 -1474 18507
11 0.011 649.71 0.789 45.686 ~30.43 49.99 -6684 19315
12 0.012 643.03 0.768 46.439 -11.12 50.74 -12116 19606
13 0.013 630.91 0.767 45.642 8.49 49.94 -17000 19298
14 0.014 613.91 0.785 43.386 27.78 47.69 -20575 18426
15 0.015 593.34 0.822 40.045 46.21 44.34 -22386 17135
16 0.016 570.95 0.877 36.126 63.35 40.43 ~22433 15621
17 0.017 548.52 0.948 32.101 78.97 36.40 -21078 14065
18 0.018 527.44 1.034 28.300 93.03 32.60 -18830 12597
19 0.019 508.61 1.133 24.898 105.63 29.20 -16164 11282
20 0.02 492.45 1.244 21.951 116.91 26.25 -13429 10143
21 0.021 479.02 1.366 19.446 127.05 23.75 -10846 9176
22 0.022 468.17 1.498 17.336 136.23 21.64 -8527 8360
23 0.023 459.65 1.638 15.562 144.59 19.86 -6516 7674
24 0.024 453.13 1.787 14.067 152.26 18.37 ~4812 7097
25 0.025 448.32 1.943 12.801 159.36 17.10 -3394 6608

Table 3 Solution to system of differential equations

 

 

 

21

The procedure used in the Excel spread sheet involved first
initializing all the variables at time t=0. Each column and row in the
worksheet represents one of the variables defined earlier at the
instant shown. For example, d2x/dt2 represents the acceleration of
the mass on the cushion and dx/dt is the velocity of the mass at a
instant in time. dT/dt is the temperature rate, T is air temperature, x
is the thickness of the cushion and p is the air pressure. G is the
dimensionless unit of acceleration defined as the acceleration at a
particular time interval divided by gravitational acceleration. This
variable is needed to specify the severity of the shock. When the
pressure p returns to its initial value of 14.7 psi, the weight is
leaving the cushion and the finite difference calculation is completed.

For the particular run in Table 4 the peak deceleration is seen
to be 50.74 G's‘and the shock duration is approximately 24 msec.
The rebound velocity is 152.26 in/sec which implies the coefficient
of restitution [4] to be e=.91. Considering the complexity of the three
governing differential equations, the solution over the duration of
the shock pulse seem to produce reasonable results. For example the
velocities and accelerations have magnitudes and direction that are
reasonable. The value of G increases and decreases and the
coefficient of restitution is less than 1. The actual value of G from
the Arcel Cushion Curves in Appendix B is 50 G's.

To check the validity of the model and the finite difference
solution, the results of the finite difference solution for 12 inch and a
36 inch drops were compared to the published cushion curves [16].
Figure 5 and Figure 6 show that the finite difference equations

predict the concavity of the curve and the results are similar to that

22

 

 

1 0 0 -

8 0 -
Q
a 6 0 ~ —¢— Arcel Data
0
E '“‘""0”"°‘°' Finite Difference
2 4 0 -
CD
0
8

2 0 -

O I l I I l I

 

static loading (psi)

Figure 5 Cushion curve for a 12 inch drop onto a 2 inch block cushion

 

23

500

400

—<>— Arcel Data
300-

------o-----'-- Finite Difference

 

Deceleration (G's)

200-

100-

 

 

 

static loading (psi)

Figure 6 Cushion curve for a 36 inch drop onto a 2 inch block cushion

 

24
of Arcel 512 (see Tables 4 and 5 for errors). Assuming that the finite
difference method solves the system of differential equations
accurately, the model seems to under estimate G's for a low drop
height and over estimate G's for a high drop height. This is partially
due to the fact that actual experimental cushion curves are in error
by 10% to 15% [27] on average and partially due to the use of
constant values for air properties. The heat transfer coefficient H in
particular is only accurate to within i50% and is known to be a
function of the temperature difference [10, 27]. The constant value
of H24 Btu/hr/ft2 used in the solution could therefore be in error by
as much as i3 Btu/hr/ftz. The experimental error in reading the
cushion curve must also be included in accessing the accuracy of the
finite difference method. At this point the important fact to note is
the ability of the model to simulate cushion performance and the
significance of incorporating heat transfer and thermodynamic terms.
The motivation now will be to simplify the complexity of the three
differential equations by modeling the cushion by a polytropic

pI'OCCSS .

25

 

Table 4 Percent error between finite difference result and Arcel 512
cushion curve data for a 12 inch drop onto a 2 inch thick
block cushion

26

 

.2
.4
.4
.3
.8
.0
.5
.6

Table 5 Percent error between finite difference results and Arcel 512
cushion curve data for a 36 inch drop onto a 2 inch thick
block cushion

27

2.7 The Polytropic Process

The solution process so far has been tedious and complicated
because the generation of G's requires lengthy calculations. An
alternative method to generate G's is to model the compression of air
as a polytropic process [25].

ka = constant eq 42
This model explicitly excludes the heat transfer and the
thermodynamic terms. Implicitly, all of these terms are included in
the polytropic constant k [25]. The motivation for using the
polytropic model is to simplify the stress-strain analysis so that
cushion geometries which are more complicated than a simple block
may be analyzed. This then allows for techniques such as the finite
element method to be used. It also simplifies the model considerably
by requiring only one physical constant, k instead of the many
thermodynamic terms previously used.

As a result of the finite difference calculations in Table 3,
values of position(thickness) and pressure at specific times are
known. When a complete table has been generated, the material
constant k in equation 42 can be found by fitting the known data to
the polytropic equation. The value of k will determine the amount of
heat transfer that occurs in the cushion. The value of k is 1.4 for an
adiabatic process: in an adiabatic process, compression is assumed to
occur so quickly that no heat is transferred to the cell walls. The
value of k is 1.0 for an isothermal process: in an isothermal process,
compression is assumed to occur so slowly relative to the heat
transfer rate. that the air remains at constant temperature as it is

compressed.

28

An adiabatic process would result if H or S were extremely
small or if the cushion is compressed quickly as in a drop. If H or S
were very large or if the cushion is compressed very slowly, then the
process would approach isothermal conditions. Since H and S are
both finite in reality, the compression process will be neither
isothermal nor adiabatic. The expected result is 1.0<k<1.4.

Predicting the polytropic constant involves fitting equation 42
to the data in Table 2. Since kazconstant and the bearing area is
constant, we can fit the data to pxkzconstant instead. Since the

pressure and temperature are initially p0 and x0, this leads to:

k .
P0 x0

The objective now is to assume that equation 43 is true and to choose
k to "best fit" the data in Table 3. This is a necessary step in order to
develop the constitutive equations for a closed cell foam.
2.7.1 mired—1

Because the air compression process is not truly a polytropic
process, equation 43 will not be satisfied exactly because k is chosen
to best fit all the data points for pressure and position
simultaneously [6]. As a result, there will be an error between the
actual value and the expected value. Defining the error at any

instant i using the numbers from the i th row in Table 1 as

k
P0 x0

we can find the variance in the error. This done by calculating the

sum of squares of the errors,

29

N o 0 k 2
SSE=Z[(-&Ii) —1] eq 45
i=1 p0 x0

The best fit involves choosing k so that the SSE is a minimum [6].
From the theory of differential calculus, minimizing equation 45
would be very difficult as it would result in a highly nonlinear
equation to be solved for k. Instead a similar method will be used.
2.7.2 _M_ethod_2
Instead of using equation 43 in its present form, take the

natural log of both sides,

k
P x _
“(Elmizl ‘0 eq 46
The error can now be defined as
' k
x
8i =ln(—pI-:—]+kln[;) eq 47
and the SSE becomes
N N . . k 2
SSE=28i2=2 ln(-&)+kln[i) eq 48
i=1 i=1 p0 x0

Minimizing equation 48 is now much simpler since k is no longer an
exponent. Taking the first derivative of equation 48 and setting the
result equal to 0 gives,

d N . . .
&(SSE) = 22[1n£'— + lei] 1n(i] = 0 eq 49

i=1 0 x0 x0

N N 2
Pi xi xi
ln—ln— +k ln— =0
EA x0] 1:21[ x0] eq 50

i=1 p0

 

_. P x
k— le 0x 20 eq 51
tit—i]
21 x.

The result for k using equation 51 and the data in Table 3 is k=1.072.
Using this k, the values of pxk for several points in Table 3 were
generated. The results are shown in Table 6. Even though pxk is not
constant, it appears to vary by no more than 9%. Next, a series of
tables were generated similar to Table 3 for five different drop
heights and the same initial conditions. This was done to examine
the effect of drop parameters on k. Using equation 51, the results
ranged from k=1.072 for a 6 inch drop height to k=1.188 for a 36
inch drop height. These results show that cushion compression is
closer to being isothermal than adiabatic and that k is relatively

insensitive to drop conditions.

31

 

 

 

 

 

 

 

 

 

 

time (ms) pxk
0 33.61
10 34.81
20 28.20
3 0 31.13
average: 31.94
standard deviation: 2.93 (9%)

 

Table 6 Variability in pxk at selected times from Table 3 (k=1.072)

 

32

2.8 Force vs. Compression Results
Figure 7 shows the force vs. compression results for a 2 inch
thick block cushion. These results were obtained using the
polytropic equation with the following definitions of force and
compression.

1. Force = (p - po)A
2. Compression: u = x0 — x

Substituting the expression for p from equation 43, the force vs.

compression relation for a polytropic process becomes,

k
x
F: A —°- —1
3. p0 [[xo—u) :|

The force vs. compression curves in Figure 7 were obtained using a
polytropic constant of k=1.072 and k=1.1882 along with the values
A2100 in2 and x022 in. The results show that the force vs.
compression curve is fairly insensitive to the variability in k.
Therefore, the remaining tests performed in this thesis will be based
on k=1.072 for Arcel 512. The important part to note here is that the
constants in the three differential equations can now be lumped into
one polytropic constant and the equations of elasticity may be

developed from this.

 

 

33

25000-

 

20000-

force (lb)

 

 

 

compression (inch)

Figure 7 The effect of k on the force vs. compression curve

 

 

34

2.9 Finite Element Method

The finite element method is a powerful tool that models a
complex domain as a collection of geometrically simple sub domains.
The finite element method can be interpreted as a piecewise
application of the variational method to solve differential equations
in which the approximation functions are algebraic polynomials and
the undetermined parameters represent the values of the solution at
a finite number of pre-selected points called nodes [17]. The
variational method is a technique used to solve a differential
equation approximately by putting the differential equation into a
variational form and taking the approximate solution to be a
combination of shape functions.

The accuracy of the finite element method depends upon the
element used to model the complex domain. Elements can have
nodes that have one, two or three degrees of freedom. The number
of nodes in the element as well as the number of elements in the
domain all contribute to the accuracy of the approximation. The type
of element and mesh used will generally require some knowledge of
the physical system to be modeled.

When one looks at a typical cushion, especially a block cushion
with a square bearing area, lines of symmetry become noticeable. A
three dimensional problem can be reduced to a one dimensional
problem and expect a reasonable degree of accuracy because the
depth dimension is much greater than the side dimensions. The
objective is to show that an impact on a closed cell cushion can be

accurately modeled by the polytropic equation and analyzed

35

accurately by the finite element method. The results will be
compared to existing experimental data.
2.10 Finite Element Model:

Since the material constant k has already been determined and
the problem has been reduced to a one dimensional problem, the
closed cell cushion may now be viewed as an elastic solid whose
modulus of elasticity [15] is a function of the strain. The governing

equation becomes [15,26]

 

d
P(x)=(o+E<x>E“)A(x> eq 52
Ol‘
dP dA( )
E70 dxx+ +dx(A( )E(x )—) “1 53

where x is the coordinate variable along the thickness, E(x)=material
modulus at x, u=the displacement of a node during compression,
P(x)=the compressive force, o=yield pressure and A(x)=cross sectional
area at x (see Figure 8). The finite element process is able to solve
the governing differential equation (eq 53) for a variable material
modulus and a variable cross sectional area. The way in which the
area varies with depth is determined \by the geometry of the cushion.
The way in which the material modulus varies with depth is not
obvious and therefore must be determined by performing a force
balance on an element of cushion and using the polytropic
compression model. During compression, the stress applied to the
cushion must overcome the buildup in air pressure over and above
atmospheric pressure,

a: —
p p, eq 54

36

node 1

 

element dx
e(X)

 

 

 

node 2

e ———>o—<>-—e
———>

Figure 8 Schematic of 1 dimensional linear element

 

37

Using the polytropic relation for p gives,

I:

AM.

0': Vk

-Po eq 55

For a particular element, V0=Ato and V=At where A is the cross
sectional area of the element, V0 is the volume of the uncompressed
element, V is the volume of the compressed element, to is the
original thickness before the stress 0 is applied, and t is the thickness

after. Substituting into equation 55 and solving for 0 gives

k
__ 5.. _1
0_p0 I . eq 56

Now since t=to-u from Figure 9, where u is the compression,

1‘ t l 1
—"—=——-Q—= =———— eq 57
t to—u 1__l:l_ 1—8

to
where e = u/to is the strain. Substituting yields,

1 k
0=p0[1—_;) —1 6q 58

Equation 58 is the stress vs. strain relation for the cushion element.

 

Dividing both sides by the strain and noting that by definition, the

modulus is the ratio of stress to strain [15] gives

 

1 p
E: —1—°—
[(l-e)" ]8 eq 59

In the small strain limit, E approaches kpo=(l.072)(14.7)=15.76 psi
and steadily increases for larger strains. Both the material modulus
and cross sectional area are now a function of x. The area is an
explicit function of x because. it varies directly with position as the
geometry of the cushion dictates. The modulus is an implicit function
of x because it varies with the strain and the strain varies with

position.

38

 

bearing area, A

 

it.

 

 

 

T...

Figure 9 General variables used for a block cushion

 

39

Because the modulus is a function of strain and the strain
distribution is not known before hand, the material modulus must be
initially guessed. A reasonable guess would be the value, E=14.7 psi
for small strains. The problem can then be solved with each element
having this modulus. The resulting strain field can then be used to
correct the modulus for each element using equation 59 and the
process is repeated. The program will iterate until the material
modulus converges upon a value for each element. Upon
convergence, the results will eventually be used to calculate the peak
G delivered to a product in a drop.

Now that the behavior of the cushion during compression has
been modeled as a purely mechanical process, it is possible to
examine cushions with more complicated geometric shapes like
ribbed cushions. Ribbed cushions are merely block cushions with
material removed at the base which gives them the appearance of a
thinner block cushion standing on legs. They have the advantage of
limiting the required amount of foam for protection while
maintaining full support of the base of the product. The method of
analysis which relies on the solution to the two coupled differential
equations has the disadvantage of being able to handle only block
shaped cushions (constant cross sectional area). The method also
requires the knowledge of several material properties and cell
dimensions. The advantage in modeling the process as a polytropic
one is that all of these properties are lumped into one k, which is
hopefully fairly constant for a particular brand of cushion. It is now
possible to use broad based numerical methods such as the finite

element method to analyze odd shaped cushions.

40

2.11 Finite Element Program and Input:

The method to solve the system of equations involves the
implementation of a computer algorithm in a Fortran computing
environment. A computer code developed with assistance from R.
Averill [1] will be used to solve equation 53 [11,17,19]. The code was
written to model each element in a closed cell cushion according to
equation 59. The equation for each element is then combined into
system of equations and solved numerically by matrix algebra. The
algorithm to solve for the stress and strain consists of three basic
parts: the preprocessor, processor and the post processor [17].

In the preprocessor, the boundary conditions, mesh geometry
and element information is read in. This includes the data regarding
how many nodes per element, the number of elements, how the
nodes are restricted (boundary conditions), the type of interpolation
functions, and the number of Gauss points to numerically integrate.

The processor is where the local stiffness matrix is generated
for each element. The local stiffness matrix uses the interpolation
functions to generate the solution to the variational approximation
for each particular element. The local stiffness matrices are then
assembled into a global stiffness matrix where the influence of
adjacent elements are accounted for [17,19]. The contribution to a
stiffness coefficient at the global level could be the result of several
local element stiffness coefficients. The global stiffness matrix is
generally denoted as [Kij]. Once the global stiffness matrix is
assembled the system of equations are set up where the stiffness
matrix multiplies the displacement vector u and is set equal to the

nodal force vectors P: [K][u] = [P].

41

The next step in the finite element method involves imposing
the boundary conditions on system of equations. The boundary
conditions include specifying displacements or forces at a node. Once
the boundary conditions are imposed, the system of equations are
solved. The solution of interest that the computer program returns
are the displacements for each node in the mesh.

The displacements are then used in the post processor to obtain
other pertinent information. In the post processor, the gradients
(strains) and stresses are calculated for each element. From this
data, the material modulus can be calculated. The preprocessor,
processor and post processor are located in the main program.
Routines that solve the differential equation are located outside of
the main program as subroutines. The details of the complete
program are shown in the appendix and a summary of the algorithm
and the routine is shown is Figure 10.

The first subroutine which is called from the main program,
MAIN, is STIF which develops the local stiffness matrix. While
generating the local stiffness matrix for each element, STIF will call
the subroutine SHAPE to develop the interpolation function, which
will be used to generate the solutions to the variational equation of
that element. The SHAPE routine contains information related to the
linear interpolation functions. The linear interpolation functions are
used in the STIF routine to solve equation 53 in variational form. As
the program treats each element in the mesh, the local stiffness

matrix is assembled into a global stiffness matrix by the subroutine

 

 

42

 

 

 

 

 

 

 

 

 

 

 

 

SUBROUTINTE SHAPE
MAIN PROGRAM
_ _ _ __ Evalutes element
'— l interpolation functions
Preprocessor at Gauss points
*Reads and
generates data I
*Initialze arrays SUBROUTINE STIF
I *Initialize variables I

Generates Element
matricies for:
I I *second order equations

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Processor '
I *Generates I
element (STIF) SUBROUTINE ASSEMBLE
matrrcres
*
I 1:883:32: I V Assembles local matrix
*1 q into banded global matrix
I 8100833 d' - I
oun ary con itions
(BNDY) SUBROUTINE BNDI
$38,531:? uations ' Imposes specified
I l values on displacement
values of the problem

 

I SUBROUTINE SOLVE

I Postprocessor

*Calculates the

I gradients and stress
*Strajns are weighted I

 

Solves the equation
by Gauss elimination

 

 

. . technique(for
and substituted into banded symmetric
modulus equation .
I I equations)

*Compare calculated
modulus to modulus

 

 

I presently stored in I
array
*If the modulus
converges than stop I
I else return to
beginning for new

I iteration I

Figure 10 Flow chart of finite element solution [17]

 

 

 

 

43
ASSEMBLE. ASSEMBLE puts the global matrix in a banded matrix
form. The subroutine BNDY then imposes the specified boundary
conditions to the displacement and force vectors The subroutine
SOLVE uses the Gauss elimination technique to solve the matrix for
the unknown displacements.

The solution to a problem is solved by iteration. This is
necessary because the material modulus depends on the strains and
the strains are the solution to the problem. The initial condition
imposed on equation 53 to begin the calculation process involves
prescribing a compressive displacement of the upper surface while
fixing the base and selecting an initial material modulus. The
program then solves for the nodal displacements. The first iteration

is completed when the displacement at each node is calculated. From

 

this solution, values of stress and strain are calculated. Since the
strains determine the modulus through equation 59 and the modulus
used to solve for these strains in the first place was initially guessed,
a check is required to verify the solution. A comparison is made
between the old modulus (guessed) and the new one based on the
calculated strain. The second iteration now uses the new modulus
and the process continues until the modulus converges for each
element. In order to stabilize the iteration so that the solution
doesn't oscillate or diverge, a weighted strain average may be used
to calculate the new modulus for each element. The weighted strain
average using 10% of the new strain plus 90% of the old strain was

used for each element. The process is shown in Figure 10.

 

44

2.12 Force vs. Compression Transformation to G's

Once the iteration process converges for a prescribed boundary
displacement, the only result of interest is the corresponding force on
the upper surface required to produce this displacement. In order to
be able to use this information to predict the shock to a product
dropped onto a cushion, the whole iterative solution process must be
performed again for various displacements. By prescribing larger
and larger displacements and solving for the corresponding forces, a
force vs. compression relationship, F vs. u, for this cushion like the
one shown in Figure 11 can be developed. The use of this force vs.
compression curve to predict G level is now straight forward. The
potential energy Wh to .be absorbed by the cushion in dropping the
weight W from a height h is,

d».
Wh= Jqu eq 60
0

where dm is the peak dynamic deflection of the cushion. In other
words, the dynamic compression corresponds to the point on the
curve where the area underneath is equal to Wh. The corresponding
force F is related to the shock through Newton's Law

G = F/W eq 61

where F is the force corresponding to the compression dm.

 

 

force (lbs)

 

45

 

25000 -
20000 '-
15000 -
10000 - "
“—F— — — -Wh — — —
5000 - y -
. dm
0 I I I I
0 05 l 15

Figure 11 Force vs.

compression (inches)

compression curve for a cushion

46

CHAPTER3
RESULTS
3.1 Simulation of Block and Ribbed Cushions

The polytropic model and finite element program will be tested
for two cushion types: a block cushion and a ribbed cushion. The
purpose of the test is to check the theory that the compression of a
closed cell foam may be modeled as a polytropic process and to solve
the problem using the finite element method. The stabilization and
convergence of the iterative solution generally takes from two
iterations for block cushions to about thirty iterations for the ribbed
cushions at large strains.

The block cushion solution will be checked quantitatively
against the Arcel 512 cushion curve data and the ribbed cushion
against a block cushion with the same overall length, width and
thickness dimensions. The ribbed cushion solution will be further
tested by comparing the results to the equivalent volume method [8].

3.2 We

The finite element setup for analyzing the block cushion is
shown in Figure 12. Four elements with two nodes per element were
used as shown. The cross sectional areas were all taken to be A(n) =
(l)(w) and the moduli B were calculated in accordance with equation
59 using k=1.072. The force vs. compression curves for the two
different block sizes were calculated using the method outlined in
Figure 10. The first block was (l)(w)(t) -_- (10)(10)(2) in3 and the
second was (8)(8)(4) in3. The details of the input data to the finite

element program are shown in Table 7.

 

4
1 A(4).E(4) . 5

47

Block Cushion
+x Configuration A

1 / .

 

 

 

 

. n
'l
2
A(2>.E<2) I3
3w
A(3),E(3)

 

f\
\r

 

 

 

 

 

 

C: I
4

Figure 12 Mesh variables for a block cushion

 

48

 

 

LInput data FEM notation “2 in. thk. block “4 in. thk. block ll

 

 

 

 

 

 

 

 

 

 

 

 

 

Nodes/element nodelm 2

Degrees of ndfnod 1 1

freedom per

element

number of Gauss ngpstif 2 2

pts to calculate

the displacement

number of Gauss ngpgrd 2 2

pts to calculate

the strain

Number of nelmsh 4 4

elements in mesh

domain domain 2 inch 4 inch

(thickness)

material constant k 1.072 1.072

material modulus (initial guess) E(x) 14,7 psi 14.7 psi

cross sectional area A(x) 100 square 64 square
inches inches

 

 

Table 7 Input data for 2 inch and 4 inch thick block cushion

 

49

The results of the finite element solution for the force vs.
compression curve for the 2 inch block are shown in Table 8. The
displacement at the upper surface was gradually increased in steps
of .2 inches and the corresponding force was determined by the
program. Since the prediction of G's requires areas under the force
vs. compression curve, a seventh order polynomial was passed
through the data points and then integrated to obtain the area as a
function of compression. The polynomial which fits the data in Table
8 is:

F = 31310.9x7 - 194322.5x6 + 490567.5x5 — 641655.5x4 +
460788.6x3-175799.1x2 + 32669.0x + 34.2

The results for the 4 inch block cushion are shown in Table 9.
The polynomial which fits the data points is:

F = 212.5x7 — 2622.7x6 + 13085.9x5 - 33568.7x4 +
46723.3x3-33878.9x2 + 11471.3x + 147.2

 

 

50

 

 

 

 

 

 

 

 

 

 

 

 

 

prescribed corresponding energy
compression force absorbed
(in) (lb) (in*lb)
0 2000.0 0
0.1 2083.1 118.48
0.2 2175.9 339.56
0.4 2397.6 803.74
0.6 2685.2 1297.22
0.8 3072.71 1879.72
1 - 3623 2549.85
1.2 4458.2 3342.80I
1.4 5878.4 4365.04
1.6 8792.7i 5792.98
1.8 17909.8 8266.68

 

 

 

 

Table 8 Force vs. compression results for the 10 in x 10 in x 2 in
block thick cushion

51

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

prescribed corresponding energy
compression force absorbed
(in) (lb) (in*lb)
0 1280.0 0
0.1 1307.51 61.89
0.2 1337.0 185.21
0.4 1399.0 492.47
0.6 1469.4 797.39
0.8 1550.1 1088.74
1 1641.0 1389.91
1.2 1747.2 1721.35
1.4 1870.1 2089.47
1.6 2014.7 2490.17
1.8 2188.2 2917.79
2 2398.1 3373.39
2.2 2658.6 3868.45
2.4 2988.8 4423.37
2.6 3420.2 5062.48
2.8 4006.4 5810.35
3 4845. 6696.56
3.2 6138.24 7779.22
3.4 8365.4 9200.00I
3.6 13030.4 11286.25

 

block thick cushion

Table 9 Force vs. compression results for the 8 in x 8 in x 4 in

52

Note that in both Tables 8 and 9 the force is not zero when the
compression is zero. This is due to the yield pressure, (5 = 20 psi,
required to deform the cell structure (less the air). In fact the force
at zero compression is easily calculated using force 2 yield pressure x
area.

Now that the force vs. compression curves for the two block
cushions analyzed are known, they may be used to find the shock
delivered to a product in a drop. An example showing how this is

done is given below.

Example: A 32 lb product is dropped from a height of 12 inches onto
the 8" x 8" x 4" cushion.

Solution: The product of weight and drop height is,

Wh=(32)(12) = 384 in-lb. Since this is the energy which must be
absorbed by the cushion and since energy absorbed is just the area
under the force vs compression curve, the corresponding force from
Table 9 is found to be F = 1377.3 lb. From Newton's Law,
(force=(mass)(acceleration)), the corresponding peak deceleration is,
G = F/w = 1377.3/32 = 43.

The result is then compared to Arcel 512 cushion curves where the
G's for a static loading of .5 psi and cushion 4 inches thick is,

G=40

This same procedure, for an 8" x 8" x 4" cushion, was applied to other
drop situations in order to test the method over a broad range of

dynamic deflections. Table 10 shows the results.

 

53

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

drop weight (lb) thickness bearing predicted predicted actual

height (in) (in) area (1112) dm (in) G's" (FEM) G's
(ARCEL)

12 32 4 64 .10 43 40

12 160 4 64 1.25 11 10

30 32 4 64 .25 47 45

30 160 4 64 2.60 20 20

48 32 4 64 1.12 53 50

48 160 4 64 3.10 38 39

Table 10 Comparison of finite element deceleration results to Arcel

512 data

 

54

The result in Table 10 show reasonable predictions of shock for
a range of dynamic deflections. The accuracy of the prediction with
the Arcel data implies that accounting for heat transfer and
thermodynamics are necessary. More important the results show
that the polytropic assumption seems valid. To qualify the
polytropic process, extreme drop conditions and odd cushion shapes
need to be tested. This qualification will help justify the polytropic
process. As will be seen, the odd cushion shape to be test will be a
ribbed cushion. The results will then be compared to the equivalent
volume method [8].

3.3 Ribbed Cushion Results

For a ribbed cushion, A(x) and E(x) will vary in the ribs and
remain constant in the plank. The schematic for the ribbed cushions
are shown in Figure 13 and Figure 14. The meshing of the ribbed
cushion as shown in Figure 15 shows the finite element model with
one rib. Since the finite element model is one dimensional, the cross
sectional area at each position x must weighted (see Table 11). The
mesh chosen for the ribbed cushion is a one dimensional uniform
mesh. The length of each element is 0.5 inches. Further mesh
refinements were found to be unnecessary since the results yielded
reasonable values. The force vs compression results from the FEM
solution are shown in Table 12. Comparing Tables 9 and 12, one
notices that at a given compression, the energy absorbed by a block
cushion is greater than that by the ribbed cushion with the same
overall size. This is reasonable since it takes more energy to
compress a block cushion because there is more material. This is also

shown by the force compression curves in Figure 16.

55

 

 

 

 

 

 

Figure 13 Three dimensional view of the ribbed cushion

 

56

 

 

> 1.5"

 

 

4"

211

 

 

 

I 1.5"I

Figure 14 Two dimensional view of the ribbed cushion

57

—I

 

 

 

 

 

 

 

>—o—<>—<>—<><>

616.11%de
iU'i-h»uil\i

 

L 3:03 /
U9

Figure 15 Finite element mesh for the rib cushion

58

 

material modulus
Element number

Cross sectional area

 

 

 

 

 

 

 

 

 

 

 

 

 

(initial guess) E(x) A(x), (in2)
_(psi) _

1 14.7 64.0

2 14.7 64.0

3 14.7 64.0

4 14.7 64.0

5 14.7 28.9

6 14.7 27.5

7 14.7 26.1

8 14.7 24.7

 

 

 

Table 11 Geometry and initial modulus for the ribbed cushion

 

 

59

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

prescribed corresponding energy
compression force absorbed
(in) (lb) (in*lb)
0.2 1313.28 147.90
0.4 1350.40 422.92
0.6 1394.56 711.84
0.8 1442.56 992.99I
1 1502.08 1277.27
1.2 1569.28 1578.43
1.4 1649.92 1902.58
1.6 1745.28 2248.80
1.8 1861.12 2614.27
2 2003.84 2999.75
2.2 2181.12 3412.38
2.4 2408.32 3865.34
2.6 2705.92 4375.12
2.8 3112.32 4959.45
3 3697.28 5640.60]
3.2 4599.68 6460.50]
3.4 6147.84 7516.37I

 

Table 12 Force vs. compression results for a 4 inch thick ribbed
cushion with cross section dimensions as specified in Table 12

60

 

 

 

10000 -
7 5 00 -
E
8
g 5000—
2 5 00 -
"' """"" Ribbed Cushion
U Block Cushion
0 {.1 l I l l
0 1 2 3 4

compression (inch)

Figure 16 Comparison between ribbed and block cushions with the
same overall dimensions

61

The FEM solution of the ribbed cushion was also compared to
theoretical drop results [8]. Drops were performed at three different
drop heights with three different weights. The test drop heights
were 18, 30 and 42 inches with weights of 37.6 lb., 65.8 lb. and 94.2
lb.. The theoretical calculation involved finding the volume of a
ribbed cushion and converting it to an equivalent block cushion of
the same thickness. The comparison is shown in Table 13. Based
upon the error, it seems that the constitutive equation derived
earlier does model the closed cell cushion accurately. Since the equal

volume method is a proven technique, it may be said that the

 

polytropic assumption is a valid one.

The development of the polytropic process has shown that heat
transfer must be incorporated into the constitutive equation to model
cushions. The development has also shown that it is possible to lump
allithe thermodynamic terms into a constant k. As a result of the
polytropic process we are now able to model a ribbed cushion based

on numerical methods.

62

 

 

 

 

 

 

 

 

 

 

 

drop weight finite element equivalent percent
height (lb) results (G's) volume results error
(in) (G'S)

18 37.6 3 7 3 3 12
18 65.8 2 3 21 10
18 94.2 17 17 0

30 37.6 39 36 8

30 65.8 25 27 7

30 94.2 21 22 5

42 37.6 42 41 2

42 65.8 29 30 3

42 94.2 26 29 10

 

 

 

 

 

Table 13 Comparison between finite element results and the
equivalent volume results [8]

 

63

CONCLUSIONS

There were several steps involved in the modeling process of a
closed cell foam during an impact. The first model directly accounted
for heat transfer between the air and the foam. A system of three
differential equations obtained from Newton's Law, the Gas Law, and
the First Law of Thermodynamics were solved numerically by the
finite difference method (FDM). The prediction for peak G-level in a
drop was found to be in close agreement with published data. Based
on the success of this model, a simplified model which incorporated
all of the thermodynamic terms into a single property, the polytropic
constant k, was developed. This simplification converted the
problem into a non-linear, one dimensional elasticity problem which
was solved using the finite element method (FEM). The prediction
for peak G-level using this model was also found to be in excellent
agreement with published data.

To verify the accuracy of both the FEM and FDM models,
extreme drop scenarios must be tested. Table 14 shows the FDM
solution and the FEM solution for extreme drop situations. Both the
FEM and FDM model produced similar results as expected since FEM
is derived from the FDM model. The results in Table 14 shows that it
is necessary to consider thermodynamics to accurately model a
closed cell foam and that it is possible to simplify the analysis using a
polytropic process to represent the three differential equations. The
discrepancy in G at the 48 inch drop for the FDM model is probably

due to propagation of errors in the finite difference method.

64

 

 

 

 

 

drop cushion static . FDM FEM cushion
height thickness loading solution solution curve
(in) (in) (psi) (G's) (G's) data (G's)
12 4 .5 43 43 41

24 2 1 36 35 3 7

48 2 1 73 69 70

 

 

 

 

 

 

Table 14 Comparison between the finite difference, finite element,

and actual results at extreme loading conditions

 

65

As the impact becomes more severe, the FDM model becomes less
accurate due to the smaller time intervals needed to capture the
deceleration accurately. Overall, it is difficult to dispute the validity
of the model based on the results.

The theory that the force required to compress a closed- cell
foam is strain rate dependent seem to be justified. The rate
dependence however is not due to any visco-elastic properties of the
foam material but to heat transfer between the air and the cell walls.
The process of compression is close to being isothermal. The
accuracy of the results for a block cushions was verified using actual
published data. The results for ribbed cushions cannot be completely
verified since only a comparison to the equivalent volume method
could be made. The comparison however does show a good

correlation.

APPENDIX A

FINITE ELEMENT PROGRAM LISTING

66

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

variable list for one dimensional closed cell cushion block/ribbed

This program solves the linear equations of a second order
equation:
-a(X)(u’)’ + b(x)u = f(x)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
V A R I A B L E S
PARAMETERS:

IUNIT - INTEGER, UNIT NUMBER OF INPUT FILE.
OUNIT — INTEGER, UNIT NUMBER OF OUTPUT FILE.

BNDMAX - INTEGER, COLUMN DIMENSION OF GLOBAL STIFFNESS:

EDFMAX - INTEGER, MAXIMUM NUMBER OF D.O.F. PER ELEMENT.

NBCMAX - INTEGER, MAXIMUM NUMBER OF BOUNDARY CONDITIONS.

NDFMAX - INTEGER, MAXIMUM NUMBER OF D.O.F. PER NODE.

NEMMAX - INTEGER, MAXIMUM NUMBER OF ELEMENTS IN THE MESH.

NEQMAX - INTEGER, MAXIMUM NUMBER OF EQUATIONS (OR MAXIMUM
NUMBER OF D.O.F. IN THE MESH).

NGPMAX - INTEGER, MAXIMUM NUMBER OF GAUSS POINTS.

NNMMAX - INTEGER, MAXIMUM NUMBER OF NODES IN THE MESH.

NPEMAX - INTEGER, MAXIMUM NUMBER OF NODES PER ELEMENT.

ELEMENT PROPERTIES:

NDFELM - INTEGER, NUMBER OF D.O.F. PER ELEMENT.
NDFNOD - INTEGER, NUMBER OF D.O.F. PER NODE.
NODELM - INTEGER, NUMBER OF NODES PER ELEMENT.

MESH:

BNDWTH - INTEGER, HALF BANDWIDTH OF GLOBAL STIFFNESS MATRIX
FOR BENDING PROBLEMS.

EX - DOUBLE PRECISION(NPEMAX), EX(I) IS THE X-COORDINATE
OF THE I-TH NODE OF THE ELEMENT UNDER CONSIDERATION.

X

QIIIIOI)(IOi7fint7(10(7flntfiriordrinr5rintorinrdrinrvr)0rd()0(or)0t5(interinrdrintdriocdriordrinrdr>0r5C)nr5

IMESH INTEGER, FLAG TO INDICATE HOW MESH WILL BE INPUT:
IMESH = 0 IF THE MESH IS INPUT BY HAND,
IMESH = 1 IF THE MESH IS COMPUTED BY THE PROGRAM.
NDFMSH INTEGER, NUMBER OF D.O.F. IN THE MESH.
NELMSH INTEGER, NUMBER OF ELEMENTS IN THE MESH.
NODE INTEGER(NEMMAX,NPEMAX), CONNECTIVITY MATRIX:
NODE(I,J) IS THE J-TH NODE OF THE I—TH ELEMENT.
NODMSH INTEGER, NUMBER OF NODES IN THE MESH.
DOUBLE PRECISION(NNMMAX), X(I) IS THE GLOBAL

COORDINATE OF THE I-TH NODE.

BOUNDARY CONDITIONS:

ISGD INTEGER(NBCMAX,2), ISGD(I,J) INDICATES THAT THE J-TH
DISPLACEMENT D.O.F. OF THE I-TH NODE IS SPECIFIED.

ISGF INTEGER(NBCMAX,2), ISGF(I,J) INDICATES THAT THE J-TH
FORCE COMPONENT OF THE I-TH NODE IS SPECIFIED.

NSGD INTEGER, NUMBER OF DISPLACEMENT BOUNDARY CONDITIONS.

NSGF INTEGER, NUMBER OF NODAL FORCE BOUNDARY CONDITIONS.

VSGD DOUBLE PRECISIONINBCMAX), VSGD(I) IS THE SPECIFIED
VALUE OF THE I-TH DISPLACEMENT BOUNDARY CONDITION.

VSGF DOUBLE PRECISION(NBCMAX), VSGF(I) IS THE SPECIFIED

VALUE OF THE I-TH NODAL FORCE BOUNDARY CONDITION.

FINITE ELEMENT MODEL:

('3000OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO

C
C
C
C
C
C
C
C
C

C
C
C
C
C
C
C
C
C
C
C
C

DSF

EF
ESTIF
GAUSPT
GAUSWT
GDSF
GF
GSTIF
NGPSTF
NGPSIG

SF.

67

DOUBLE PRECISION(NPEMAX), DSF(I) IS THE LOCAL (XI)
DERIVATIVE OF THE I-TH SHAPE FUNCTION.

DOUBLE PRECISION(EDFMAX), ELEMENT FORCE VECTOR.
DOUBLE PRECISION(EDFMAX,EDFMAX), ELEMENT STIFFNESS.
DOUBLE PRECISION(NGPMAX,NGPMAX) GAUSPT(I,J) IS THE
I-TH GAUSS POINT OF J-TH GAUSS RULE.

DOUBLE PRECISION(NGPMAX,NGPMAX) GAUSWT(I,J) IS THE
I-TH GAUSS WEIGHT OF J-TH GAUSS RULE.

DOUBLE PRECISION(NPEMAX), GDSF(I) IS THE GLOBAL (X)
DERIVATIVE OF THE I-TH SHAPE FUNCTION.

DOUBLE PRECISION(NEQMAX), GLOBAL FORCE VECTOR,
WHICH, ON RETURN, CONTAINS THE SOLUTION VECTOR;
ALSO USED AS A DUMMY VECTOR IN EIGENPROBLEMS.
DOUBLE PRECISION(NEQMAX,BNDMAX), GLOBAL STIFFNESS.
INTEGER, NUMBER OF GAUSS POINTS TO BE USED TO
INTEGRATE THE STIFFNESS MATRIX.

INTEGER, NUMBER OF GAUSS POINTS AT WHICH THE
STRESSES ARE TO BE EVALUATED.

INTEGER(NPEMAX), SF(I) IS THE I-TH SHAPE FUNCTION.

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

68

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
Master Thesis Project Department of Packaging
By Andrew Chen
October 20, 1993
This program solves the linear equations of a second order
equation:

THIS PROGRAM WILL FIND THE STRESS, STRAIN AND YOUNGSMODULUS

OF A CLOSED CELL CUSHION BASED ON THE AIR MODEL AND THE
PROPERTIES OF THE MATERIAL SOLVED SEPERATELY

Note: This program will stabilize oscilating systems
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

00000000000

C
C
C
c
C
C
-a(x)(u’)' + b(x)u = f(x) C
C
C
C
C
C

C *** Put declaration, dimension and open statements here.

implicit double precision (a-h,o-z)

integer neqmax, bndwth, nemmax, npemax, nodelm, ndfnod,
node, nelmsh, ngpstf, ngpgrd, bndmax, ndfelm, count,
elmnum, edfmax, ndfmax, ngpmax, nnmmax, nbcmax, ires,
ni, i, j, nodmsh, ndfmsh, nsgd, nsgf, isgd, isgf, finish,
k

+ + + +

double precision vsgd, domain, elenth, gauspt, gauswt, pa,
sf, dsf, gf, estif, ef, gstif, jacobn,
elx, x, xi, dsx, gdsf, vsgf, nodsep,
x0, al, a2, a0, bl, b2, b0, fl, f2, f0, area, e, strain,
kon, stress,strntp,strnag

+-++ +

parameter (bndmax=20, edfmax=3, nbcmax=20, ndfmax:l, nemmax=20,
neqmax=20, ngpmax=4, nnmmax=20, npemax=3, pa=14.7,
finish=10)

dimension node(nemmax,npemax), elx(nnmmax), isgd(nbcmax,2),

+ vsgd(nbcmax), isgf(nbcmax,2), vsgf(nbcmax),
+ gauspt(ngpmax,ngpmax), gauswt(ngpmax,ngpmax),
+ gf(neqmax), gstif(neqmax,bndmax), ef(edfmax).
+ gdsf(npemax), estif(edfmax,edfmax),area(nemmax),
+ sf(npemax), dsf(npemax),a0(nemmax),e(nemmax),
+ strain(nemmax), stress(nemmax),
+ strntp(nemmax), strnag(nemmax)
open (unit = 10, file = ’dtl', status = ’old’)
open (unit = 11, file = 'sdtl', status = 'unknown')
C ____________________
C --— PREPROCESSOR --—
C ____________________

write(ll,*)’This is the data output for evaluating the strains &’
write(ll,*)’stresses of a 1 dimensional material that obeys the'
write(ll,*)’following equation’

write(ll,*)’ —a(x) u” + b(x) u = f(x)'

C *** Read input data here.
a1:
a2:
b0:
bl:
b2:
f0:
fl:
f2:

00000000

1 continue
read(10.*) kon

69

read(10,*) nodelm, ndfnod, ngpstf, ngpgrd, nelmsh
read(10,*) (e(n), n=l , nelmsh)

read(10,*) (area(n), n=1 , nelmsh)
read(10,*) a1, a2

read(10,*) b0, b1, b2

read(10,*) f0, fl, f2

read(10,*) nsgd, nsgf, domain

read(10,*) ((isgd(i,j), jzl, 2), i=1, nsgd)
read(10,*) (vsgd(i), i=1, nsgd)

read(10,*) ((isgf(i,j), j=1, 2), i=1, nsgf)
read(10,*) (vsgf(i), i=1, nsgf)

c --------- writing the input data ———————————————————————————————————
write(ll,*)’+++++++++++++++++++++++++++++++++++++++++++’
write(ll,*)‘the nonlinear constant used is ’, kon
write(ll,*)'the number of nodes per element is ',nodelm
write(ll,*)‘the number of elements per mesh is ',nelmsh
write(ll,*)‘the length of the domain is ',domain
do 2 n=1, nelmsh
write(ll,*)’initial modulus=’,e(n),'for element’,n
write(ll,*)‘initial cross sectional area =', area(n)
write(ll,*)' ............................................ '

2 continue

c-——- section to calculate the number of nodes in a mesh, the global

c coordinate and the boolean connectivity matrix

c ___________________________________________________________________

do 10 i=1, nemmax
do 20 jzl, npemax
node(i,j) = O
elx(i) = 0.0
strntp(i) = 0.0
strnag(i) = 0.0
20 continue
10 continue
nodmsh = nelmsh*(nodelm - 1)+l
ndfmsh = nodmsh*ndfnod
ndfelm = nodelm*ndfnod
bndwth = ndfelm
c—-——calculation of the coordinate system
do 30 i=1, nodmsh
elenth = domain/nelmsh
nodsep = domain/(nodmsh-l)
elx(i) = nodsep*float(i-l)
write (ll,*)'Global node',i,'at',elx(i)
30 continue
c write(ll,*)’nodmsh’,nodmsh,'elenth',elenth

c———boolean connectivity calculation

35
c

50

40
c
c

do 35 i=1, nodelm
node(l,i) = i
continue
write(ll,*)' nodelm is’,nodelm
do 40 i=2, nelmsh
do 50 j=1, nodelm
node(i,j)= node(i-1,j) + (nodelm - 1)
continue
continue
write(ll,*)‘the connectivity matrix local node vs global element'
do 55 i=1, nelmsh
write(ll,*) (node(i,i). i=l,nodelm)

 

7O

55 continue

C --- INITIALIZE THE ARRAYS WHICH CONTAIN THE GAUSS POINTS AND WEIGHTS.

CALL GQUAD (NGPMAX, GAUSPT, GAUSWT)

C *** Main processing operations or call to processor subroutine goes here.

do 5 count=1, finish
write(ll,*)':::==::::=:=;::===:::::::::::::;:::::::::::::::::::::'
write(ll,*)'this is iteration step number ',count

write (11 I *) ’ =:==:=:::::=========:=:====::::::::::::::::::=::::::: '

a1:
a2:
b0:
b1:
132:
f0:
fl:
f2:
X0 : 0.0

CDOCDCDOCDCDO

do 65 i =1, nemmax
a0(i)=0.0
65 continue
do 8 i = 1,edfmax
do 6 j = l,edfmax
estif(i,j)=0.0
ef(i)=0.0
6 continue
8 continue

16 continue
18 continue

do 60 elmnumzl, nelmsh
n = elmnum

a0(n) = e(n)*area(n)

call stif (elx, ndfelm, ngpmax, gauspt, gauswt, elenth,

+ nodelm, estif,ef,a0,b0,f0,a1,a2,bl,b2,fl,

+ f2, elmnum, nodmsh, const, x, jacobn, xi, ngpstf,nemmax)

call ASMBLE (NEQMAX, BNDMAX, EDFMAX, NEMMAX, NPEMAX, NODELM,
NDFNOD, NODE , ESTIF , BF , GSTIF , GF ,
. ELMNUM) '

60 continue

do 23 i = l,nodelm
do 26 j = 1,nodelm
write (11,*) 'estif of node,’,i,',',j,'is’,estif(i,j)
26 continue
23 continue

call fbndry (neqmax, nbcmax, ndfnod, nsgf, isgf, vsgf, gf)
call ubndry (neqmax, bndmax, nbcmax, ndfnod, ndfmsh, bndwth,
+ nsgd, isgd, vsgd, gstif, gf)

7 1

ires = 0
call solve (neqmax, bndmax, ndfmsh, bndwth, gstif,
+ gf, ires)
do 70 k=1, nodmsh
write (ll,*)'displacement of ',gf(node(k,l)),’at node',k,
+ ' at x coord ’,elx(k)
70 continue
C .....................
c —-— POSTPROCESSOR —-—
C _____________________

C *** Printing of solution and computation of gradients goes here.

do 80 j = l, nelmsh
do 90 hi = l, ngpgrd

xi = gauspt<ni,ngpgrd)

call shplld (npemax, nodelm, xi, elenth, sf, dsf, gdsf,
+ jacobn)

x = 0.5*elenth*(l+xi)+elx(j)

dsx=0.0

do 99 i = l,nodelm

dsx: dsx + gdsf(i)*gf(node(j,i))
99 continue

dsx
(strain(j)+strntp(j))/2

strain(j)
strnag(j)

 

90 continue
write(ll,*)' ________________________________________ ,

e(j) = (pa/strnag(j))*(((1-strnag(j))**(-kon))-1)
write(ll,*)‘youngs modulus is’, e(j),’ for element',j
stress(j) = e(j)*strnag(j)

write(ll,*)’ stress is ',stress(j)
write(ll,*)’ strain is ',strain(j)
80 continue
do 82 n=1, nemmax
strntp(n)=strain(j)
strnag(n)=0.0
strain(n)=0.0
stress(n)=0.0
82 continue
5 continue

STOP

END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c c
c subroutine: stif c
c purpose: calculates the stiffness matrix of the variational c
C statement c
c c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

subroutine stif (elx, ndfelm, ngpmax, gauspt, gauswt, elenth,
+ nodelm, estif,ef,a0,b0,f0,a1,a2,bl,b2,fl,
+ £2, elmnum, nodmsh, const, x, jacobn, xi, ngpstf,nemmax)

integer i, j, ndfelm, ngpmax, ngpstf, nodelm, elmnum

double precision elx, xi, ni, const, gauspt, gauswt,
+ elenth, ef,

+ jacobn, sf, estif, gdsf, dsf, x, ax,
+ bx, fx, ex,

72

+ a1, a2,°a0, bl, b2, b0, fl, f2, f0
dimension ef(3), gauspt(4,4), gauswt(4,4), estif(3,3),
+ gdsf(3), elx(SO), dsf(3), sf(3),a0(nemmax)
c ----- initialize all arrays

do 10 i=1, ndfelm
ef(i) = 0.0
‘do 20 j=1, ndfelm
estif(i,j) = 0.0
20 continue
10 continue
c ----- do loop on the number of gauspts

do 40 ni = 1, ngpstf
xi = gauspt(ni,ngpstf)

ax = 0.0
c ————— passing the value of x(i) to ex(i)
ex = elx(elmnum)
c ----- call shape function to evaluate the interpolation functions

call shplld (npemax, nodelm, xi, elenth, sf, dsf, gdsf,

+ jacobn)
const = (jacobn)*gauswt(ni,ngpstf)
x=0.5*elenth*(l.0 + xi) + ex
c
c ----- define the coefficient of the differential equation
c
ax = a0(elmnum) + al*x + a2*x**2
bx = b0 + bl*x + b2*x**2
fx = £0 + f1*x + f2*x**2
c
c ----- compute the coefficent matrix and the column vector for the deq
C .
do 50 i=1, ndfelm
ef(i) = ef(i) + const*sf(i)*fx
do 60 j=l, ndfelm
estif(i,j) = estif(i,j) + ax*const*gdsf(i)*gdsf(j) +
+ bx*const*sf(i)*sf(j)
60 continue
50 continue
40 continue
do 23 i = l, ndfelm
do 26 j = 1, ndfelm
c write (11,*) 'estif of node,',i,',',j,’is',estif(i,j)
26 continue
23 continue
return
end
c—c ----------------------------------------------------------------------
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C SUBROUTINE: ASMBLE C
C PURPOSE: ASSEMBLE THE GLOBAL MATRICES IN BANDED FORM. C
C C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE ASMBLE (NEQMAX, BNDMAX, EDFMAX, NEMMAX, NPEMAX, NODELM,
NDFNOD, NODE , ESTIF , EF , GSTIF , GF ,
ELMNUM)
IMPLICIT LOGICAL (A-Z)

INTEGER NEOMAX, BNDMAX, EDFMAX, NEMMAX, NPEMAX, NODELM,

73

NDFNOD, NODE , ELMNUM, INOD , JNOD , IDOF ,
JDOF , ROW , COL , COLl , EROW , ECOL
DOUBLE PRECISION ESTIF , EF , GSTIF , GF
DIMENSION NODE(NEMMAX,NPEMAX) , ESTIF(EDFMAX,EDFMAX), EF(EDFMAX),

GSTIF(NEQMAX,BNDMAX), GF(NEQMAX)

DO 40 INOD = l, NODELM
ROW = NDFNOD * (NODE(ELMNUM,INOD) - 1)

DO 30 IDOF = l, NDFNOD
ROW = ROW + l
EROW : (INOD - l) * NDFNOD + IDOF
GF(ROW) = GF(ROW) + EF(EROW)

DO 20 JNOD = l, NODELM
COLl : NDFNOD * (NODE(ELMNUM,JNOD) - 1)

DO 10 JDOF = 1, NDFNOD
COL : COLl - ROW + JDOF + l
ECOL : (JNOD - l) * NDFNOD + JDOF
IF (COL .GT. 0) THEN

GSTIF(ROW,COL) = GSTIF(ROW,COL) + ESTIF(EROW,ECOL)
END IF
10 CONTINUE
20 CONTINUE
3O CONTINUE
4O CONTINUE
RETURN
END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C SUBROUTINE: FBNDRY C
C PURPOSE: IMPOSE THE PRESCRIBED FORCE B.C.’S ON THE C
C BANDED SYMMETRIC STIFFNESS MATRIX AND MODIFY C
C THE FORCE VECTOR. C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE FBNDRY (NEQMAX, NBCMAX, NDFNOD, NSGF , ISGF , VSGF ,
GF )

10

IMPLICIT LOGICAL (A-Z)

INTEGER NEQMAX, NBCMAX, NDFNOD, NSGF , ISGF , IBC ,
IDOF

DOUBLE PRECISION VSGF , GF

DIMENSION ISGF(NBCMAX,2), VSGF(NBCMAX), GF(NEQMAX)
DO 10 IBC = l, NSGF
IDOF : NDFNOD * (ISGF(IBC,1) - l) + ISGF(IBC,2)
GF(IDOF) = GF(IDOF) + VSGF(IBC)
CONTINUE
RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

00000

SUBROUTINE: GQUAD
PURPOSE: ASSIGN VALUES TO THE ARRAYS "GAUSPT" AND "GAUSWT“
USED IN GAUSSIAN QUADRATURE.

00000

74

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C ___

10

20
C -__
C --_
C ---
C -_-

SUBROUTINE GQUAD

IMPLICIT LOGICAL (A-Z)

INTEGER

NGPMAX,

DOUBLE PRECISION GAUSPT,

DIMENSION

GAUSPT(NGPMAX,NGPMAX),

INITIALIZE ARRAYS.

DO 20 I = l,NGPMAX

DO 10 J = l,NGPMAX
GAUSPT(I,J) = 0.000
GAUSWT(I,J) = 0.000

CONTINUE

CONTINUE

(NGPMAX, GAUSPT, GAUSWT )

GAUSWT

GAUSSIAN QUADRATURE OF ORDER 1.

GAUSPT(l,l)
GAUSWT(1,1)

0.000
2.000

GAUSSIAN QUADRATURE OF ORDER 2.

GAUSPT(1,2)
GAUSPT(2,2)
GAUSWT(1,2)
GAUSWT(2,2)

H I! II II

1.000
GAUSWT(1,2)

—l.OOO/SQRT(3.000)

GAUSSIAN QUADRATURE OF ORDER 3.

GAUSPT(1,3)
GAUSPT(2,3)
GAUSPT(3,3)
GAUSWT(1,3)
GAUSWT(2,3)
GAUSWT(3,3)

0.000

5.000/9.000
8.000/9.000
GAUSWT(1,3)

~SQRT(3.000/5.000)

GAUSSIAN QUADRATURE OF ORDER 4.

GAUSPT(1,4)
GAUSPT(2,4)
GAUSPT(3,4)
GAUSPT(4,4)
GAUSWT(1,4)
GAUSWT(2,4)
GAUSWT(3,4)
GAUSWT(4,4)

RETURN
END

—O.86ll363116
-O.3399810436
-GAUSPT(2,4)
-GAUSPT(1,4)
0.3478548451
0.6521451549

GAUSWT(2,4)
GAUSWT(1,4)

GAUSWT(NGPMAX,NGPMAX)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

000000

SUBROUTINE:

PURPOSE:

SUBROUTINE SHPLlD (NPEMAX, NODELM, XI
JACOBN)

SHPLlD
EVALUATE THE 1-D LAGRANGIAN INTERPOLATION FUNCTIONS C
AND THEIR GLOBAL DERIVATIVES AT THE GAUSS POINTS.

GDSF ,

ELENTH, SF

C
C

C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C -.._

C ___

C .___

10

75

IMPLICIT LOGICAB (A-Z)

INTEGER NPEMAX, NODELM, I
DOUBLE PRECISION XI , ELENTH, SF , DSF , GDSF , JACOBN
DIMENSION SF(NPEMAX), DSF(NPEMAX), GDSF(NPEMAX)

LINEAR LAGRANGE INTERPOLATION FUNCTIONS FOR Z-NODED ELEMENTS.

IF (NODELM .EQ. 2) THEN

SF(l) = 0.5000 * (1.000 — XI)
SF(Z) = 0.5000 * (1.000 + XI)
DSF(l) = -0.5000
DSF(Z) = 0.5000

QUADRATIC LAGRANGE INTERPOLATION FUNCTIONS FOR 3-NODED ELEMENTS.

ELSE IF (NODELM .EQ. 3) THEN

SF(l) = ‘0.5000 * XI * (1.000 - XI)

SF(Z) = 1.000 - XI * XI

SF(3) = 0.5000 * XI * (1.000 + XI)

DSF(I) = '0.5000 * (1.000 - 2.000 * XI)

DSF(Z) : -2.000 * XI

DSF(3) = 0.5000 * (1.000 + 2.000 * XI)
END IF

COMPUTE THE GLOBAL DERIVATIVES OF SF(I).

JACOBN : ELENTH * 0.5000
DO 10 I = 1, NODELM
GDSF(I) = DSF(I) / JACOBN
CONTINUE
RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C
C
C
C
C
C

SOLVE:

THIS PROGRAM SOLVES A BANDED SYMMETRIC SYSTEM OF EQUATIONS.
THE BANDED MATRIX IS INPUT THROUGH BAND(NEQNS,NBW), AND

RHS IS THE RIGHT HAND SIDE (FORCE VECTOR) OF THE EQUATION.
NEQNS IS THE NO..OF EQUATIONS AND NBW IS THE HALF BAND WIDTH.
IN RESOLVING, IRES .GT. 0, LHS ELIMINATION IS SKIPPED.

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE SOLVE (NRM , NCM , NEQNS , NBW , BAND ,
RHS , IRES )

IMPLICIT DOUBLE PRECISION (A-H,O-Z)

INTEGER NRM , NCM , NEQNS , NBW , IRES
DOUBLE PRECISION BAND , RHS

DIMENSION BAND(NRM,NCM), RHS(NRM)

MEQNS:NEQNS-l

IF (IRES.GT.0) GO TO 40

DO 30 NPIV=1,MEQNS

NPIVOT:NPIV+1

LSTSUB=NPIV+NBW—l

IF (LSTSUB.GT.NEQNS) LSTSUB=NEQNS
DO 20 NROW=NPIVOT,LSTSUB

C
C
C
C
C
C
C
C
C

76

C --- INVERT ROWS AND COLUMNS FOR ROW FACTOR.

NCOL=NROW-NPIV+1
FACTOR:BAND(NPIV,NCOL)/BAND(NPIV,1)
DO 10 NCOL=NROW,LSTSUB
ICOL=NCOL-NROW+1
JCOL=NCOL—NPIV+1
10 BAND(NROW,ICOL):BAND(NROW,ICOL)—FACTOR*BAND(NPIV,JCOL)
20 RHS(NROW):RHS(NROW)-FACTOR*RHS(NPIV)
30 CONTINUE
GO TO 70
40 DO 60 NPIV=1,MEQNS
NPIVOT=NPIV+1
LSTSUB=NPIV+NBW~1
IF (LSTSUB.GT.NEQNS) LSTSUB=NEQNS
DO 50 NROW=NPIVOT,LSTSUB
NCOL=NROW-NPIV+1
FACTOR:BAND(NPIV,NCOL)/BAND(NPIV,l)
50 RHS(NROW)=RHS(NROW)-FACTOR*RHS(NPIV)
60 CONTINUE

C --- BACK SUBSTITUTION.

70 DO 90 IJK:2,NEQNS
NPIV=NEQNS-IJK+2
RHS(NPIV):RHS(NPIV)/BAND(NPIV,1)
LSTSUB=NPIV-NBW+1
IF (LSTSUB.LT.1) LSTSUB=1
NPIVOTzNPIV-l
DO 80 JKI:LSTSUB,NPIVOT
NROW=NPIVOT-JKI+LSTSUB
NCOL=NPIV-NROW+1
FACTOR:BAND(NROW,NCOL)

80 RHS(NROW)=RHS(NROW)-FACTOR*RHS(NPIV)

90 CONTINUE
RHS(1)=RHS(l)/BAND(1,1)

 

RETURN

END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C SUBROUTINE: UBNDRY C
C PURPOSE: IMPOSE THE PRESCRIBED DISPLACEMENT B.C.‘S ON THE C
C BANDED SYMMETRIC STIFFNESS MATRIX AND MODIFY C
C THE FORCE VECTOR. C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE UBNDRY (NEQMAX, BNDMAX, NBCMAX, NDFNOD, NDFMSH, BNDWTH,
. NSGD , ISGD , VSGD , GSTIF , GF )

IMPLICIT LOGICAL (A-Z)

INTEGER NEQMAX, BNDMAX, NBCMAX, NDFNOD, NDFMSH, BNDWTH,
NSGD , ISGD , IBC , IDOF , BWMl , ROW ,
COL , K
DOUBLE PRECISION VSGD , GSTIF , GF , UVALUE
DIMENSION ISGD(NBCMAX,2), VSGD(NBCMAX),
GSTIF(NEQMAX,BNDMAX), GF(NEQMAX)
DO 40 IBC : l, NSGD
IDOF : NDFNOD * (ISGD(IBC,1) - l) + ISGD(IBC,2)
UVALUE : VSGD(IBC)

77

BNDWTH - 1
IDOF '- BNDWTH

BWMl
ROW

DO 20 K = l, BWMl
ROW : ROW + 1
IF (ROW .GT. 0) THEN
COL : IDOF - ROW + 1
GF(ROW) : GF(ROW) - GSTIF(ROW,COL) * UVALUE
GSTIF(ROW,COL) : 0.000

END IF
20 CONTINUE
GSTIF(IDOF,1) = 1.000
GF(IDOF) : UVALUE
ROW : IDOF

DO 30 K = 2, BNDWTH
ROW : ROW + 1
IF (ROW .LE. NDFMSH) THEN
GF(ROW) : GF(ROW) ‘ GSTIF(IDOF,K) * UVALUE
GSTIF(IDOF,K) = 0.000
END IF
30 CONTINUE
40 CONTINUE
RETURN
END

APPENDIX B

ARCEL 5.12 CUSHION CURVES

79

  
   

ARCO Chemical Company CARGO <>i

DMsion o! AllanlicRichlieldCanpany

 

Dynamic Cushioning
Performance

ARCEIE

Moldable Polyethylene Copolymer

 

 

80

 

 

Dynamic Cushioning Performance

 

 

Following are dynamic 12' DROP
cushioning performance 723.3... 1.. m... 0......- “ch 72%;. 2.5 mm. 0......- 2.0Pcr
curves for foam molded : .

from

ARCELG' 512

Moldable Polyethylene
Copolymer

Deceleration. G‘s
Deceleration. G‘s

 

0.5 1 5 2.0 2.5 3.0 0.5

1 .0 . 2.0 2.5 8.0
Static Stress. psi

1 .O t .5
Static Stress. psi

 

18" DROP
Flour- 3 Figure 4
18" Drop. tat Impact Denalty = 2.0 PCF t8" Drop. 2~5 Impacts Denalty = 2.0 PCP

 

1 1

Deceleratlon. G's
Deceleration. G's

 

 

0.5 1.0 t .5 2.0 2.5 3.0 0.5 1.0 1.5 2.0 2.5 3.0
Static Stress. pal Static Streaa. pol
24" DROP
Flour. 5 Flow 6
24" Drop. let Impact Oenalty =- 2.0 POP 24" Drop. 2-5 Impacts Denalty I 2.0 PCF

1 1

Deceleration. G's
Deceleration. G'a

 

0.5

2.0 2.5 3.0 0.5 2.0 2.5 3.0

1 .0 1 .6
Static Stress. pal

t .0 t .6
Static Stress. pat

© ‘Al .1 —- Ll‘ a4“
19$"..- .- ._ ,

 

81

 

ARCEL 512 Resin

 

 

30" DROP 48" DROP
Home 7 Home I figure 13
30" Drop. 1st Impact Density - 2.0 POP 30" Drop. 2-5 Impacts Denslty - 2.0 POP 48" Drop. 1st Impact Denatty - 2.0 PCF

1.5"

 

 

 

 

 

 

 

.- .- .a

o o o

. :- 6

§ 2 2

‘2' E E

2 2 2

O O 0

0 0 o

0 0 0

O O o

0.5 1.0 1.5 2.0 2.5 3.0 0.5 1.0 1.5 2.0 2.5 3.0 0.5 1.0 1.5 2.0 2.5 3.0
Static Stress. psi Statlc Stress. psl Statlc Stress. ps1

36" DROP
Flame 9 no.“ 10 Figure 14
36" Drop. 1:: Impact Density = 2.0 PCF 36" Drop. 2~s Impacts Density 2 2.0 PCF 48“ Drop. 26 tmpocts Donslty = 2.0 PCP

1 1 1::

1.0" k /3—°'
l UV \J u.
/ w
.. .0 9 W I '
o o o /
c' c' c'
2 2 2
‘g’ E E 110-1 w
2 2 2
3 3 a / Lu,-
5 8 a —/ /
A g

1
o;

 

 

 

 

 

 

 

 

 

 

 

0.5 1.0 1.5 2.0 2.5 3.0 0.5 1.0 1.5 2.0 2.5 3.0 0.5 2.0 2.5 3.0
Static Stress. psl Static Stress. psi

1 .0 1 .5
Statlc Stress. ps1

 

42" DROP
Home 11 Figure 12
42“ Drop. 1st Impact 000th - 2,0 pop 42” Drop. 2-5 Impacts Density x 2.0 PCF

1

Deceleratlon. G's
Deceleration. G's

 

0.5 1

2.0 2.5 3.0 0.5 1 5 2.0 2.5 3.0

.0 1.5 1.0 .
Static Stress. ps1 Static Stress. psi

82

LIST OF REFERENCES

l. Averill, Ronald, Assistant Professor, Michigan State University,
Lecture material from Material Science and Mechanics Course
MSM 805. September 1992.

2. Burgess, G.J. " Some Thermodynamic Observations on the
Mechanical Properties of Cushions ",J. Cell. Plast..,24, 57-69 (Jan.
1988)

3. Burgess, G.J." Consolidation of Cushion Curves
and Sci... 17, 189-194 (1990)

4. Branenburg, Richard and Julian Lee Fundamentals of Packaging
Dynamics, Skaneateles, NY: L.A.B (1991)

5. Chen-Wai-Fah and Atef F. Saleeb. Constitutive Equations For
Engineering Materials, Volume 1 Elasticity and Modeling, New
York: John Wiley and Sons (1982)

6. Devore, Jay L. Probability and Statistics for Engineering and the
Sciences, Second Edition, Montery, CA: Brooks/ Cole Publishing
Company (1987)

7. Gerald, Curtis, and Patrick 0. Wheatley. Applied Numerical
Analysis, Sixth Edition, Reading, MA: Addison-Wesley Publishing
Company (1989)

8. Granthen, Gary "Predicting Shock Transmission Characteristics for
Ribbed Expanded Polypropylene Cushions Using Standard
Curves for Flat Plank Cushions", MS Thesis, School of Packaging,
Michigan State University (1991)

9. Hilyard, M.C. Mechanics of Cellular Plastics,, New York: Applied
Science Publishers (1982)

10. Holman,J.P. Heat Transfer, Sixth Edition, New York: McGraw-Hill
Book Company (1986)

11. Koffman, Elliot B., and Frank L. Friedman. Fortran 77, Third
Edition, Reading, MA: Addison-Wesley Publishing Company
(1987)

12. Kreyszig, Erwin. Advanced Engineering Mathematics , Fourth
Edition, New York: John Wiley and Sons (1979)

13. Larson, Roland E., Robert P. Hostetler,CalculusWith Analytic
Geometry, Lexington, MA: D.C. Heath Company (1982)

14. Mantonis, V. (1964) SPE J., 20, 1024

, Packaging Tech.

 

15.

16.

17.

18.

19.

20.

21.

22.

23.

24.

25.

26.

27.

Popov, Egor P. Engineering Mechanics of Solids, Englewood Cliffs,
NJ: Prentice Hall (1990)

Product Data: "Dynamic Cushioning Performance", Arcel Moldable
Polyethylene Copolymer, ARCO Chemical CO, Division of Atlantic
Richfield CO, Philadelphia, PA 19101

Reddy,J.N. An Introduction to the Finite Element Method, New
York: McGraw-Hill Publishing Company (1984)

Robertson, John A., and Clayton T. Crowe. Engineering Fluid
Mechanics, Third Edition, Boston: Houghton Mifflin Company

(1985)

Segerlind, Larry J. Applied Finite Element Analysis, Second
Edition, New York: John Wiley and Sons (1984)

Shuttleworth, R., V. Shestopal and P. Gross." Open Cell Flexible
Polyurethane Foams: Comparison of Static and Dynamic
Compression Properties", J. Appl. Polymer Sci,30, 333-343
(1985)

Throne, J. and R Progelhof." Closed Cell Foam Behavior Under
Dynamic Loading-LStress-Strain Behavior of Low Density
Foams." Journal of Cellular Plastics (Nov./Dec. 1984).

Throne, J. and R Progelhof."Closed Cell Foam Behavior Under
Dynamic Loading-II. Dynamics of Low Density Foams." Journal
of Cellular Plastics (Jan/Feb. 1985).

Tipler, Paul A. Physics,, Second Edition, New York: Worth
Publishers Inc. (1982)

Totten, Troy L, G.J. Burgess, and S. Paul Singh."The Effects of
Multiple Impacts on the Cushion Properties of Closed-Cell
Foams",Packaging Tech. and Sci.,3, 117-122 (1990)

Wark, Kenneth, Jr. Thermodynamics, Fifth Edition, New York:
McGraw-Hill Book Company (1988)

Willems, Nicholas, John T. Easley and Stanley T. Rolfe. Strength of
Materials, New York: McGraw-Hill Book Company (1981)

"Standard Test Method for Shock Absorbing Characteristics of
Package Cushioning Material", ASTM Stan‘dards, D1596, Selected
AST M Standards on PACKAGING, Second Edition, Philadelphia,
PA (1987)

 

 

 

F.

‘LIBRQRIE
ll l ll
1 1 l

lj

“M
if

TQTE UNIV

ml *
312M301l3l55957

‘ ll l

N
9
C
II
H
C
.1
m