THESES

SITY LIBRARIES

Illlllllllllllllllllllllllll lol Will

3 12930

 

 

 

 

 

 

 

 

 

This is to certify that the

thesis entitled

DYNAMIC MODELING OF CLOS- CELL CUSHIONING MATERIAL

presented by

DHITRI V. POKUDIN

has been accepted towards fulfillment
of the requirements for

MASTERS degree in PACKAGING

 

 

GARY BURGESS QM Hwyk

Major professorq

Date 5722/98

 

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

 

 

 

LIBRARY

Michigan State
University

 

 

 

PLACE iN RETURN Box
to remove this checkout from your record.
To AVOID FINE-5 return on or before date due.

 

MTE DUE -. MTE DUE DATE DUE

 

 

 

NJ 2000

 

 

 

 

with» '

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

DYNAMIC MODELING OF CLOSED CELL
CUSHION ING MATERIAL

By

Dmitri V. Pokudin

THESIS

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

MASTER OF PACKAGING SCIENCE

School of Packaging

1998

ABSTRACT
DYNAMIC MODELING OF CLOSED-CELL CUSHIONING MATERIAL
By

Dmitri V. Pokudin

Closed cell cushioning materials are clearly a multi-phase material system, thus it is
not surprising that all attempts to describe their behavior by Hooke’s Law and the Gas
Law which idealize a single phase material response have failed. The approach taken in
this paper separates the polymer matrix and gas contributions, applying appropriate
idealizations to each of them. The model does not try to preconceive the simplest
idealizations and parameters but rather is used to study the system to find out which
phenomena are the most significant and to measure the values of the associated variables.

The study showed that representation of the trapped air as an ideal gas with some
heat transfer to the walls and the polymer matrix as a plastic material with an additional
rate dependent resistance gives the best results. The model’s search algorithm allows all
the required system parameters to be found from a single experimental acceleration
profile. These are used to build a modeled acceleration profile aimed at replacing the
Fourier filtering technique and answer all of the other questions about the cushion
behavior. Since the measured coefficients are fairly stable for a given material, cushion
thickness, drop height and weight, more drops can be processed to obtain and tabulate
their average values. Obtained data is similar in its description of the material properties to
the cushion curves which predict the peak G, but it is capable of providing the entire

impact profile and complete information about it.

ACKNOWLEDGEMENTS

This thesis would not have been possible the great input from Dr. Gary Burgess.
His knowledge of packaging dynamics, physical concepts, patience and miraculous ability
to dismiss my wrongful but very convincing arguments has led to the successful
completion of this thesis.

iii

TABLE OF CONTENTS

LIST OF TABLES
LIST OF FIGURES
CHAPTER 1 INTRODUCTION AND LITERATURE REVIEW
CHAPTER 2 MATERIALS AND METODS
2.1 Philosophy and justification of the modeling approach
2.2 Model assumptions
2.3 Derivation of the model’s physical foundation
2.4 Study and verification of the model’s assumptions
2.4.1 Heat transfer considerations
2.4.2 Defining the time of real cushion contact
2.4.3 Fv considerations
2.4.4 FX and X0 considerations
2.4.5 Justification of the model’s use as a filtering tool
2.4.6 Verification of suggested cushion recovery time
2.4.7 Consideration of model’s vertical leading edge
CHAPTER 3 RESULTS
3.1 Corrections and additions to the model assumptions
3.2 Model use for tabulating the material properties
3.3 Model use for fit filtering
CONCLUSIONS
APPENDIX

BIBLIOGRAPHY

PAGE

vi

l2

16

19

20

21

3O

31

32

34

35

36

38

39

48

LIST OF TABLES

Table 1 Fitted results for Elastic Fx model

Table 2 Fitted results for Plastic P" model

Table 3 Fitted results for Plastic F" model with adjusted X0
Table 4 Fitted results for Elastic F" model with adjusted Xo
Table 5. Results of fitting original and noisy signals

Table 6. Interpolated parameters for an arbitrary drops

PAGE

25

26

28

29

31

35

LIST OF FIGURES

Figure 1. Typical acceleration profile
Figure 2. Under-filtering of the accelerometer signal
Figure 3. Over-filtering of the accelerometer signal
Figure 4. Basic system set up
Figure 5. Adiabatic and Isothermal model checks
Figure 6. Introduction of heat transfer coefficient h=22 and h=43
Figure 7. Heat transfer models acceleration vs. X
Figure 8. Addition of Fv to explain the reaction force

at the time of contact
Figure 9. F" adjusted for the best fit
Figure 10. First estimation of possible Fx (Elastic)
Figure 11. Second estimation of possible 1:“x (plastic)
Figure 12. Fitted model curve with elastic Fx
Figure 13. Worst fit for the model curve with elastic Fx
Figure 14. Non physical result for plastic 1:"x
Figure 15. Fit filtering of a noisy signal
Figure 16. Model Predictions for {17.8 kg 0.557m} and

{13.25 kg 0.6611m} drops

Figure 17. Unknown signal fit filtering

PAGE

16

17

19

20

21

22

24

25

26

27

3O

36

37

CHAPTER 1
INTRODUCTION AND LITERATYRE REVIEW

Various types of cushioning materials are used in the packaging industry to protect
fragile products from damaging levels of acceleration. The mechanism of theirllactions is
to provide separation and some degree of independent movement to the product inside its
individual shipping container. By doing so, it is ensured that the shocks experienced by the
container are transmitted to the product through the cushion which both absorbs and
spreads part of the shock energy over time. Depending on the severity of the expected
shocks and the product fragility, different amounts of different cushion materials provide
the lowest cost for an adequate protection level. One of the most popular and frequently
used materials are cellular polymers.

Cellular polymers used as cushion materials are multi—phase systems that consist
of a polymer matrix and a fluid phase. The fluid is generally gas (air) and is either trapped
or continuous in the polymer matrix [3]. Thus two major types of polymer cushions exist:
closed-cell and opened—cell. The two types of material differ mainly in the gas phase
contribution to the cushion reaction force during its compression. Closed-cell air is
compressed and reacts back with increased pressure whereas an open-cell air escapes and
provides little or no resistance. There is much variation in the material properties even
within the types themselves. Cushion properties depend on the exact composition of the
polymer matrix, average cell size and wall thickness. Hundreds of different cushion grades
are available.

In order to utilize the advantages of one or the other material effectively and

minimize the cost, packaging engineers have to consider and evaluate a lot of different

package configurations. It is hardly ever practical to make and test prototypes for all the
possible choices. Instead, preliminary estimations are made to narrow the field of cushion
materials. Because of the complexity of closed-cell polymers, which I will now
concentrate on, none of the attempts to describe their properties in terms of a few
idealized phenomena and associated constants were successful. Hooke’s Law, which is
useful for most solid state materials, does not work [4]. Neither does purely adiabatic or
isothermal compression of the trapped air [1]. So in their calculations engineers have to
rely on the material properties supplied by the cushion manufacturer in the form of cushion
curves. Cushion curves relate the peak G transmitted by the cushion to its thickness,
product weight and bearing area for different expected drop heights. Since engineers are
faced with hundreds of similar tables and because no automation of the process can be
implemented, it becomes very difficult to use them.

As described in [5], each point on the cushion curve has to be obtained from an
average of 5-10 drops of a flat massive plate (drop head) directed squarely by guide rods
onto the sample cushion. In all, thousands of drops have to be performed and drop head
acceleration profiles analyzed. 80 the cushion curves are very difficult to generate too.

A typical acceleration profile recorded by the accelerometer mounted on the drop

