FEA of diffusion-reaction in tires.

There is considerable interest on the part of tire manufacturers to improve the long term performance of tires. One consequence of this desire is the recent development of the 80,000 mile passenger tire. Accelerated wheel tests such as stepped-up speed and stepped-up load are common approaches to determine tire life. However, by virtue of their accelerated nature, these tests cannot totally reflect the long-term serviceability of tires. Diffusion of oxygen through the innerliner and other plies supports oxidative reactions in the body of the tire (ref 1) gradually degrading the material properties of the internal components. A critical area in this respect is die belt edge due to its flex cycling and heat history.

Researchers refs. (2 and 3) have studied the aging characteristics of tire compounds, particularly the belt skim, and related it to the durability of the tire. Tokita, (ref 4) et al, developed a finite difference model to predict service life based on the oxidative characteristics of the tire components. He used a finite difference approach to solve the plane problem of steady state diffusion oxidation reaction. His work is based on directly solving the steady state diffusion-oxidation equation using concentration as the primary unknown and assuming a first-order oxidation reaction. In his work, the problem was discretized using finite difference approximations in plane, polar coordinates to characterize the curvature of the sidewall of the tire. The concentrations are, in general, discontinuous across ply layer boundaries due to differences in solubilities. However, for the steady state problem any difference in the fluxes across the ply boundaries can be attributed to oxidation. Tokita used this condition to write flux balance equations at the nodes along ply boundaries.

In the present work, to avoid the complications of a discontinuous unknown, the governing equations are reformulated to yield partial pressure as the primary unknown. This article will discuss the method used to reformulate the governing equations and their discretization using the finite element method as well as some preliminary results. A finite element program called Femoxi was developed. It is capable of solving transient, plane, two-dimensional, and two-dimensional axisymmetric problems of diffusion with reaction. Additional software was developed to permit the use of a digitizer and a commercial finite element model builder and post processor program. The entire system operates on a personal computer and is currently being used to study oxygen diffusion and reaction in tires. A model to predict tire life is currently in the early part of validation and will be discussed briefly. This software system can also be used to study the diffusion of air to predict intra-tire pressures and to study the diffusion of moisture and subsequent reaction with the steel cords.

Governing equations

For generality in application, the governing equations have been developed for both the plane and two dimensional axisymmetric coordinate systems as shown in figure 1.

In the following, we develop the governing differential equations for the plane problem in X-Y coordinates shown in figure 1. The development of the axisymmetric capability shown in figure I is discussed during the presentation of the finite element formulation. Looking forward to a finite element formulation, we assume that diffusivities are constant over any given finite element (a finite element may be thought of as a small piece of the structure of unit thickness in the case of the plane problems and a small ring-shaped piece of the tire cross-section in the axisymmetric problem.

Because certain tire plies have strong directional properties, we treat them as orthotropic in the ply direction and the direction perpendicular to it. Since the finite element solution is referenced to the X-Y or R-Z planes, we must provide for the off diagonal terms of the diffusivity tensor which occurs in other than principal coordinates. We further assume that the oxidation reaction is governed by a first order rate equation. With these assumptions the governing equation of diffusion-oxidation in the tire is given by

[Mathematical Expression Omitted]

In this equation C is the concentration, [D.sub.x], [D.sub.xy], [D.sub.yx] and [D.sub.y] are the components of the diffusivity tensor, k is the first order reaction rate constant and t is time. We assume both diffusivity tensor and reaction rate constant to be independent of concentration. Van Amerongen (ref. 5) justifies this assumption for the diffusion constants on the basis that at normal pressures only small amounts of the gas are taken up by the rubber. It is assumed that the same argument can be extended to the reaction rate constant.

Direct implementation of equation 1 using the finite element method poses a problem. The concentrations have discontinuities across the boundaries of the plies of the tire due to the different solubilities of the plies. Therefore, we reformulate equation 1 looking forward to finite element discretization. The finite element model will be constructed so that each finite element will contain only one rubber component. On this basis assume that the solubility is constant within any individual element. If we take account of Henry's Law, P = c/s, where P is partial oxygen pressure and s is solubility, we find after dividing each term of equation 1 by element solubility

[Mathematical Expression Omitted]

For most problems in tires, the time to equilibrate via diffusion after initial inflation is very small relative to long term durability considerations. In such cases we may take (Mathematical Expression Omitted] = 0. Thus, the steady form of equation 2 is

[Mathematical Expression Omitted]

With the assumptions discussed above, equations 2 and 3 are the reformulation of equation 1 which is valid over any finite element volume and suitable for implementation using the finite element method. Since the pressure variable is continuous across plies of the tire, these equations apply over the entire problem domain. That is, to the assemblage of all the individual finite elements which represent the tire.

In order to calculate the degree of reaction at any point in the tire cross section we use the pressure output of the finite element program. For the transient equation 2, the pressure is output by Femoxi for a series of time steps at which the problem is solved. The extent of reaction is given by the following equation

ks [integral of][.sub.o].sup.t] P(t) dt = A which is integrated numerically.

