Part A of this work discussed in detail a rheological constitutive model to describe the stress and strain behavior of time dependent materials (1). Part B explains the results of some example problems solved using the proposed rheological model and curvilinear form of the equation of virtual work. The problems considered are a polymeric beam modeled as a:
1. Generalized Maxwell Fluid subject to only its own weight.
2. Generalized Maxwell Fluid subject to a surface traction load.
3. Generalized Maxwell Solid subject to a surface traction load.
These problems have many applications in the plastics industry, particularly in the thermoforming process. In thermoforming, a plastic sheet is clamped in a frame, heated until softened, and then shaped into a mold to form a plastic part via vacuum and/or air pressure. In many applications, the sheet is pre-stretched by inflation into a bubble before it is formed into a part. The thermoforming process takes place in the rubbery region of the material above the glass transition temperature ([T.sub.g]) and below the melt temperature ([T.sub.m]), where the material has enough hot strength capability to support itself when heated. Parts typically thermoformed are cups, food containers, trays, packaging containers, etc. In addition, spas, tubs, shower units, electronic part and machinery covers, etc., as well as certain automotive items are thermoformed.
The material studied was a commercial sheet stock manufactured by the Firestone Synthetic Rubber & Latex Company of Akron, Ohio. The material consists of a blend containing 50% of a multiblock styrene-butadiene copolymer (Stereon 900) and 50% of a polystyrene (BASF 2524). The copolymer contains 73 wt% bound styrene, and has a weight average molecular weight ([M.sub.w]) of 141.000 and a number average molecular weight ([M.sub.n]) of 122,000. The copolymer imparts flexibility to the material which is used to make drinking cups. packaging trays, containers, dishes, and other products. The glass transition temperature ([T.sub.g]) of the polystyrene is 112.9[degrees]C. The blend begins to soften at this temperature. The thermoforming temperature range for this material is approximately 130[degrees]C to 150[degrees]C.
GENERALIZED MAXWELL FLUID AND SOLID MODELING
A time dependent constitutive model requires: 1. an elastic component so that the integrity of the sheet (object) is retained; 2. a viscous component to enable the material to flow and deform to the larger depths. In addition, thermal expansion/ relaxation and contraction effects and the temperature dependence of the material parameters must be considered in the context of a large deformation theory.
A multi-element generalized Maxwell fluid and solid describes the response of the material. The "method of reduced variables" was used to combine the data at several different temperatures into a single master curve (2-4). The Levenberg-Marquardt method (5) was used to fit the master curve to the experimental data. The styrenic material studied here was also modeled previously using a one-dimensional formulation with only the body force present (6).
Dynamic mechanical testing was performed to determine the mechanical response of the material. A parallel-plate viscometer (Rheometrics Mechanical Spectrometer) with a radius of 12.5 mm and a gap of 1.4 mm was used to measure the storage G' and loss G" shear moduli as a function of frequency [omega] at different temperatures (140[degrees]C, 150[degrees]C, and 160[degrees]C).
The data at the three temperatures were combined into a master curve using the formulas (2)
[G'.sub.r]([omega], [T.sub.o]) = G'([omega], T)[T.sub.o]/T (1a)
[G".sub.r]([omega], [T.sub.o]) = G" ([omega], T)[T.sub.o]/T (1b)
The temperatures in Eqs 1 are absolute temperatures. The reference temperature ([T.sub.o]) was taken as the middle temperature. The variables [G'.sub.r] and [G".sub.r] represent the reduced storage and loss shear moduli, respectively.
The moduli are shifted horizontally by log([a.sub.T]) where [a.sub.T] is the shift factor that is given by (2)
[a.sub.T] = [[eta].sub.o] (T) [T.sub.o]/[[eta].sub.o] ([T.sub.o])T (2)
In the above expression, [[eta].sub.o] is the zero shear rate viscosity. The numerator is related to the time required for a particular response at a temperature T, while the denominator is related to the time required for the same response at the reference temperature [T.sub.o]. The moduli values for 140[degrees]C and 160[degrees]C were shifted, using 150[degrees]C as the reference temperature. These reduced values were then used to compute a temperature dependent generalized Maxwell model.
The zero shear rate viscosity could not be experimentally obtained. These values were derived by formulating the storage and loss moduli curves for each temperature and then calculating the zero shear rate viscosity from the relationship (2):
[[eta].sub.o] = [lim.sub.[omega][right arrow]0] G"/[omega] = [summation over (N/p=1)] [G.sub.p] [[xi].sub.p] (3)
for the case of the Maxwell Fluid. In Eq 3, N is the number of Maxwell elements, while [G.sub.p] and [[xi].sub.p] are the elemental shear modulus and relaxation times, respectively. The variable N is replaced by N-1 for the case of a Maxwell solid. The results were then related to the absolute temperature by means of an Arrhenius dependence [[eta].sub.o] = Aexp ([E.sub.act]/(RT)), where A is a pre-exponential factor, [E.sub.act] is an energy of activation constant, R is the universal gas constant, and T is again the absolute temperature. These fitted zero shear rate viscosities were then used to obtain the shift factors.
The values for the shear modulus [G.sub.p] and the relaxation time [[xi].sub.p] for each element are obtained from the relationships (2):
[G.sub.p](T) = [G.sub.[r.sub.p]]([T.sub.o]) (T/[T.sub.o]) (4a)
[[xi].sub.p] (T) = [[xi].sub.[r.sub.p]]([T.sub.o]) [a.sub.T] (4b)
The shift factor was calculated to be 3.0275 and 0.3475 for 140[degrees]C and 160[degrees]C, respectively, for the case of both the Maxwell fluid and solid. The parameter values for [G.sub.p] and [[xi].sub.p] at 150[degrees]C are shown in Tables 1 and 2. A comparison of the reduced experimental moduli with the modeled values is shown in Fig. 1a and b.
The styrenic sheets studied contained residual stresses due to prior extrusion processing operations. Upon heating, the sheets were observed to distort as the residual stresses were removed in the absence of any loads. Thermal relaxation measurements were made by cutting 25.4 mm (one inch) squares and heating the material in the oven to the proper temperature. The samples were left in the oven for an additional two minutes to equilibrate. Each sample was then removed from the oven. The dimensions of each sample were measured using hand-held dial calipers. The thermal expansion results are shown in Table 3. Thermal expansion of the material was found to be anisotropic at the elevated temperatures. This finding plays a major role in determining how the sheets deform when subjected to sagging under their own weight or subjected to a traction force such as air pressure. The effects of thermal expansion were incorporated in an approximate manner through the thermal expansion tensor (1).
METHODOLOGY OF SOLUTION AND FINITE ELEMENT FORMULATION
The constitutive equation (Eq 33, Ref. 1) is used along with the equation of virtual work (Eq 13, Ref. 1) and the Finite Element Method (FEM) to analyze the three cases presented in the Introduction. Nine-noded quadratic isoparametric elements were used with biquadratic shape functions employed to capture the curvilinear nature of the displacements. Tables 1-3 contain the material parameters with Poisson's ratio set to 0.35. The original coordinate system was Cartesian; however, the [x.sub.1] and [x.sub.2] displacements are both functions of [[zeta].sup.1] and [[zeta].sup.2]. The resulting set of equations were nonlinear. A globally convergent method was utilized to solve for the nodal displacements (7).
The element length was varied geometrically along the length of the sheet. Varying the element length allows one to "concentrate" more smaller elements where they are required. A great deal of refinement was required near the end points to best capture the behavior at each end. In this work, a total of 70 elements were utilized along the length with a symmetric distribution about the center-line. Meanwhile, two elements were employed through the thickness. Four-point Gaussian quadrature was used to integrate along the thickness and length of each element. Techniques for integration in curvilinear coordinates described by Narasimhan (8) were employed. The stresses used in the hereditary integral were calculated using the superconvergent patch recovery (SPR) method explained by Zienkiewicz and Zhu (9, 10). This method involved calculating the metric tensor components at the Gauss points for two consecutive elements along the length (a patch), and then calculating a "best curve fit through these points. The str esses and strains are then recalculated using this "best fit" curve. A cubic curve fit was used to calculate the shear, axial, and thickness stresses and strains. Calculations were made at each set of Gauss points along the length. Since two elements were used through the thickness, SPR calculations were done for eight sets of Gauss points, with each set representing the behavior at a thickness point.
The hereditary integrals were evaluated using Simpson's rule when an even number of time steps existed. In the case of an odd number of time steps, the trapezoidal rule was used to evaluate the integral for the first time step between t' = 0 and t' = [DELTA]t. The subsequent time intervals were evaluated using Simpson's rule. This approach tends to minimize errors of integration since the first time step (0 [less than or equal to] t' [less than or equal to] [DELTA]t) involves the most significant fading of the memory function.
SOLVING THE PROBLEM OF SIMPLE EXTENSION
Initially, in order to validate the finite element formulation and implementation, the problem of simple extension discussed in Part A (1) was solved using the finite element geometry defined in the previous section. The physical components of the stress and strain satisfied the analytical results (1). A graph of the variation of the physical components of strain [[epsilon].sub.11] and [[epsilon].sub.22] with the dimensionless physical component of elongation stress [[sigma].sup.11]/[lambda] exactly matched the analytical results. Figure 2 contains the variation of the physical components of strain with the applied load assuming a Poisson's ratio of 0.3 and a shear modulus equal to the sum of the shear moduli of each Maxwell element in Table 1 (460,076 Pa). These results serve to verify the validity of the finite element solution.
POLYMERIC BEAM EXAMPLE PROBLEMS
The two-dimensional curvilinear finite element formulation was solved at 150[degrees]C for the three different cases described in the Introduction. For all cases, the sheets were clamped at both ends and the temperature expansion was assumed to be the same through the thickness and length. For Case 1, the specific gravity of the styrenic material was assumed to be 0.96. The applied surface traction load was 5.171 Pa (7.5 x [10.sup.-4]psi) for both the Case 2 Maxwell fluid and Case 3 Maxwell solid analysis. This load was applied normal to the sheet on the upper surface. A time step of 0.0025 second was used to simulate the response at 150[degrees]C for all three cases. A few sample figures representing the displacements, moments, stress, and strain are presented here.
NODAL DISPLACEMENTS RESULTS
The center-line displacements for all three cases are shown in Figs. 3 and 4. The shape changes dramatically at the endpoints because of the clamped boundary condition. Application of surface traction gives a circular displacement pattern, whereas the displacement pattern of the sheet with only the body force acting is more U-shaped. Surface traction forces are always applied perpendicularly to the sheet, at all times forcing the material to displace outward away from the center-line at x = 0. An example of the deflection curve at various times for the case of surface traction is shown in Fig. 5.
A deflection curve was obtained at t = 0 because the Maxwell model gives an initial elastic response under the instantaneously applied surface traction. Use of the Jaumann rate derivative would necessitate integration through many pseudo time steps to compute the elastic response. The approach presented in this paper allows one to go directly to the solution.
STRESS AND DEFORMATION RESULTS
Elongation (axial) strain results are found in Figs. 6 and 7. The results shown were calculated at the Gauss point just above the center-line. Application of a body force load yields the result that the material near the end of the beam is stretched more than at the middle ([[zeta].sup.1] = 0) because the material at the clamped end is more in line with gravity (Fig. 6). The case of application of the surface traction load is shown in Fig. 7. The larger strains are produced at the middle of the sheet ([[zeta].sup.1] = 0), where the deflection is the greatest for the case of surface traction. In this case, the surface traction load is normal to the sheet (beam). This stretches the material away from the center-line ([[zeta].sup.1] = 0), resulting in a larger change in shape (Fig. 7).
Variation of the beam moment along the length is shown in Fig. 8. Moment is a measure of the resistance to bending. It is highest at time zero, when the initial response of a Maxwell material is elastic. As viscous effects dominate with time, the stresses in the material relax and the resistance to bending decreases with time.
The variation of the physical component of shear strain ([[epsilon].sub.12]) at the Gauss point Just above the centerline Is shown in Fig. 9. There is wide variability in the shear strain results near the clamped ends. The shear strain is largest near the clamped ends and close to zero at the center ([[zeta].sup.1] = 0). The abrupt change in shape due to the clamped boundary condition causes the largest amount of shear deformation.
END EFFECTS AND STRESS SINGULARITIES
In a two- or three-dimensional model, a stress singularity exists at the clamped end at the top and bottom corners. This can cause sharp waviness near the ends (Fig. 9). Fluctuations are also seen in the variation of the elongation stress ([[sigma].sup.11]) with material location for the instantaneous elastic response of the Maxwell Fluid model (Fig. 10).
As time elapses, the variability ("noise") became more significant. Use of the hereditary integral to compute the stress requires storage of the strain history; therefore, any initial variability would never disappear and would be added because of the hereditary integral. The increased variability at the end points with increasing time is also likely to be a function of the change in shape of the sheet at the ends.
The effect of stress singularities (discontinuities) at the boundaries has received some attention in the literature (11-19). Attention mainly focuses on elastic problems in a small deformation context. A case can be made that many ignore and do not report these problems. Here the unfiltered finite element solutions are reported.
To obtain a rigorous solution to the problem, one generally represents the displacements near the singular point in polar coordinates. The domain at a singular point is shown in Fig. 11 (16, 17). The region [[omega].sub.R] is closest to the singularity, the region represented by [[omega].sub.o] is removed from the singularity, and the region [[omega].sub.p] is a transition region between [[omega].sub.R] and [[omega].sub.o]. The region [[omega].sub.o] is represented by the usual linear or quadratic isoparametric elements. The region [[omega].sub.R] surrounds the singular point. The displacement near the singularity is usually represented by a known function in polar coordinates:
u = [[sigma].sup.[infinity].sub.i=0] [D.sub.i][f.sub.i](r, [theta], [[gamma].sub.i]) (5)
where [D.sub.i] is an expansion coefficient, [f.sub.i] is a known function, and [[gamma].sub.i] is a constant. The transition region [[omega].sub.p] provides a smooth match between the asymptotic expansion used over [[OMEGA].sub.R] and the usual isoparametric element functions used over [[OMEGA].sub.o]. Alternatively, an approach based upon the theory of boundary eigensolutions (18, 19) can be employed. This approach also requires some knowledge concerning the form of [f.sub.i].
The function [f.sub.i] in Eq 5 has been derived for some boundary value problems in the small deformation elastic case. Typically, [f.sub.i] has a dependence of [r.sup.[gamma]i]. In many cases, [[gamma].sub.i] is complex. An expression for [f.sub.i] needs to be derived for the case of finite deformations in a convected coordinate system. This is beyond the scope of this paper and will be the subject of a subsequent work.
CONCLUSIONS AND RECOMMENDATIONS
The convected rheological model developed in Part A was used with the curvilinear equation of virtual work to successfully solve a practical problem of a viscoelastic beam allowed to sag under its own weight. This has applications in the thermoforming process. The effect of thermal expansion was also incorporated into the generalized Hooke's law using a material (convected) coordinate system. The model predicted that the material would produce larger strains and have smaller stresses as time elapsed. This is physically intuitive. The stresses must decrease (relax) to allow the material to flow to the deeper depths. The computational finite element solution of the problem of simple extension confirmed the validity of the convected coordinate system approach to describe large deformation problems.
Use of adaptive refinement would reduce but not eliminate the "noise" in the stress and strain results close to the damped ends. Singularity functions in a finite deformation context using a convected coordinate system need derivation. This would complete the solution for two- and three-dimensional problems.
It is believed that the convected constitutive model would work for a variety of cases, including three-dimensional (3-D) formulations. A 3-D model would show deformation changes not only along the length and thickness, but also along the width.
In summary, the convected Maxwell rheological model formulated from the generalized form of Hooke's law was used in conjunction with the curvilinear form of the equation of virtual work to solve several very practical model problems. This formulation is fairly straightforward to implement and is applicable for the simulation of many practical processes that involve use of viscoelastic materials.
[Figure 1a omitted]
[Figure 1b omitted]
[Figure 2 omitted]
[Figure 3 omitted]
[Figure 4 omitted]
[Figure 5 omitted]
[Figure 6 omitted]
[Figure 7 omitted]
[Figure 8 omitted]
[Figure 9 omitted]
[Figure 10 omitted]
[Figure 11 omitted]
Table 1 Generalized Maxwell Fluid Model Parameters ([T.sub.o] = 150[degrees]C). Element Number [G.sub.i] (Pa) [[xi].sub.i] (s) 1 6410 19.823 2 21,281 2.486 3 61,046 0.319 4 85,434 0.041 5 291,905 0.003 Table 2 Generalized Maxwell Solid Model Parameters ([T.sub.o] = 150[degrees]C). Element Number [G.sub.i] (Pa) [[xi].sub.i] (s) 1 10,108 9.978 2 30,459 1.223 3 65,926 0.191 4 81,499 0.027 5 312,543 0.002 6 1172 N/A Table 3 Thermal Expansion Coefficient of the Styrenic Sheets. Machine Transverse Temperature Direction Direction 140[degrees]C -0.1419 0.00918 150[degrees]C -0.2180 0.04821
(1.) M. J. Stephenson and G. F. Dargush, Polym. Eng. Sci., this issue.
(2.) R. B. Bird, R. C. Armstrong, and O. Hassager, Dynamics of Polymeric Liquids. 2nd Ed., John Wiley & Sons, New York (1987).
(3.) S. L. Rosen, Fundamental Principles of Polymeric Materials, John Wiley & Sons, New York (1993).
(4.) A. C. Papanastasiou, L. E. Scriven, and C. W. Macosko, J. Rheol., 27, 387 (1983).
(5.) D. W. Marquardt, J. Soc. Indust. Appl. Math., 11, 431 (1963).
(6.) M. J. Stephenson, M. E. Ryan, and G. F. Dargush, Polym. Eng. Sci., 39, 2199 (1999).
(7.) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes In Fortran, Cambridge University Press, London (1992).
(8.) M. N. L. Narasimhan, Principles of Continuum Mechanics, John Wiley and Sons, New York (1993).
(9.) O. C. Zienkiewicz and J. Z. Zhu. Int. J. Num. Meth. Eng., 33, 1331 (1992).
(10.) O. C. Zienkiewicz and J. Z. Zhu. Int. J. Num. Meth. Eng., 33. 1365 (1992).
(11.) M. L. Williams, J. Appl. Mech., 12, 526 (1952).
(12.) R. Wait and A. R. Mitchell, J. Comp. Phys., 8, 45 (1971).
(13.) B. G. Prakash and R. Ramchandra, J. Math. Phy. Sci., 10, 333 (1976).
(14.) R. Ramchandra and B. G. Prakash, J. Struct. Mech., 6, 133 (1978).
(15.) T. M. Hrudey and M. M. Hrahock, J. Eng. Mech, 112, 666 (1986).
(16.) Z. Yosibash and B. Schiff. Int. J. Fract, 62, 325 (1993).
(17.) Z. Yosibash and B. Schiff, Finite Elem. Anal. Des., 26, 315 (1997).
(18.) A. R. Hadjesfandiari and G. F. Dargush, Int. J. Num. Meth. Eng., 50, 325 (2001).
(19.) A. R. Hadjesfandiari and G. F. Dargush, J. Appl. Mech., 68, 101 (2001).
M.J. STEPHENSON (1)(*) and G.F. DARGUSH (2)
(*.) Corresponding author.
(1.) Research and Development Group Aristech Acrylics LLC Florence, KY 41042
(2.) Department of Civil Engineering State University of New York at Buffalo Buffalo, NY 14260
|Printer friendly Cite/link Email Feedback|
|Title Annotation:||Development of a curvilinear viscoelastic constitutive relationship for time-dependent materials, part B|
|Author:||Stephenson, M.J.; Dargush, G.F.|
|Publication:||Polymer Engineering and Science|
|Article Type:||Statistical Data Included|
|Date:||Mar 1, 2002|
|Previous Article:||Development of a curvilinear viscoelastic constitutive relationship for time dependent materials. Part A: Theoretical discussion.|