head is shown in Figure l.

 

Acceleration [G]

 

 

 

 

Tlme [ms]

Figure 1. Typical acceleration profile.

The only parameter of interest for the cushion curves is peak acceleration.

Although the general shape of the signal is mainly ignored, much attention is given to its
fine structure (ripples). Their presence makes it difficult to determine the exact location
and value of the peak acceleration. Moreover since the real acceleration of the massive
drop head is expected to be smoother than that, the ripples are attributed to the equipment
noise and are usually filtered out before any of the impact parameters are measured. Some
of the recorders used in industry have hardware filtering and others do it in their software.
Both methods are very similar in their physical nature. Software algorithms expand the
real profile into a Fourier series and drop the trailing terms starting from some preselected
filtering frequency [6]. Hardware attenuates high frequencies by analog filters. The result
produced is very similar and obviously depends greatly on the choice of the filtering

frequency. The mathematical apparatus of Fourier filtering does not have anything to do

with the physical phenomena of the process at hand which consequently does not provide
any information on the proper filtering frequency to use. Excessively high frequencies do
not change the signal very much and may leave ambiguity in peak G readings (Figure 2.)
Whereas low filtering frequencies distort the signal too much and may render its other
parameters such as impact duration, restitution coefficient and area under the curve

needed by some more advanced applications useless (Figure 3.)

 

Acceleration [G]
G

 

 

 

 

 

Tlme [ms]

Figure 2. Under-filtering of the accelerometer signal.

 

B It

Acce|eret Ion [G]
a

n—
a
I

 

 

 

 

 

Time [ms]

Figure 3. Over-filtering of the accelerometer signal.

As a result, data on the cushion material properties provided by existing signal analysis
techniques is incomplete and possibly unreliable.

The questions which prompted this thesis are whether or not the process of
collecting and presenting the data on cushion material properties can be simplified.
I wanted to avoid conventional techniques of signal processing and look the entire
unfiltered profile in order get more information about the cushion properties, to possibly

develop a better method of filtering and loss of important signal features.

CHAPTER 2

MATERIALS AND METHODS

2.1 Philosophy and justification of the modeling approach

 

In earlier investigations of closed-cell materials, cushions were modeled as a
volume of trapped air [1], and as solid elastic and plastic materials [7] and [8]. Theoretical
formulas to determine peak accelerations and impact durations were worked out. Obtained
results and generated cushion curves were reasonably close to the experimental data. It
should be pointed out that ill spite of the apparent success, neither of the models in fact
accurately described the process and ind it was pointed out that some complicated
approach involving continuous heat transfer has to be taken.

In my approach, I am separating the polymer matrix and gas contributions and
applying appropriate idealizations to each of them. To maximize the amount of
information obtained from an experimental signal, I am looking not only at the peak 6’8
but at its entire profile. The model dynamically traces all the cushion parameters during the
entire impact duration. In order to do this a computer model utilizing a small time step
approach was created.

The philosophy behind modeling approach consists of the following: the model
attempts to describe the entire process of impacting a cushion with a drop head in terms of
known physics laws and the cushion physical characteristics. In order to do this, the
physical system is reasonably simplified and defined by a set of model assumptions.
Logical guesses and assumptions are made about the unknown phenomena and some
characteristics which are unavailable or difficult to measure. Subsequent running of the

model with one or another set of unknown variables and comparing the result to the real

shock profile verifies whether the assumptions were right and provides reasonable ranges
for the unknown values of model variables. Finally, when the model is well studied and is
believed to predict real system behavior with sufficient accuracy, a search algorithm is
employed to find a set of variables values which minimizes the mean square difference
between the model and real profile. The obtained curve is essentially a least square model
fit and can be used in place of the filtered curve to analyze the impact. On the other hand,
the model parameters fully describe everything and the curve is not really needed. The
approach is similar to the linear regression where some presumably linear data is fitted a

line to obtain its slope which is then used in place of the original data.

2.2 Model assumptions

1) the basic system set up is shown in Figure 4.

 

 

 

 

 

 

 

 

 

 

 

 

Drop head
M /
Drop Bearing Area
height /i/ _
Cushion To, Po 7 X_apparent

 

 

 

 

Figure 4. Basic system set up

2)

3)

4)

5)

6)

The object of investigation is a closed cell cushioning material. Its force producing
components are trapped air and matrix of PE bubble walls.

Cushion cross section area stays constant during compression, thus reducing the
problem to one dimension.

There is no loss of thermal energy to the environment. The low air-to-PE convection
heat transfer coefficient (order of lOW/mzK) and the fact that average duration of the
impact event is only 20- 40 milliseconds makes the loss negligible.

Initially, internal air is in thermal equilibrium with its environment and has normal
atmospheric pressure.

All trapped air is considered to be a single volume of ideal gas. Van der Waals
corrections for initial and maximum compression states result in less than 0.1%

difference compared to the Ideal Gas Law.

7) The heat transfer to the PE walls is not negligible because of the big contact area
between internal air and walls[ 1].

8) The internal energy and thus the temperature of the walls changes as a result of the
heat transfer and work done by the friction forces. The walls are much more massive
than the air and have greater heat capacity, thus the change is small. (Even if all the
energy of the falling mass is absorbed by the PE only its temperature will be raised
only by approximately 1°C. For the purpose of air behavior, this can be disregarded,
but is kept later to estimate cushion recovery times.)

9) The matrix of PE walls resist change of their height with forces proportional to the
rate and the rate squared during the entire impact, designated by Fv(t). They also resist
with some additional force as a function of relative compression itself, Fx(t).

10) The guide rods of the testing equipment can also have rate dependent friction, which is
indistinguishable from the force inside the cushion and is taken in account by the same
Fv(t) term.

11) Some of the cushion thickness consists of damaged cells. This is taken into account by
decreasing effective thickness.

12) For the entire time of the impact drop head face and the cushion surface stay together.
The end of an impact is defined by the moment of drop head and cushion separation. If
some of the friction comes from equipment, head acceleration at that moment will not
be equal to zero. So the cushion condition for the end point is P(t,nd) = P0.

13) Model input data consists of system constants describing major parameters of the
physical set up including test cushion dimensions:

0 Xwn - Cushion thickness.

Po - Initial equilibrium air pressure.

Area — Cushion area.

M - Mass of the drop head.

m - Cushion or PE mass.

H - effective drop height calculated from experimental impact velocity.
T. - Initial equilibrium temperature.

C., p - Heat capacity and density for air.

Cp- Heat capacity for PE.

14) The search algorithm needs AW vs time - experimental drop head

acceleration profile.

15) The following model variables are inputted or searched for by the algorithm for each

individual test:

Fx, Fv - Resistance forces of the PE matrix.

h - PE-to-Air convection heat transfer coefficient.

16) For given system constants and model variables, the model generates:

A(t) Predicted acceleration profile. Time positions of any point of interest
(peak G, duration, etc.) can be precisely read right from the graph.

Fit coefficient - mean square deviation between A(t) and Ammo).
Gm - Peak acceleration in G.

t(G..,)- Time position of Gm

Xn- Drop head position above the base plate at the moment of maximum

compression.

10

0 X“...- drop head position at the moment of its separation from the cushion.

0 t(T.=T.) - time when temperature of the air is equal to the temperature of the
walls. Walls start to cool after that.

0 t(V=0)- time of the moment when direction of velocity reverses (same as
maximum compression).

0 dTm- given for air and walls as the maximum temperature change reached by
them over the initial conditions.