In this equation, A represents the degree of reaction which has taken place. For example, depending on units used, the value of A at time t might be expressed in terms of moles of oxygen reacted per cubic millimeters of rubber. If a critical value of the amount of oxygen reacted can be determined from tests, equation 4 can be used to predict tire failure due to the oxidation reaction.

As mentioned earlier, for most tire problems it appears that the problem may be considered as steady state. For the steady state problem, the degree of reaction, A, is given by

A = ksPt (5)

Using a critical value of say oxidation, we can rewrite equation 5 as [t.sub.c] = [A.sub.c]/ksP (6) where [A.sub.c] is the critical value of A and [t.sub.c] is the time to reach this critical value.

Determination of diffusivity constants

Diffusivity constants are sometimes determined using permeability tests. In the case of an orthotropic material, the principal diffusivities [D.sub.i] may be found from [D.sub.i] = [P.sub.i]/S i = 1,2 (7) where [P.sub.i] are measured permeabilities and subscripts 1 and 2 represent principal directions. It is assumed that equation 7 is valid when reaction rates are low relative to diffusion rates. This has been found to be die case in our studies. Figure 2 shows the finite elements implemented in the finite element program Femoxi. In these elements the small numbered dots are called nodes. The X-Y and R-Z coordinate frames are called the global coordinates. The x-y coordinate frames are called local or element coordinates. By aligning the elements along ply layers in the tire, the local x-y coordinate system shown in figure 2 will coincide with the principal directions of the diffusivity tensor. Since material properties are entered separately for each element and since diffusivities are entered in terms of the local element coordinate system, only the principal diffusivities, which can be determined from equation 7, are required as input to the program. The direction of the local x coordinate is determined by the nodal ordering for the element. This is called the Boolean connectivity. For example, in figure 2, the connectivity sequence 1-2-3-4 determines the local x direction (1-2 direction). The local y direction is 90 degrees counter-clockwise to the local x direction.

Variational formulation of governing equation

We illustrate the procedure of finite element formulation by first developing the variational form of equation 3. The temporal problem of equation 2 was solved using standard finite difference techniques for time dependent equations (ref. 6). It is not discussed further in this article since such techniques are well known and so far only steady diffusion-oxidation problems have been studied.

Using standard finite element procedures (ref. 6) we multiply equation 3 by a test function v and integrate by parts over the element area Ar (we have assumed the element thickness to be unity). This serves to transfer one order of differentiation from the variable P to the variable v, yielding the variational form of equation 3 as

[Mathematical Expression Omitted]

Here, [q.sub.n] is the flux around the boundary of the element

[Mathematical Expression Omitted]

where [n.sub.x] and [n.sub.y] are the components of the outward normal vector n shown in figwe 2. Here the symbol ds is the increment of boundary along the edges of the element. The variable [q.sub.n] is the outward flux of the diffusing fluid of the problem being solved, e.g., oxygen flux. Where elements are joined during assembly all such fluxes except those on the outer boundary will cancel.

Finite element formulation in x-y

Next, we expand the unknown pressure, P, in terms of shape functions (ref 6) which span the finite elements. These shape functions [N.sub.i] have the property that they take on the value 1 when i = n, where n is the node number, and [N.sub.i] = 0 for the other nodes. By way of illustration, figure 3 shows the four shape functions for the four node isoparametric (ref. 7) quadrilateral which is one of the elements implemented in Femoxi.

Expressing P in this manner, we have

[Mathematical Expression Omitted]

where P(x,y) is the interpolated value of P within the element, [P.sub.i] is the value of pressure at node j, [N.sub.j] is the value of the shape function [N.sub.j] at node j, and n is the number of nodes in the element.

As mentioned earlier, we assume that the diffusivity tensor is symmetric, i.e., [D.sub.xy] = [D.sub.yx]. With this assumption, and, substituting equation 10 into the variational form of equation 8 we have equation 11 shown in table 1.

We may also write equation (11) as

[Mathematical Expression Omitted]

where

[Mathematical Expression Omitted]

Isoparametric finite element formulation

We next discuss the transformation from the global X-Y coordinate [integral of] system to a natural coordinate system [xi]-[eta] which has its origin at the center of the element and is scaled so that [xi] and [eta] range between + 1 and - 1.

We now convert equation 11 from an integral over x-y to an integral in terms of the natural coordinates [xi]-[eta]. This is necessary since it would be difficult to integrate equation 11 if the element had irregular boundaries as in figure 2. This is accomplished by a mapping of the element from the global X-Y frame to die master element in [xi]-[eta] space as shown in figure 4. The master element is square and has comer coordinates between +1 and -1. This facilitates the subsequent numerical integration.

To accomplish this transformation we expand the x and y coordinates in terms of the same shape function we used for the pressure p in equation 10. Thus we have

[Mathematical Expression Omitted]

The term isoparametric derives from the fact that the same shape functions are used for interpolation of both pressure and geometry. Equation 16 permits expression of the partial derivatives of the shape functions Ni with respect to the [xi]-[eta] coordinates instead in equation 11. The elemental area term in that equation is now given by dXdY = Jd[xi] d[eta] where J is called the Jacobian of the transformation between X-Y and [xi]-[eta]. It is the determinant of the Jacobian matrix [J],where

[Mathematical Expression Omitted]

([J.sup.*].sub.11/[J.sup.*].sub.21]) ([J.sup.*].sub.22]) = [J].sub.-1] (16) Next, we define the elements of the inverse of [J] using the asterisk. ([J.sup.*].sub.11][J.sup.*].sub.21]([J.sup.*].sub.22]) = [J].sup.-1] (18) With this, the derivatives of the shape functions may be expressed in terms of [xi]-[eta] as in equation 19 (top of next page. Using equations 16 and 19, equation 13 can now be integrated in terms of the local [xi]-[eta] coordinates as equation 20 in table 1, where Ar is now the area in [xi]-[eta] space. At this point the integration is carried out mumerically over the unit rectangle in the [xi]-[eta] space.

Boundary conditions 

The primary variable of the problem is the pressure. This is the essential (ref 7) boundary condition and must be specified at least at one point (node) in the model. The secondary variable is the flux [q.sub.n] which need not be specified. The flux is the natural (ref. 7) boundary condition of this problem. In applying Femoxi to the oxidation problem, internal tire partial pressure of oxygen is specified along the inside of the tire and atmospheric partial pressure of oxygen is specified along the outside boundary. A tapered pressure transition is used at the rim. A zero flux boundary condition is used along the symmetry line.

Axisymmetric formulation

The axisymmetric version of equation 8 is developed by considering a one radian segment of the three elements discussed previously. By way of illustration, such axisymmetric elements are shown in figure 5. Also shown in figure 5 are the integration points at which the area integral of equation 11 is evaluated.

The integration is usually carried out using a numerical integration technique due to Gauss. Hence, these integration points are often referred to as Gauss points. For these elements we must now replace the unit thickness of in the plane problem characterized by equation 8 by the thickness at the Gauss integration points. For a one radian segment, this is simply the respective radii at the various integration points. The products of the derivatives of the shape functions [N.sub.i], which are evaluated at the Gauss point, are multiplied by the radii at their respective points.

For the three node triangular element, the integration is carried out in closed form rather than by numerical integration. In this case, the appropriate terms of equation I I are multiplied by the average radius of the element. This is found by averaging the radii of the three nodes defining the element.

Fnite element program Femoxi

Using the procedure discussed, a special purpose finite element program called Femoxi was coded in Fortran. Other than the element formulation, assembly, application of constraints, solution and recovery of secondary variables, all involve standard finite element technology. The program solves both the transient and the steady problem of diffusion reaction. The user may employ eight node quadrilaterals, four node quadrilaterals and three node triangular elements in either a plane Cartesian or axisymmettic coordinate frame.

Software was also developed linking Femoxi to a digitizer and a commercial PC finite element preprocessor and postprocessor (ref 8). The digitizer is used to locate the nodes of the finite element model using a large drawing or photograph of the tire cross section. The preprocessor is used to develop the finite element model of the tire. It aids inputting both geometrical and material data. For example, each element may have different properties such as

principal diffusivities and reaction rate constants. The postprocessor is used to display results. Typical display results are contours of pressure throughout the tire cross-section.

Although the software provides for solution of transient problems, this feature has not yet been fully explored.

Development of experimental data

To use the model, diffusion coefficients and oxidation rate constants are needed for the tire components. Initial work has focused on calculating the diffusion coefficients for air as a function of temperature for several model tire compounds. For validation purposes, the model is initially being used to predict the pressure profile through the cross section of the tire under static, non-oxidative conditions to detemiine equilibrium breaker and carcass pressure.

The diffusion coefficients are being calculated using the time lag approach (ref 9). Rubber membranes are subjected to an upstream pressure of 330 KPa, and the downstream pressure is monitored as a function of time until a steady state pressure slope is reached. The steady state slope is extrapolated to the time axis of a pressure versus time plot. The intersection is called the time lag, [theta], and the corresponding diffusion coefficient is given by D = [L.sup.2]/6[theta] (21) where L is the membrane thickness. Using this approach, typical air diffusion coefficients have been found to be 0.2 x 10-6 [cm.sup.2]/sec for a 100-phr butyl tube and a halobutyl innerliner and 2.5 x [10.sup.-6] [cm.sup.2]/sec for a black sidewall. For the model validation currently in progress, air diffusion coefficients for several model tire compounds have been calculated for various temperatures.

The calculation of oxygen diffusion coefficients is currently in progress using an alternative approach. Additionally, since air is approximately one-fifth oxygen, the oxygen and air diffusion coefficients should be related by a common factor.

An appropriate technique to calculate both the oxidation constants and critical oxidation values is still being developed. Initial success has been found with a microtome, infrared technique. Tire tests have shown that failures result when an extent of oxidation equal to about twenty-four percent change in a certain wavelength occurs in the belt-edge region. Studies are currently in progress relating loss of tensile strength properties to the change in this particular wavelength.

Analysis of a typical automobile tire

The development of tire models for analysis followed the procedure outlined by R. Probhakaran (ref. 10). Two-dimensional tire layout drawings either carefully hand sketched or taken from prints were used for entry of nodes and coordinates. Depending upon the size and the complexity of the cross section, models have varied in size between one hundred and two hundred, four-node quadrilateral elements. A typical problem executes in a few minutes.

To date, work has focused on using the model to predict the pressure profile in tires under static, room temperature conditions. Figure 6 shows two such profiles. Figure 6 is that of a half tire cross section with a 100 - phr halobutyl innerliner and a general-purpose rubber innerliner. In both cases, the inflation pressure was 210 KPa. The halobutyl innerliner results in lower intra-tire pressures because of its ability to better resist the diffusion of air. Such models predict important tire properties such as carcass and breaker pressures, both of which have been related to tire durability. Work still in progress is directed toward modelling of the full diffusion-oxidation reaction for various tire designs.

Conclusions

Problems of diffusion-reaction in tires are amenable to solution using the finite element method. Such problems typically arise in studies of oxidation of the belt edge and corrosion of steel reinforcement. Adopting a finite element approach, a new computational tool called Femoxi has been developed. Femoxi runs on a PC and solves such problems in several minutes. Using this tool, design alternatives aimed at alleviating such problems can be efficiently and quickly evaluated. Current effort is directed toward development of reaction rate constants and further validation of the model.

References

[1.] D.M. Coddinton, Rubber Chem. Technol 52, 905 (1979). [2.] A. Ahagon, M. Kida and H. Kaidore, Rubber Chem. Technol 63, 683 (1990). [3.] H. Kaidore and A. Ahagon, Rubber Chem. Technol. 63, 698 (1990). [4.] N. Tokita, W.D. Sigworth, G.H. Nybakken and G.B. Ovyang, "Long term durability of tires," Int. Rubber Conf., Kyoto, Oct. 15-18, 1985, pp. 672-679. [5.] G.T Van Amerongen, Rubber Chem. Technol. 37, (5), (Dec. 1964). [6.] J.N. Reddy, "An introduction to the finite element method, McGraw-Hill Book Company, New York (1984). [7.] O. C. Zienkiewicz and R.L. Taylor, "The finite element method, 4th Ed., Vol. 1, McGraw-Hill Book Company, New York (1988). [8.] H. Kardestuncer, Ed., "Finite element handbook," Chaps. 3 and 4, McGraw-Hill Book Company, New York (1987). [9.] G. T Van Amerongen, "The permeability of different rubbers to gases and its relation to diffusivity and solubility," J. of App. Physics, p. 972 (1946). [10.] R. Probhakaran, "Interactive graphics for the analysis of tires," Tire Science and Technology, TSTCA Vol. 13, No. 13, pp. 127-146 (July-Sept. 1985).