Study of Hakiel's Nonlinear Model for Wound Roll Stress Implementations

Comparing Mathematica, Matlab, and COMSOL Script

Alexei V Azarov and Chien Liu

Imaging Science Lab
ZINK Imaging, Inc.

Background

Zink Imaging www.zink.com

Zink == "Zero Ink" 
Printing full color pictures without ink cartridges or ribbons

The Magical Fax Paper
Three dye layers activated independently by single pass thermal printing

The Problem in Winding Paper

Long sheets of paper ("web") develop undesirable wavy patterns
when wound into a roll

Nomenclature

WEB ≠ Spider Web
WEB == Long Sheet of Paper

The Problem in Winding Paper

Two (Big) Steps to Model The Paper Winding

Calculate Equilibrium Stress Distribution (This Talk)

Z. Hakiel, Nonlinear model for wound roll stresses, Tappi Journal, May 1987, pp .113 - 117

Model Assumptions

1. Roll approximated as concentric loops of web (no slip assumption)

2. Web assumed perfectly uniform in thickness and width (flatness assumption)

3. The stresses and strains assumed to be functions of radius only and not of axial or circumferential position (symmetry assumption)

4. The stresses in axial direction are equal to zero (2D plane stress assumption)

Use Perturbation Method to Check for Instability (Not Today)

Equations for Equilibrium Stress Distribution

Winding one layer at a time: at each step inter-layer pressure changes by δP(r), r ≡ radius

(r^2)(((d^2)δP)/(d(r^2)))+3 r ((dδP)/(dr))-((g^2)-1)δP=0, (g^2)≡((Subscript[E, φ])/(Subscript[E, r]))

with two boundary conditions:

(Subscript[r, c]) ((dδP)/(dr))=(((Subscript[E, φ])/(Subscript[E, c]))-1+ν) δP          at r=(Subscript[r, c]) (radius of the core)

δP=(h/(Subscript[r, 0]))(Subscript[T, w])          at r=(Subscript[r, 0]) (radius of current outmost layer)

where (Subscript[E, φ]), (Subscript[E, r]), and (Subscript[E, c]) are the elastic moduli of the web (circumferential), roll (radial), and core respectively, h ≡ layer thickness, (Subscript[T, w]) ≡ web winding tension. In general (Subscript[E, r]) (and thus (g^2)) is a function of the current pressure distribution P(r) because (Subscript[E, r]) depends on the thin air layer between wraps.

After solving for the incremental pressure δP, the incremental in-roll tension δT(r) can be obtained:

δT=-δP-r ((dδP)/(dr))

Analytical Solution in Special Case

To compare candidate software and implementation we have studied a special case with analytical solution

If (Subscript[E, r]) (and thus (g^2)) is a constant then the solution for δP can be found with DSolve function:

equ=dp''[r]+3*dp'[r]/r-(γ^2) dp[r]/(r^2)==0;
bcc=dp'[1]==-(α^2) dp[1];
bco=dp[r0]==(β^2);

Sol1=DSolve[{equ,bcc,bco},dp[r],r];

Sol1=Simplify[Sol1,Assumptions->(Im[γ]==0)&&(Im[α]==0)&&(Im[β]==0)]

where the radius r is normalized to the core radius (Subscript[r, c]), (β^2)≡(h/r0)(Subscript[T, w]), r0 is the current outer radius and h is the increment of r0.

ddp1[r_,r0_,α_,β_,γ_]=dp[r]/.Sol1[[1]];

If h approaches to 0 then the pressure distribution is

P(r)=(∫_r^(Subscript[r, out]))ddp1(r,(Subscript[r, 0]),α,β,γ)d(Subscript[r, 0])

The last integral can be calculated in Mathematica and expressed with the Hypergeometric2F1 built-in function. For a particular set of α, β and γ values the expression looks like this:

Integrate[ddp1[r,r0,α,β,γ],{r0,r,4},Assumptions->(Im[r]==0)&&(r≥1)&&(r<4)]

(1/(-29+79 (Sqrt[6])))165 (2^(1-2 (Sqrt[6]))) (r^(-1-2 (Sqrt[6]))) (89-10 (Sqrt[6])+(89+10 (Sqrt[6])) (r^(2 (Sqrt[6])))) (-4 (r^(Sqrt[6])) Hypergeometric2F1[1,(1/12) (6-(Sqrt[6])),(3/2)-(1/(2 (Sqrt[6]))),(((2^(-4 (Sqrt[6]))) (-8521+1780 (Sqrt[6])))/7321)]+(4^(Sqrt[6])) r Hypergeometric2F1[1,(1/12) (6-(Sqrt[6])),(3/2)-(1/(2 (Sqrt[6]))),(((-8521+1780 (Sqrt[6])) (r^(-2 (Sqrt[6]))))/7321)])

Finally the tension distribution can be found as follows:

T=-P-r ((dP)/(dr))

Computation Algorithm

Starting from the 1st layer:
1. Find δP and δT from the equations above
2. Find the new pressure distribution as P=P+δP
3. Find the new tension distribution as T=T+δT
4. Calculate (Subscript[E, r]) and (g^2)
5. Increment r to next layer
6. If all layers wound than exit calculation loop, if not go to 1.

P and T are obtained as sums of δP and δT => computational errors are accumulative

Need to solve boundary problem for δP for each layer. Average roll consist of several thousands of layers => calculation time is important

Algorithm Implementations

The numerical algorithm was implemented using Mathematica, Matlab, and COMSOL in several ways:

Method 1. Matlab and COMSOL script. The finite difference method was used to solve the boundary problem at each step. The tri-diagonal matrix was inverted using custom-built forward and backward substitution formulas.

Method 2. Matlab and COMSOL script. The same implementation as in Method 1. but the tri-diagonal matrix was inverted using the "\" operator.

Method 3. COMSOL Script 1.1 script. The COMSOL 1-d general form stationary application mode of was used to solve the boundary problem on each step.

Method 4. Mathematica Module. The Mathematica built-in boundary solver was used to solve the boundary problem at each step. The Mathematica program may not be fully optimized.

Parameter Values

Typical values were used in the numerical computation

tn=330;(* tension at the outer surface 330 psi *)
mr=100000;(* elastic modulus radial direction 100000 psi *)
mt=600000;(* tangential direction *)
mc=890000;(* core *)
rc=6;(* radius of core 6 inches *)
th=0.003;(* thickness of web 0.003 inches *)
Nw=6000;(* number of layers 6000;*)

Corresponding values used in the analytical solution

γ=(Sqrt[(mt/mr)-1]);
α=(Sqrt[1-(mt/mc)]);
β=(Sqrt[tn]);

Timing Result

The test calculations were performed on an example of a roll with 6000 layers. 

A laptop computer with Intel Centrino Duo 2.0GHz and 2GB memory was used for test calculations.

Method 1. 4sec with Matlab. 4min 35sec with COMSOL Script.

Method 2. 46min with Matlab. ~6h with COMSOL Script.

Method 3. ~8h with COMSOL Script.

Method 4. 20min 34sec in Mathematica 6.0

Pressure Distribution Result

Tension Distribution Result

Negative tension may lead to onset of instability

Small difference in tension may make a difference between stable and unstable equilibrium

Conclusions

Mathematica is the best tool to tackle problems starting from analytical formulation to numerical computation

The Mathematica 6.0 nonlinear boundary problem solver is slower than the direct finite difference method custom-built with the Matlab script

Nevertheless it holds the first place among other built-in solvers such as the Matlab "\" operator

The accuracy of the Mathematica solver is about the same as the accuracy of the direct finite difference method