0 dTu- given for air and walls as the temperature change over the initial
conditions at the end of the intact.

0 Restitution Coefficient— Final velocity of drop head divided by its initial
velocity.

15) The model uses (It - variable time step in order to find a trade off between the

model accuracy, execution time and round off error.

11

2.3 Derivation of the model’s physfial fourgation

The physics of the small step model is based on starting with initial conditions and
propagating them in time by adding changes incurred on them over the small period of
time. The major parameters of the system are the position, velocity and acceleration of the

cushion drop head boundary. The initial values are:

A(O) = -g
V(0)=-,/2-g-H (1)
X(0)= xwm

Where g is acceleration due to gravity. Due to extreme smallness of chosen time step (it,
the total force on the drop head over this period is considered to be constant.
Acceleration, velocity and position of the drop head at any moment of time t can be
expressed using Newton’s second law and its basic integral forms together with
acceleration, velocity , position and total force at the t-dt moment of time. It is convenient

to denote position t-dt with an integer n and t as (n+1).

 

Halal
ADM-I = M
Vn+l=Vn+An.dt (23, b,C)

 

-dt2
Xn+l= Xn+Vn-dt+ A"2

Repeating the calculations the necessary amount of times, we can obtain position, velocity
and acceleration as a discrete functions of t=0,dt,2dt,3dt,...ndt for the entire impact
duration. The total upward force is

Fl“ = 1'7." + EV + F." - 1%- Area (3)
where F‘" is the force due to the compressed air and F" , F" are due to the resistance of

PE walls as a function of rate and position. According to the assumptions made:

12

F: = Cl" -V,, + 62” -V,3
I",x = Some positive fimction of X

(4 a, b)
where Clv, C2v and Fx are unknown variables of the model. Later some specific forms of
Fx are introduced.

Let me next address Ema), the last calculable unknown needed to complete the
model:

Ff" = —P, - Area (5)
Pu is propagated similarly from its initial value P0. The necessary iterative relationship
similar to (2) is based on the assumption that air compression follows the Ideal Gas Law.

In one time interval dt, air undergoes a transition between the two states:

{pa ,Vf', 7.7"}‘9 {19 V“’ 7:; } For any process involving ideal gas, the relationships

n+l’ Ml’

between absolute temperature, pressure and volume are [9]:

air
P P Tn+l XII

at. = .. ' Tna, . X..+t (6)
Here the assumption of constant crossectional area and thus V‘i'=Area*X is used. Xn
could be replaced by its expression from (2,c) or left like it is noting that it should be
evaluated before PM. is needed. The expression for '1‘" is worked out below and presented
in (12 a).

Using the assumption that the system does not lose energy to the environment, the
first law of thermodynamics applied to the whole cushion can be used to write an energy
balance for every transition from n to the n+1 states:

de’ + de'" + dEfd" = de' + duff“ (7)

Here dUn is a change of internal energy of the air and thermal energy in the walls and de

l3

is the work done by external forces on the air and walls air and walls. dB“ is the
mechanical energy associated with the wall compression (later in the case of elastic
compression it will be equal to the work of Fx and zero for plastic deformation) At the
same time, some thermal energy dQn was exchanged between air and walls. So (7) can be
separated in two equations:

du:ir =dwߢir _dq-Mv
dU,;'°"’ = deW’“ - 4153“" + dQ‘H"

(8 a. b)

The heat term dQn is proportional to the convection heat transfer coefficient, contact
area, process duration and average wall-to-air temperature difference. For small dt, both
'1‘“ and TW" can be considered to be constant and equal to their respective values in the
beginning of the transition.

dQ:—Mv = h . Arr-NV . (flair _ 1:.de ) . dt
All-0w = 6' X‘PP‘N’" .

D

(9 a. b)

Area

Here air to wall contact area A“"“' is estimated from consideration of the PE average cell

dimension, D = 0.00114 m, representing them as cubes or spheres. It stays constant
through out the entire compression in spite of the change in cushion geometry. Estimation
of the convection heat transfer coefficient h represents the most difficulty. One estimate on
it can come from considering the average dimensions of the PE cells and static heat
transfer coefficients for PE and air . Careful consideration of the situation is presented in
[1]. It gives h=44.8W/(m21(). Another estimation is based on experimental data for free
convection in spherical cavities and is also given in [l] but actually comes from [10]. After
consideration of the average cell size and a reasonable range of temperature differences

gives h=21.9W/(m2K). In view of such differences and the fact that heat transfer is likely

14

to be accelerated by the rapid change of the cells geometry, b will be left as an unknown
model parameter.

Using similar considerations of small (it and constancy of pressure and resistance
forces, the work done by the air is dW°"= P*dV and by the walls is dW“"‘ = (F'+F")dx..

din" = I: -(V:"' -V"’)= P, -(X,, - XM)-Area

n+1

(10 a, b)
4W?” = (5' +F..‘)'(X. - Xau)
The changes of internal energy of the air and thermal energy of the walls are proportional

to the changes in their absolute temperatures.

dU:"' = Cv -(T.“" - mil-Area - Xo-p

dUt‘th = PE '(7L"°"'-T.."S“')'m

(11 a, b)
Specific heats Cv and Cpp, could be considered constant for the given temperature and
compression ranges. Their respective values are 717.3 and 1003.8 J/(kg*K) [11].

Equations (8),(9), (10) and (l 1) can be put the together to obtain iterative expressions for

1*" and T"“"‘,

air air Mair _dQ:‘”V
7in] = 1; -
CV -Area- Xo-p
dW'””‘ day”: + def’”

72:1“! = 1;de _ n

( 12 a, b)
C PE - m
It should be mentioned that (12 b ) is worked out for the most general case of some part
of friction work being converted to heat. Later when the nature of friction is discussed,
additional conditions of dlawalls =F"*dx (elastic) and dEW“=o (plastic) will be taken into
consideration.

Now all the formulas can be implemented in the computer program shown in

Appendix.

15

2.4 Study aad verific_ation of the model’s assumptions

2.4.1 Heat transfer considerations

 

Running the model with no Fv or FX and setting the heat transfer coefficient h to
zero gives the adiabatic model described in [l] and shown below by Model 1. Setting h =
1000000 (a reasonable number for infinity in computer applications) gives the isothermal
model. It is shown below by Model 2 in Figure 5. Maximum air and wall temperatures

there are raised only by 2.43 and 1.07 0C above To = 295 K (room temperature) which is

 

 

negligible.
so
— M
40 «- — aha-ll
— ill-am
g 3|! -i
I:
9
‘5 2° -
2
o
0
2 lo -
0 fl: r .L i v _
ii 5 lo rs 20 as an as

 

 

 

-lo
Tlrne [ms]

Figure 5. Adiabatic and Isothermal model checks

Both peak accelerations agree well with the theoretical formulas derived for them in [1].
This fact proves the validity of the small time step approach. However, both models give
an overestimation of peak G when in [l] adiabatic was an over estimation and isothermal

underestimation of cushion curve data. This happens first of all because my real signal is a

16

profile of a first drop on the cushion, not the average of five consecutive drops as it is
done for the cushion curves and secondly it is very likely that the manufacturer put some
safety factor in their curves.

The major differences between the model signal profiles and the real signal do not
leave any doubt that neither isothermal nor adiabatic processes are dominant contributors
in cushion reaction forces, especially during the first half of it. To see whether the
introduction of heat exchange improves the situation, the model was run with h = 22 and

43 W/(mZK), shown in Figure 6 by Model 1 and Model 2 respectively.

 

$888818

Acceleration [6]

~
I. B
L 4
r

 

 

 

 

in

 

Time [ms]

Figure 6. Introduction of heat transfer coefficient h=22 and h=43.

The good thing is that the peak G and the general shape of the signal does not appear to
be very sensitive to the changes in h. A two fold change in h produced only a 5%
difference in peak G and the results are just somewhere between the limiting cases of the

purely adiabatic and isothermal models (not that it was unexpected). But introduction of

heat transfer did not make the model profile more like the real signal. Its peak G’s and
general shape are still wrong.

But heat transfer it is not altogether insignificant. It helps to address the
mechanism of residual cushion compression at the end of the impact and consequently its
recovery times. Up to the moment of maximum compression both the air and walls are
heated up by compression heating and heat transfer. Then on the way up, the air expands
and cools down, but for some time it is still warmer than the walls and continues to beat
them. At some moment later, the air temperatyre catches up with the wall temperature and
the walls start to cool as well. The condition of the impact end P(t,,.d) = P0 is reached well
before the walls have time to reach, but enough for the air to cool below To. In the models
shown above with h = 22 and 43 W/(mzK), the maximum temperature reached by air was
respectively 49 and 31 0C above To. In the end, it cooled down to 31 and 22 0C below To.
The walls reached 0.82 and 0.99 0C at 23.3 and 18.6 ms respectively, which is 82% and
62% of the impact time, after which it had enough time to cool only to 0.58 and 0.43 °C
above To. The fact that the air is colder than the initial conditions means that the cushion
volume and thickness are reduced, referred to as “permanent set” in [5]. Starting from
0.054m thickness the models ended up at 0.0481 and 0.0501. To observe the situation, the
same plots shown above can be integrated to produce acceleration as the function of

cushion thickness.

18

 

8
+

Acceleration [
_
u
I

 

 

 

 

011m 0M 0M 0836 D 0845 0.1150 0.1156

 

Distance [m]

Figure 7. Heat transfer models acceleration vs. X

Here it looks like only Model 1 intersects the X axis at 0.481m exactly right where it is
supposed to, whereas Model 2 is off its 0.0501 m mark mentioned above. This is caused
by the fact that actual axis conversion has been made on the basis of Model 1. Repeating

the same for Model 2 gives the right result too.

2.4.2 Defining the real time of cushion contact

Further investigating the shock profile, one should notice that for the first 1.5 ms
or 0.007m of travel, the real acceleration is considerably lower than predicted by any of
the models. This can not possibly be remedied by the addition of FV or Fx as they would
only increase the model response force. Careful consideration of the possible causes for
the phenomenon leads to the realization that this section is due to the viscous forces
resulting from displacement of thing layer of air by the rapidly oncoming drop head. Both

the drop head and cushion are really experiencing this force. As a result, right before the

19

actual impact, the cushion is slightly presqueezed and the drop head decelerated. In order
to take these affects into account, the real time of impact should be placed at the
intersection of x-axis and a linear extension of the impact leading edge, the deceleration of
drop head calculated and the effective cushion thickness should be measured under slight

pressure.

2.4.3 FV considerations

Now the real cushion response at time of contact is not equal to zero. There is no
air compression and Fx has hardly started yet so the force had to have come from rate
dependent forses. This condition can be used to estimate C1 or C2 or their combination

needed for FV in (4 a). Making such estimations for C1 and taking C2=0 in Model 1.

 

   

 

 

 

4o
36 -P — M
— w
an a» _ r"
y
E15 -
C
0
‘5 an -
E
9
o is -
U
2 Fit Cost: 97.68
m -
5 .-
o ’ . . . ' *‘ .
5 lo l5 20 I? H“ is
.5

 

Time [ms]

Figure 8. Addition of Fv to explain the reaction force at the time of contact.

20

 

Fv pushes Model in the right direction towards better fitting the real signal, shown in

Figure 8. It also explains the negative acceleration on the trailing edge of the real signal

profile. But the effect is not nearly powerful enough. Estimating C2 or the C1,C2

combination gives very similar results. So the contact time force condition is probably not

useful. Letting the computer find the best fit in terms of Cl and C2 gives Figure 9.

 

L I l

u-
u
1

Acceleration [G]
S

u
l
r

—m
—M
-F"

Fit Coef = 19.41

 

 

nr

10

n . a *1 a,

Time [ms]

 

Figure 9. F" adjusted for the best fit

This model behaves much more like the real signal. The fit coefficient went down to 19.4,

meaning that the average error in the predicted and actual accelerations is 4.3 G. The fact

that the leading edge of the model is vertical is discussed later.

2.4.4 Fx and X0 considerations

The next variable which can improve the model is Fx mentioned in (4 b). Since FX

is unknown, the aproach is to first look at what Fx profile is required to make the model

21

behave exactly as the real signal does. Since FX is a function of position, it helps to look at

acceleration vs. position as in Figure 10.

 

  

 

 

 

 

36
m 1
-- M
25 .. — this]
a! -1L - F‘
(E, 15 .
5 lo -l
8
3 5 4
o i i i V1 :_.
arias rum 01135 um was nso miss
.5 a.
-|D ‘-
45

 

Distance [m]

Figure 10. First estimation of possible Fx (Elastic)

I should disregard the leading edge up to 0.045m there Fx is negative as it happens
only because it is compensating for too much Fv there. Otherwise it looks almost linear,
increasing as compression grows . Linear extension of the profile intersects the X axis at
position of contact. This is exactly the behavior one would expect from any elastic
material. So one way to represent Fx is as an elastic force that obeys Hooke’s Law: Fx=
Area*Y*dx/Xo where Y is Young’s modulus for the material: unfortunately it is not
available for the complex structure at hand, but I can let the model find its value for the
best fit. Remember that model should be adjusted to reflect the fact that elastic

compression does not heat the walls: all of the work goes to the mechanical energy

22

(condition dWx = dEx in (12 b)). There dWwalls is not separated in dWx and de but it is
easily done.

Another Fx form which I think could be applicable to such a material as PE
especially at high compression ratios is elastic - perfectly plastic . Its behavior is elastic up
to some critical value of relative compression (strain). Then the force stays constant
(perfectly plastic) all the way until maximum compression and becomes elastic on the way
back. Young’s modulus and thus the stress vs. strain curve slopes stay the same during
expansion as it was during compression. This approach can only improve the model

because it includes the previous one in it. If the model does not prefer the elastic -
perfectly plastic approach, it will find during the search that the length of plastic plateau is
equal to zero thus returning to the elastic representation. This by the way is exactly what
happens.

While changing FV in the process of looking for best fit (and varying KW“.
usefuhless of which will be shown later) I saw that another shape of an arbitrary Fx can

work too, as in Figure 11.

23

 

Acceleration [G]

 

 

 

 

 

Distance [m]

Figure 11. Second estimation of possible Fx (plastic)

This tells me that another form of Fx as a constant (perfectly plastic) should not be

ignored.

Introduction of an elastic FX improved both the general appearance and the fit

coefficient.

24

 

Acceleration [G]
a

    

Fit Coef= 1.85

 

 

 

Time [ms]

Figure 12. Fitted model curve with elastic F"

Fitting all five experimental sets of data and tabulating the results obtained gives Table 1.

Table l. Fitted results for Elastic F" model

 

 

 

 

 

 

 

 

Mass, Drop Height Fit coefficient Clv C2V Young’s modulus
[N*s/m] [N*s2/m2] [kN/mzL
1 17.8 kg 0.4577m 1.9 706 0 79950
2 17.8 kg 0. 2543m 1.1 756 0 1669095
3 17.8 kg 0.661 1m 2.9 600 0 38502.38
4 8.7 kg 0.6611m 6.72 564 0 2245167
5 8.7 kg 0.4577m 9.32 538 0 3100167

 

 

 

 

 

Even the worst fit coefficient of 9.32 looks fairly good as shown in Figure 13.

25

 

 

Acceleration [G]

Fil Cosf= 9 2

 

 

 

 

 

Time [ms]

Figure 13. Worst fit for the model curve with elastic Fx

So apparently the model with an elastic Fx has enough flexibility to be used for filtering.
But the fact that Young’s modulus varies so much means that it is not physical and thus
not acceptable. Trying the elastic-plastic approach gives identical results because the
search algorithm finds best fits when the length of plastic plateau is equal to zero leaving
only elastic part.

Going back to the possibility of perfectly plastic Fx and fitting it to all five
experimental sets of data and tabulating the results gives Table 2.

Table 2. Fitted results for Plastic FK model.

 

 

 

 

 

 

 

 

 

Mass, Drop Height Fit coefficient C1v C2v FX[N]
[N*s/m] [N*s2/m2]
1 17.8 kg 0.4577m 14.56 709 0 892.6472
2 17.8 kg 0.2543m 28.08 780 0 2280.511
3 17.8 kg 0.6611m 6.34 580 0 523.854
4 8.7 kg 0.661 1m 72.6 612 0 869.6859
5 8.7 kg 0.4577m 120 700 0 1043.794

 

 

 

 

26

 

This is not an improvement at all. As a matter of fact, fits to data 2,3 and 5 are clearly not
physical. Their shape is not what one would expect which is reflected in the fit coefficient

too. This is shown below in Figure 14.

 

 

    

Acceleration [G]

Fit Coef= 120.5

 

 

 

 

 

Time [ms]

Figure 14. Non physical result for plastic Fx

But the fact that the two drops with the most energy are still in semi acceptable shape says
that the factor that is not taken into consideration by the model is connected to the energy
of the drops. One of the possible things this factor can be is Xo, the effective thickness of
the cushion. As of now it is taken as the apparent cushion thickness minus the PE
Volume/Area (X0 = 0.054m). Physical observations of cushion compression under static
load shows evidence that the central regions of the cushion experience smaller relative

compression than the ones next to the cushion faces. Moreover, smaller than average size

27

cells if they happen to have thicker than average walls get compressed less or hardly at all.
All this justifies the reduction of cushion thickness in the model, but does not give any
realizable estimates of how much. It essentially introduces another unknown variable in
the model. Xo can not be greater than the apparent thickness and is probably not less than
10-20% of it.

Conducting a search for the best X0 using the elastic and elastic plastic F" does
not lead any where. Fit coefficients there are too small already and the search sets
X0=Xmm to produce the results identical to the ones shown in Table 1. The plastic Fx

approach gives the results in Table 3.

Table 3. Fitted results for Plastic Fx model with adjusted X0

 

 

 

 

 

 

 

Mass, Drop Height Fit ClV c2V Xo/xm FX[N]
coefficient N*s/m] [N*s2/m2]

1 17.8 kg 0.4577m 2.58 625 0 0.7822 668.7

2 mgg 0.2543m 2.53 633 0 0.53 590.2

3 17.8 kg 0.66llm 4.21 545 0 0.88 576.2

4 8.7 kg 0.66llm 13.66 544 0 0.65 483.0

5 8.7 kg 0.4577m 14.42 508 0 0.543 507.8

 

 

 

 

 

 

 

 

Now the fit coefficients are much better and there is some consistency in Fx, ClV and
good logic in X0 (it is decreasing for smaller weights and smaller heights). Fit coefficients
and the general shape of the curves are similar to the ones shown in Figures 12 and 13.
Relative success of the plastic Fx does not mean that elastic can not do better. The
difficulty there is that X0 is acting as a parameter by itself and affects Young’s modulus
too. It so happens that the latter effect is stronger and pulls the model towards Xo=Xmm

and away from a consistent Young’s modulus. A massive search algorithm which takes

28

several sets of data, puts first priority on Young’s modulus consistency and looks for the
minimum square deviation across all the data simultaneously can separate that
dependency. Unfortunately this is not straightforward and has to be left for later. For now
I can insist that the same type of X0 reduction should happen during both model
approaches. I can verify the possibility by forcing Xo to be equal to its values in Table 3.
Running a set of fittings for the elastic Fx again gives the results in Table 4. The numbers
marked by a * in the table were given the value to keep overall consistency of entire table.
They were estimated using an iteration procedure of putting initial values for X0 to get
best Young’s moduli, then using its average to estimate new best values for X0 and over

again.

Table 4. Fitted results for Elastic Fx model with adjusted Xo

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Mass, Fit C1v C2V Xo Young’s
Drop Height coefficient [N*s/m] [N*s2/m2] Xm modulus

1 17.8kg 0.457m 4.28 720 O 0.78 34013.6
2 17.8kg 0.254m 3.49 618 92 0.53 32074.8
3 17.8kg 0.661m 7.78 572 30 0.88 340136"
4 8.7kg 0.661m 15.89 471 37.5 0.65* 35430.8
5 8.ng£.457m 17.94 585 2.87 0.50* 44977.3

 

 

So Young’s modulus can be constant with out introducing too much disturbance in other
parameters. It means that the elastic Fx can be is as good as plastic. Individual differences
between different drops of the same weight from the same height are greater than the
difference between the two models. As a matter of fact, in reality it is most likely that a
combination of the two is working. Knowing the limiting cases is enough to get a good fit

and estimate all the other parameters. The plastic Fx model has slightly better fit

29

coefficients and is much more stable during the automatic searches, so later it is more
convenient to use it for filtering and tabulation of average values for the model

coefficients.

2.4.5 Justification of the model’s use as a filtering tool
The stability of the filtering technique by fitting the model to the real shock pulse is

verified by checking if it can filter out artificial noise in the form of a sine wave.

 

 

 

 

3s
30-
25..
—n.al
920- -Model
C
.9
E is a
.9
§lu -
<
5‘ FilCoef= 383
l s n u 20 filQ/Wllv at
-s

 

Time [ms]

Figure 15. Fit filtering of a noisy signal.

In the example shown in Figure 15, a 1000 Hz sine wave (noise) was added to the
acceleration profile # 1 from Table 3. Running the model’s search algorithm and bringing
the coefficients together in Table 5, we can see that the results of the fit filtering are

practically identical.

30

Table 5. Results of fitting original and noisy signals

 

 

 

 

 

 

 

 

 

 

Signal fitted Fit coefficient Clv Xo/Xmm FX[N]
[N*s/m]

1 Original 2.58 625 0.7822 668.7

2 Noisy 4.64 620 0.7834 671

 

The fit coefficient was expectedly higher, but over all, the model was not confused by the
additional noise. Increasing the noise frequency gives even better results. Highter or lower
noise amplitudes appears not to have any affect at all. Frequencies below 500 Hz should
be avoided, but those usually have low amplitudes and do not matter any way. So the
model appears not to be sensitive to sine wave noise too much and probably it is as

effective in filtering out natural noise in the signal as well.

2.4.6 Verifigtion of suggested cushion recovery tm

There is still no real handle on the heat transfer coefficient h because the difference
in fit coefficients for h=22 and h=43 with all the other coefficients kept the same is on the
order of 5%, which is easily absorbed in the uncertainties of all other coefficients. Heat
transfer coefficients taken outside the 5 to 50 range consistently give worse fit
coeflicients which cannot be compensated by any other parameters. For the absence of
any better choice, I will take h=32 for the fit filtering and results tabulation. Now for this
fixed value of h the model gives final temperatures for the air and walls for example for
the #1 drop in Table 3 of 0.56 °C and -23.24 °C above initial temperature respectively. It
will take heat transfer and work of air expansion against Po only 18.7 ms to warm up the
air and cool the walls to reduce the difference by a factor of ten. Unsteady state heat

transfer is dominating the work factor in the energy loss process and it is an exponential

31

process, so it its not surprising that the next factor of ten reduction is reached at 37.5ms
(roughly double the time). Thus by the cushion recovery time recommended in [5] of one
minute between the consecutive test drops is over the air-to-wall the temperature
difference would be reduced by about 3000 orders of magnitude (overkill). On the other
hand, both temperatures will converge on the value which is 0.5 °C higher than the initial
temperature. The time to dissipate this energy can be estimated. Because of poor heat
conductivity of the cushioning material, it will clearly be more than a few minutes. The
energy of consecutive drops will accumulate even if the one minute recommendation is
followed. After a dozen or so drops, the elevated cushion temperature could be felt by
hand. Fortunately, raising the initial temperature by even 10°C does not introduce too
much difference. For the above mentioned example, peak G would increase by 0.5 G

which is only 1.5 %.

2.4.7 Consideration of model’s vertical leading edge
The fact that model has vertical leading edge and real signals never do can be
explained by several factors:
0 Cushions can not possibly have perfectly uniform thickness so they meet the drop head
with the high point first and then rapidly increase the area of contact.
0 The top layer of cells is no doubt damaged more than the center so the FY coefficient is
smaller while these are compressed and is rapidly reaching a constant value. This
explanation has a supporting evidence of slight and not explained by anything else

increase in cushion response right before the end of impact.

32

0 Some of the restricting forse comes fi'om the slider rods. It is practically zero when the
head is freely sliding. As the head comes into contact with the cushion, the reaction
force is likely not to be passing exactly through its center of mass, gradually creating a
rotational moment which jams the rods and increases their friction.

No matter which factor or what combination of them is really working, this can be
taken into affect by adding just one more variable. My estimations show that it will
drastically improve the fit coefficient but does not change any other important parameters

by more than l-2%, so I will not complicate the model.

33

CHAPTER 3
RESULTS AND CONCLUSIONS
3.1 Corrections and additions to th_e model assumptions.

As a result of the model study, the following corrections and additions to the

original model assumptions have been made:

Xo -effective cushion thickness is introduced as a model unknown.

C2V Should be completely dropped from the model since it is consistently found that
best fits have it equal to zero indicating no Fv dependence on second order of
compression rate. This does not contradict any physics of the system because there is
no considerable air flow with which the second order of rate resistance is usually
associated.

Plastic Fx is preferable to elastic not only because of the ease of use but because of
better fits, stability of the model coefficients, and the plastic nature of PE. The
condition dEn=0 is used in (12 b). It should be mentioned that the latter condition
practically affects only the temperature of the PE walls and not that much either since
FV > Fx for the most of time.

The heat transfer coefficient is taken to be 32 W/(mzK), the middle of its reasonable
range for the lack of better choice. Fortunately, because of the model’s extreme
insensitivity to h, this does not render the model inaccurate.

The vertical leading edge of the modeled profile is ignored for the lack of verifiable

physical explanation and its little overall effect on the model.

3.2 Model use for tabulating the material promrties

 

Now with the model studied and complete, given a new cushion it is enough to
make several drops to get averages on Fx, X0, and Clv. Within reasonable limits, they
could be interpolated for any other real situation. Using the interpolated values and
running the model we can receive an estimated acceleration profile and predict the peak G,
impact duration, restitution coefficient, maximum compression as well as almost anything
else one can think of. For example, using Table 3 we can interpolate model coefficients for

a {17.8 kg 0.557m} and {13.25 kg 0.6611m} drops. They are given in Table 6.

Table 6. Interpolated parameters for an arbitrary drops.

 

 

 

 

 

 

 

 

 

Mass, Drop Height Cl" Xo/Xmm F"[N]
[N*s/m]

1 17.8 kg 0.557m 585 0.83 615

2 13.25 kg 0.6611m 545 0.765 529

 

Predicted acceleration profiles together with some other useful information are shown in

Figure 16, Models 1 and 2 respectively.

35

 

- “I! — “12
50 Gmalr= 34.8415 45.8297
t(Gm)= 12.2 10.5
40 .. Xmin= 0.02157 0.01%
Xfinal = 0.0419 0.0384
trro=rw) =17.9 '55
so -. t(V=0)=14.0 12.0
Restitytion 0.545% 051725

Acceleration [G]
B

-

a
l
t

 

a
.1

Air PE Walls Air
- all=32.473 0.8525 34.8879 0.
- fin =-20.211 0.5877 -21.0702 0.

 

PE Wall

787
2

 

 

 

 

Time [ms]

Figure 16. Model Predictions for {17.8 kg 0.557m} and [13.25 kg 0.6611m} drops.

3.3 Model use for fit filtering

The model can be used to fit itself to the real signal profile essentially filtering it

and finding the coefficients for this particular drop. The only input parameters required by

the model are drop height, drop head mass and cushion thickness which are usually

known. The result of such an approach is shown in Figure 17.

36

 

813888

 

 

 

 

E w 1'

c

9 ~ /

g .

3 15 «-

U

3 lo ..
5 1’ j Fit Coer= 2.21
o l i :

s W . 30 35

6 \

 

Time [ms]

Figure 17. Unknown signal fit filtering

The real signal in Figure 17 is the result of {25.6kg 0.508m] drop on the same cushion but
the which was beaten more than cushions used for the data in Table 3. So the model not
only gives a good fit but predicts a decrease in Fx and Clv, an increase in X0 , peak G,
restitution coefficient and other logical changes in temperatures reached due to the fact

that the cushion is worn out.

CONCLUSIONS

The model of the cushion compression process in conjunction with the best fit
search algorithm allows the user to process experimental data on cushion tests more
accurately and efficiently than using existing techniques.

Because fit filtering as opposed to Fourier filtering is based on the real physics of
the process, it can recognize and partially cancel some experimental systematic errors. It
also deals with random errors (noise) at least as well as the Fourier filtering, but does not
cause a dilemma of filtering frequency choice. The process of filtering simultaneously
provides information on all useful drop parameters.

The type of material constants generated as the result of such data analysis carries
more information about the cushion properties and is much more compact than the
standard cushion curves. A new ASTM standard for generation and application of these
properties can be easily designed. It will require the use of a computer and a standard
algorithm, but in our time it is fairly easy to set up. The advantage of using the computer
is that the search for the material capable of providing appropriate product protection
levels at the minimum cost can be automated.

Extension of the model search capabilities and abilities to work with multiple sets
of experimental data is not very difficult to implement. It most likely will lead to better
approximation of the heat transfer coefficient and incorporation of a different apparent
cushion thickness dependence, which in turn will make the model more accurate and

further simplify the material properties tabulation.

38

APPENDIX

39

Computer Program Code

Public Sub Build_model(lndex As Integer)
’Index gives ability to use the same function for different variable sets
col = (Index + 3) * 2
Sum = 0
VnotDone = True
TnotDone = True
Temp_FofX = 0
dt = Tstep / Niterations
XN = Xo(Index)
XaN = XN
VN = -Sqr(2 * g * DropHeight)
AN = 0
PN = P0
TN_Air = InitialTemp + 273.15
TN_wall = TN_Air
CurTime = TContact(Index)
M_air = 1.225 * Xo(Index) * Area ’For my volume
N = M_air/ 28.9644 ’number of moles
Cp = 1003.8 "‘ M_air ’for my mass
Cv = 717.29875 * M_air 'For my volume
Cpe = 1901 * Mpe 'For my mass
h = Ktherm(Index) * 5263 * Area * X_apparent
Convective Heat Transfer coefficient for my area of contact
With FFiltering.Graphl .DataGrid
ContactTimeInteger = IntRangeStart + CInt('IContact(Index) / Tstep) + 1
Will start impact there
If TContact(Index) > 0 Then
Aend = 0
For i = 1 To ContactTimeInteger ’ no cishion compression yet

.GetData i, 2, temp, 0

.SetData i, col, 0, 0

.SetData i, 10, 0, 0

Abeg = Aend

Aend = temp * g

For j = 1 To Niterations

'before the contact Xn=Xo

VN = VN + (Abeg + (Aend - Abeg) * j / Niterations) * dt
’Velocity of the mass
Next j
Next

AN = Aend
End If
VContact(Index) = VN

counter = IntRangeStart + CInt(T_modelStart(Index) / Tstep) + 1
’Will start impact there
If (counter > IntRangeStart + 1) And (counter > ContactTirnelnteger) Then
Temp_FofV = FofV (Index, VN)
Temp_FofX = AN * mass - Temp_FofV - Area * (PN - Po)
Abeg = AN
For i = ContactTirnelnteger To counter
If FFiltering.View_FofV(Index).Checked Then .SetData i, 10, FofV (Index, VN)/
g / mass, 0
If FFiltering.View_FofX(Index).Checked Then .SetData i, 10, (AN "' mass -
Temp_FofV - Area * (PN - Po)) / g / mass, 0
.SetData i, col, 0, 1
.GetData i, 2, temp, 0
Abeg = Aend
Aend = temp * g
For j = 1 To Niterations
step (Index)
AN = Abeg + (Aend - Abeg) * j I Niterations
Next j
Next
End If
CurTime = T_modelStart(Index)

Temp_FofV = FofV (Index, VN) ’ It is used later to show the graphs
Temp_FofX = FofX(Index, XaN, VN, FofX_type_Number(Index))

VN_modelStart(Index) = VN
AN_modelStart(Index) = AN
Temp_FofV _mode18tart = Temp_FofV
Temp_Fofl(_modelStart = Temp_FofX

Do While (PN >= Po) And (counter <= IntRangeEnd)
.SetData counter, col, AN / g, 0 Recording next point
.GetData counter, 2, temp, 0
Sum=Sum+(temp-AN/g)"2
If FFiltering.View_FofV(Index).Checked Then .SetData counter, 10, Temp_FofV / g
/ mass, 0
If FFiltering.View_FofX(Index).Checked Then .SetData counter, 10, Temp_FofX I g
/ mass, 0
For i = 1 To Niterations
step (Index)
AN = (Area * (PN - Po) + Temp_FoD( + Temp_FofV) I mass
Next 1

41

counter = counter + 1
Loop

For i = counter To IntRangeEnd

.SetData i, col, 0, 1 ’clearing the rest of the points

.SetData i, 10, 0, 1
Next i
FFiltering.XfBox = Format(XN, "#00W”)
FFiltering.TfBox = Format(TN_Air - 273.15 - InitialTemp, "HOW?
FFiltering.TmaxBox = Format('T_Air_Max - 273.15 - InitialTemp, "#0.0###")
T_Air_Max = 0
FFiltering.TwallmaxBox = Forrnat(T _wa11_Max - 273.15 - InitialTemp, "##0.0W")
T_wall_Max = 0
FFiltering.TwallFinalBox = Format(TN_wall - 273.15 - InitialTemp, "#0.0M")
FFiltering.GmaxBox = Format(Amax / g, "##0.0HF#")
Amax = 0
FFiltering.FitCoefBox = Format(Sum/ (counter - (IntRangeStart +

CInt(T_modelStart(Index) / Tstep) + 1)), ”W0.0#")

FFiltering.RestitytionBox = Format(-VN / VContact(Index), "0%”

End With
End Sub

Public Sub step(Index As Integer) ‘ Is called from the main body of Build_model

dX=VN*dt+AN*dt*dt/2

XaN = XN + dX

VN = VN + AN * dt

IfVN > 0 And VnotDone Then
Xmin = XN
FFiltering.XminBox = Format(XN, ”0.0#ii##")
FFiltering.TimeOfVisOBox = Format(CurTime * 1000, "##0.0HH#")
VnotDone = False

End If

'work done by gas
dW_air = PN * Area * dX 'exact adiabatic PN * Area * XN / (gam - l) * (l -
(XN / XaN) " (gam -1))

'work by forses of friction done on the walls

dW_Friction = -dX * (T emp_FofV + Temp_FofX) ' no Temp_FofX if it was a
spring

' heat given to the air

dQ = h * (TN_wall - TN_Air) * dt

42

TaN_Air = TN_Air + (dQ - dW_air) / Cv TN_Air * (XN / XaN) " (gam - l)
TN_wall = TN_wall + (dW_Friction - dQ) / Cpe

PN = PN * (XN I XaN) * (TaN_Air/ TN_Air)

Temp_FofV = FofV (Index, VN) ’ It is used later to show the graphs
Temp_FofX = FoD((Index, XaN, VN, FofX_type_Number(Index))
CurTime = CurTime + (it

XN = XaN

TN_Air = TaN_Air

If TN_Air < TN_wall And TnotDone Then
FFiltering.TimeOf’TisTwallBox = Forrnat(CurTime * 1000, "#00W")
TnotDone = False

End If

If AN > Amax Then
Amax = AN
VofAmax = VN
XofAmax = XN
FFiltering.TimeOfArnaxBox = Forrnat(CurTime * 1000, ”##0.0##")

End If

If TN_Air > T_Air_Max Then
T_Air_Max = TN_Air

End If

If TN_wall > T_wall_Max Then
T_wall_Max = TN_wall

End If

End Sub

Private Sub Minimize_Click(Index As Integer)’ used for variable search and fit
minimization.

Dim Smalest_variable As Double

Dim VarStep As Double

Dim Error As Double

Dim Smalest_Error As Double

Dim N_Intervals As Integer

Dim Error_Valyes(0 To 20) As Double

Dim max, min As Double

Dim Best_Smalest_Error As Double

Dim best_Xo As Double

Dim best_FofV_Cl As Double

Dim best_FofV_C2 As Double

Dim best_I(therm As Double

Dim best_FofX_Va1yes(0 To 4) As Double
Dim best_FofX_PosFrac(0 To 4) As Double

43

 

If T_modelStart(Index) = TContact(Index) Then

Randomize ’ Initialize random-number generator.
CoeficientsOK_Click (Index) ’ get the numbers in and run the model first time to get the
constants
Best_Smalest_Error = 1E+300 ’ get first best error to be big
For k = 1 To N_RandomStarts ’N of random (10 to 20) numbers of intervals applied.
N_Intervals = Int((N_Intervals_above_5 * Rnd) + 5) ’ Generate random value
between 10 and 20.
For 1 = 1 To N_variable_sycles ’N times to cycle all variables

’Xo minimization

If XoAsFractionBox(Index).BackColor = 12632256 Then ’variable is marked for
minimization
Smalest_Error = 1E+300 ’ get first refference error to be big
max = Xo_max(Index)
min = Xo_min(Index)
For m = 1 To N_Step_Reductions 'N times to reduse the step and interval
VarStep = (max - min) / N_Intervals
For p = 0 To N_Intervals
Xo(Index) = min + p * VarStep

Error = Quick_model_Sum(Index)
If Error < Smalest_Error Then
min_number = p
Smalest_Error = Error
Smalest_variable = Xo(Index)
End If
Next p
If min_number > 0 Then min = min + VarStep * (min_number - 1)
If min_number < N_Intervals Then max = min + VarStep * 2
Next m
Xo(Index) = Smalest_variab1e
End If

’Same thing for FofV_Cl
IfF1(Index).BackColor = 12632256 Then Variable is marked

Smalest_Error = 1E+300 ’ get first refference error to be big

max = FofV_Cl_max(Index)

min = FofV_Cl_min(Index)

For m = 1 To N_Step_Reductions ’N times to reduse the step and interval
VarStep = (max - min) / N_Intervals
For p = 0 To N_Intervals

FofV_C1(Index) = min + p * VarStep
Error = Quick_model_Sum(Index)
If Error < Smalest_Error Then
min_number = p
Smalest_Error = Error
Smalest_variable = FofV__Cl(Index)
End If
Next p
If min_number > 0 Then min = min + VarStep * (min_number - 1)
If min_number < N_Intervals Then max = min + VarStep * 2
Next m
FofV_C1(Index) = Smalest_variable
End If

’Same thing for FofV_C2
If F2(Index).BackColor = 12632256 Then ’variable is marked
Smalest_Error = 1E+300 ’ get first refference error to be big
max = FofV_C2_max(Index)
min = FofV_C2_min(Index)
For m = 1 To N_Step_Reductions ’N times to reduse the step and interval
VarStep = (max - min) / N_Intervals
For p = 0 To N_Intervals
FofV_C2(Index) = min + p * VarStep
Error = Quick_model_Sum(Index)
If Error < Smalest_Error Then
min_number = p
Smalest_Error = Error
Smalest_variab1e = FofV_C2(Index)
End If
Next p
If min_number > 0 Then min = min + VarStep * (min_number - 1)
If min_number < N_Intervals Then max = min + VarStep * 2
Next m
FotV_C2(Index) = Smalest_variable
End If

’Same thing for Ktherm

If K_therm(Index).BackColor = 12632256 Then ’variable is marked
Smalest_Error = 1E+300 ’ get first refference error to be big
max = h_max(1ndex)
min = h_min(Index)
For m = 1 To N_Step_Reductions ’N times to reduse the step and interval

45

VarStep = (max - min) / N_Intervals
For p = 0 To N_Intervals
Ktherm(Index) = min + p * VarStep
Error = Quick_model_Sum(Index)
If Error < Smalest_Error Then
min_number = p
Smalest_Error = Error
Smalest_variable = Ktherm(Index)
End If
Next p
If min_number > 0 Then min = min + VarStep * (min_number - 1)
If min_number < N_Intervals Then max = min + VarStep * 2
Next m
Ktherm(Index) = Smalest_variable
End If

’Same thing for 5 FofX_Valyes and FofX_PosFrac
For r = 0 To 4
If BValuesBox(r + 5 "' Index).BackColor = 12632256 Then ’variable is marked
Smalest_Error = 1E+300 ’ get first refference error to be big
max = FofX_Valyes__max(Index, r)
min = FofX_Valyes_min(Index, r)
For m = 1 To N_Step_Reductions ’N times to reduse the step and interval
VarStep = (max - min) I N_Intervals
For p = 0 To N_Intervals
FofX_Valyes(Index, r) = min + p * VarStep
Error = Quick_model_Sum(Index)
If Error < Smalest_Error Then
min_number = p
Smalest_Error = Error
Smalest_variable = FofX_Valyes(Index, r)
End If
Next p
If min_number > 0 Then min = min + VarStep * (min_number - 1)
If min_number < N_Intervals Then max = min + VarStep "' 2
Next m
FofX_Valyes(Index, r) = Smalest_variable
End If

If BPositionBox(r + 5 "' Index).BackColor = 12632256 Then ’variable is marked
Smalest_Error = 1E+300 ’ get first refference error to be big
max = FofX_PosFrac_max(Index, r)
min = FofX_PosFrac_min(Index, r)
For m = 1 To N_Step_Reductions ’N times to reduse the step and interval

 

VarStep = (max - min) / N_Intervals
For p = 0 To N_Intervals
FofX_PosFracandex, r) = min + p * VarStep
Error = Quick_model_Sum(Index)
If Error < Smalest_Error Then
min_number = p
Smalest_Error = Error
Smalest_variable = FofX_PosFracandex, r)
End If
Next p
If min_number > 0 Then min = min + VarStep * (min_number - 1)
If min_number < N_Intervals Then max = min + VarStep * 2
Next m
FofX_PosFracandex, r) = Smalest_variable
End If
Next r
Next 1
If Best_Smalest_Error > Smalest_Error Then
best_Xo = Xo(Index)
best_FofV_C l = FofV_C 1 (Index)
best_FofV_C2 = FofV_C2(Index)
best_Ktherm = Ktherm(Index)
For r = 0 To 4
best_FofX_Valyes(r) = FofX_Valyes(Index, r)
best_FofX_PosFrac(r) = FofX_PosFrac(Index, r)
Next r
End If
Next k
Xo(Index) = best_Xo
FofV_C 1 (Index) = best_FofV_C1
FofV_C2(Index) = best_FofV_C2
Ktherm(Index) = best_Ktherm
For r = 0 To 4
FofX_Valyes(Index, r) = best_FofX_Valyes(r)
FofX_PosFrac(Index, r) = best_FofX__PosFrac(r)
Next r
Write_To_Boxes
Build_model (Index) ’run the model to update graph
Else
MsgBox ("Time of modelstart should be the same as time of contact to run
minimization")
End If

End Sub

47

BIBLIOGRAPHY

48

. Burgess, G.J. “ Some Thermodynamic Observations on the Mechanical Properties of
Cushions”, Journal of Cellular Plastics, 24, 57-69 (Jan. 1988)
. Burgess, G.J. “ Consolidation of Cushion Curves” , Packaging Technology and
Science, 17, 189-194 (1990)
. Hilyard, M.C. Mechanics of Cellular Plastics, New York: Applied Science Publishers
(1982)
. Andrew W. Chen, “ Finite Difference and Finite Element Modeling of Closed-cell
Cushions With Block and Ribbed Geometry Based on Strain Dependent Terms”, MS
Thesis, School of Packaging, Michigan State University (1991)
. ASTM, Selected ASTM Standards on Packaging, 4m edition, ASTM D1596-
91,American Society for Testing and Materials, Philadelphia, PA, (1994)
. Handbook of digital signal processing, Academic Press, San Diego,(1987).
. Throne, J. and R. Prodelhof. “Closed Cell Foam Behavior under Dynamic Loading-I.
Stress-Strain Behavior of Low Density Foams.” Journal of Cellular Plastics,
(NovJDec. 1984)
. Throne, J. and R. Prodelhof. “Closed Cell Foam Behavior under Dynamic Loading-II.
Dynamics of Low Density Foams.” Journal of Cellular Plastics, (Jan/Feb. 1985)
. Zemansky, M. and R. Dittrnan. Heat and Thermodynamics, McGraw-Hill, New York,
(1981)

Kreith, F. Principles of Heat Transfer, Intex Educational Publishers, New York,

(1973)

49

 

“111111111111111: