CONTINUUM THEORY FOR GAS - SOLID- LIQUID MEDIA A Duserfahcn Ear the Degree of DH. D. MECHIGAN STAR UNWERSIT‘! Robert John Gustafson €974 "hi: <3 vs University This is to certify that the thesis entitled CONTINUUM THEORY FOR GAS-SOLID-LIQUID MEDIA presented by Robert John Gus tafson has been accepted towards fulfillment of the requirements for PhD degree in M Michigan Start: 0": r. - ..__’ ‘u. : ABSTRACT CONTINUUM THEORY FOR GAS-SOLID-LIQUID MEDIA By Robert John Gustafson The objective of this work was to develop a continuum model which adequately describes the mechanical behavior of a biological product. The medium was assumed con- structed of two sets of interconnected pores separated by a solid. One set of pores contained a gas; the other a liquid. Constitutive and stress-strain relations were derived through energy considerations for the gas-solid- liquid medium. A series of loading conditions was outlined for determination of the material constants. Experimental attempts were made to determine the compressibility of apple parenchyma using three loading conditions which involved the internal gas pressure and a hydrostatic pressure. Apparatus limitations prevented successful determination of the compressibilities. Suggestions for improvement of the apparatus were made. Robert John Gustafson The finite element method was used to obtain numerical solutions to the axisymmetric boundary value problem for the gas—solid-liquid medium. The stress distribution in a Spherical body without a skin (un- restrained), with a skin (restrained), and a restrained body subjected to flat plate compression were studied. Internal liquid pressure was varied between 689.5 kPa (100 psi) and 2758 kPa (400 psi) using 689.5 kPa (100 psi) increments while a zero gas pressure was used. Material properties of an apple were used in the contitutive equations. An unrestrained homogeneous body with the liquid under pressure expanded without developing stress. The restrained body was found to have a hydrostatic (com- pressive) stress in the parenchyma and tension stress in the skin. Flat plate compression combined with the liquid pressure produced shear stresses and a hydrostatic stress in the parenchyma. I M. Approveo Department Chairman CONTINUUM THEORY FOR GAS-SOLID-LIQUID MEDIA BY Robert John Gustafson A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Agricultural Engineering 1974 ACKNOWLEDGMENTS The author sincerely appreciates the assistance of all who have aided in this study. He is especially appreciative, however, of the counsel and guidance pro- vided by his major professor, Dr. Larry J. Segerlind (Agricultural Engineering), during this graduate program. To the other members of the guidance committee, Dr. G. E. Mase (Metallurgy, Mechanics, and Material Science), Dr. G. E. Merva, and Dr. J. B. Holtman (Agricultural Engineering), the author expresses his deepest gratitude for their time, professional interest, and constructive suggestions. Appreciation is also expressed to the Agricultural Engineering Department for the assistantship and other financial support given. 11 TABLE OF CONTENTS LIST OF FIGURES ABBREVIATIONS AND SYMBOLS I. INTRODUCTION II. REVIEW OF LITERATURE 2.1 Mechanical PrOperties of Plant Materials . . . . 2.2 Biot Theory of Elastic Porous Media III. DEVELOPMENT OF QUASI-STATIC THEORY OF A GAS-SOLID-LIQUID SYSTEM . . . 3.1 Model Description 3.2 Theoretical Development 3.3 Linear Stress-Strain Relations 3.4 Governing Equations for Transient Phenomena 3.5 Determination of Elastic Coefficients 3.5.1 Jacketed Compressibility Test - No Liquid Flow . 3.5.2 Jacketed Compressibility Test - No Gas Flow . . 3.5.3 Jacketed Compressibility Test - No Liquid Flow, No Gas Pressure . iii 10 13 13 15 21 29 33 34 36 37 Page 3.5.4 Jacketed Compressibility Test - No Gas Flow, No Liquid Pressure . 38 3.5.5 Jacketed Compressibility Test - No Liquid Flow, No Gas Flow . . 40 3.5.6 Alternate Procedure for Determination of M . . . . . . . 41 3.5.7 Alternate Procedure for Determination of N . . . . . . . 43 3.6 Closure . . . . . . . 45 IV. EXPERIMENTAL DETERMINATION OF COEFFICIENTS a AND M-FOR APPLE FLESH . . . . . . 48 4.1 Equipment and Procedure . . . . 48 4.2 Results of Experimentation . . . . 53 4.3 Discussion of Apparatus Limitations . . 54 V. FINITE ELEMENT FORMULATION . . . . . 56 5.1 Development of the Variational Equations . 59 5.2 Implementation of the Method . . . 65 5.3 Results of Finite Element Analysis . . 69 5.3.1 Unrestrained Sphere . . . . 69 5.3.2 Restrained Sphere . . . . 70 5.3.3 Flat Plate Contact . . . . 70 5.4 Closure . . . . . . . 74 VI. SUGGESTIONS FOR FURTHER STUDY . . . . 78 BIBLIOGRAPHY . . . . . . . . 81 iv Figure 301 01010101009- (J'lIvD-OJNHN LIST OF FIGURES Differential Element of Medium Equipment for Measurement of Bulk Compressibility Bulk Compression Apparatus Simplex Triangular AxiSymmetric Element Outline of Finite Element Computer Program Grid for Finite Element Application Hydrostatic Stress in Parenchyma of Body Maximum Principal Stress and Shear Stress in the Skin Flat Plate Contact - Distribution Curves for Tmax (kPa) Flat Plate Contact - Distribution Curves for Principal Stress (kPa) Suggested Bulk Compressibility Apparatus Page 14 50 51 57 66 67 71 72 75 76 80 8.1, 8.2, 33 b b2, b3 1’ Cl, 02, C3 ABBREVIATIONS AND SYMBOLS area ratio of gas pressure to applied hydro— static pressure coefficients for shape functions matrix defined by {e} = [Q [U] material prOperty compressibility of the liquid compressibility of the gas degree Centigrade centimeter material properties matrix material property material property - modulus of elasticity dilatation of bulk material dilatation of the liquid dilatation of the gas dilatation of the solid strain tensor strains in the coordinate directions vi traction on the liquid liquid porosity tractions on the solid element force matrix traction on the gas gas porosity material property strain invariants inch bulk modulus of elasticity - ”open system” bulk modulus of elasticity — ”closed system” liquid permeability gas permeability element stiffness matrix kilopascal = (N/m2) x 103 matrix of shape functions material property milimeter material property shape function unit normal vector force matrix material property hydrostatic pressure pressure of liquid pressure of gas vii W1 or {W} a 0.1, 0.2, (13 81, 82, 83 51, 52, 53 $1! ¢21 ¢3 pounds per square inch displacement of liquid stress component in r direction stress component in z direction surface strain energy displacement matrix displacement of gas displacement in the radial direction displacements displacements of solid matrix volume volume of bulk material volume of liquid volume of gas displacement in z direction gas flow relative to solid strain energy per unit volume work by concentrated forces work by applied stresses work by internal and applied loads liquid flow relative to solid material prOperty coefficients for approximating polynomials of the displacement components viii Yx: Y Y y, Yi Of 01.1 {r} TIJ material property shear strain components gas flow parameter liquid flow parameter jacketed compressibility - no gas flow jacketed compressibility - no liquid flow Kronecker delta, 0 for i f j, and 1 for i = J strain matrix relative gas flow unjacketed compressibility jacketed compressibility no liquid pressure no gas flow, jacketed compressibility no gas pressure no liquid flow, material property - Lamé elastic constant for "open system" material property - Lamé elastic constant for "closed system" material property - Poisson's ratio viscosity of liquid viscosity of gas total potential energy density stress component on liquid stress component on gas stress components on the solid stress matrix total stress components - bulk material ix normal stress components - bulk material shear stress components - bulk material effective stress volume of gas relative liquid flow jacketed compressibility - no liquid flow, no gas flow I. INTRODUCTION Extensive loss is often incurred in the production of various fruits; notably sweet cherries, tomatoes, and certain varieties of apples; as a result of cracking of the skin and fleshy tissues some time prior to harvest. Cracks affect the appearance of the fruit, encourage mold and insect contamination, cause trimming losses for the canner, and result in consumer dissatisfaction. Rapid deterioration usually follows the exposure of the ruptured tissue to air, and the injured fruit becomes worthless or of inferior grade. Cracked fruit, therefore, represent a considerable loss in income to both processing and fresh market industries. Most of the previous research on fruit cracking has been concerned with the environmental conditions conducive to fruit splitting and breeding crack-resistant varieties. It was the author's contention that present mechanics models of materials were inadequate for representation of the behavior of many biological products. A model which incorporates the effects of the liquid and gas elements was needed. The objective of this study was the develop- ment of such a model. 2 The work reported in this thesis may be divided into three parts: 1. The development of a model containing parameters and variables necessary to more adequately represent the mechanical behavior of fruit tissue. The experimental determination of the parameters necessary for use of the model developed. Use of the finite element method for solution of the continuum equations yielding stress distribu- tions for a body approximating the shape of a fruit. II. REVIEW OF LITERATURE Two areas must be considered when developing a continuum mechanics model for a biological product such as a fruit: the plant physiology and horticultural area and the material science and mechanics area. This review of literature is divided into two sections: 1) that work related to plant materials and ii) the mechanics theory applicable to the present approach. 2.1 Mechanical Properties of Plant Materials Several authors have recently reviewed the available literature related to fruit anatomy pertinent to mechan- ical properties and splitting of fruits. Tennes (1973) reviewed literature available on fruit structure of tomatoes and cherries as related to fruit splitting. He also reviewed the cultural practices which have been used in an attempt to alleviate the problem. Brusewitz (1969) reviewed in detail the anatomy of plant material which might be pertinent to the mechanical properties of the material. Akyurt (1969) also discussed available literature on the cellular structure of plant material and modeling attempts which have been recorded. Because of these 3 4 studies, the present review of literature was restricted to works which have direct bearing on the proposed method of modeling the elastic action of a system of cells. Lit- erature on the growth of plants was not reviewed unless it is pertinent to the elastic action of the material. Many researchers pursuing the phenomenological approach to the mechanics of biological materials have assumed the tissues are continuous, homogeneous, and iso- tropic. Only a few researchers have attempted to describe the effect of such variables as cell dimensions and turgor pressure. Of those attempting to include such variables, the mechanical properties of plant materials have been studied on three levels: the cell wall, the cell, and tissue. The structure of the cell wall has been studied in great detail. A number of references, such as Frey- Wyssling (1952), describe the make-up of the cell wall. Probine and Prestone (1962) concluded that the anisotropic nature of the cell wall affects its mechanical properties and cell growth. Cleland (1971) concentrated on certain aspects of cell wall extension including the mechanical properties of primary walls and their relation to cell enlargement. He presented two conclusions from his rheological studies. First, the mechanical properties of all primary cell walls are probably qualitatively similar. Second, the difference 5 that exists between cell elongation and mechanical proper- ties of the walls are sufficient to indicate that cell wall extension does not involve simple physical stretching of the wall. Frey-Wyssling (1952) pointed out that mechanical stresses on tissues not only involve elastic or plastic alteration of the cell wall, but they also bring about morphological deformations of cell shape, which may be even more important than the rheological behavior of the cell wall. Most attempts at describing cell action have been based on the assumption of some regular geometric pattern and shape for the cells. Matzke and Duffy (1955), how- ever, stated that rigid conformity to pattern does not exist. They found the average number of faces per cell close to fourteen and that the faces vary in shape from triangular to nonagonal. Haines (1950) found that the relation between cell extension and turgor pressure for spherical isotropic cells obeying Hooke's Law is not linear but hyperbolic. He further states, "There can be no approximate linear relationship to any cell dimension satisfactory for purposes of calculating turgor or osmotic pressure." Broyer (1952) defined a coefficient of distention of a boundary as a = Z%;A and a coefficient of enlargement A P He calculated relations between and for 6 several simple geometric shapes and found that relative volume and changes in volume can be used to determine e as defined for actual pressure changes. Work required for change in relative volume was calculated assuming constant pressure. Philip (1958b) developed a dynamic theory of the classical osmotic plant cell in quantitative form and extended it to the case where diffusible solute is present. He assumed change in cell volume to be. linearly proportional to change in turgor pressure, and developed a first-order expression for the change in relative volume with respect to time for several initial and boundary conditions. In a second paper, Philip (1958a) developed relations for propagation of turgor pressure and other pr0perties through cell aggregates describing them mathematically as diffusion phenomena and assuming the elastic modulus of the material defined by l/E = l/T(V/VO -1) where T is the turgor pressure. He stated further study of the cell-wall stress distribution was needed to refine the propagation equations. Building on the work of Broyer and Philip, Slayter (1967) concluded that for given cell dimensions and permeabilities, both the elastic properties of the cell and the internal osmotic pressure influence the rate of swelling and shrinking. 7 Mela (1967a), using microscopic photography, developed a method for studying Young's modulus by extension of mitochondria cells measured for different concentrations of salt water. He assumed membranes of uniform thickness which do not have pores large enough to affect stress- strain calculations. In a second article, Mela (1967b) reported Young's modulus to be a non-linear function of temperature with a minimum at 12—1300. Studies must be done on the level of a multi-cell structure, or tissue, to include the effect of interaction between cells. These studies have been performed by both removing a segment of tissue for analysis and by studying the action of the whole body such as a fruit. Falk §t_al. (1958) studied the relations between turgor pressure and Young's modulus using the resonant frequency of potato tuber parenchyma. They concluded that the cell wall material follows Hooke's law and that there are changes in elasticity of the whole parenchyma due to turgor pressure which are in turn reversible thanks to the ideal cell wall material. Nilsson et_al. (1958) studied the dependence of Young's modulus of potato tuber parenchyma on turgor pressure using a simple theoretical model. The cells of the parenchyma were approximated by regular geometric cell-forms (Spheres or polyhedra), each cell being bounded by an elastic membrane and filled with an 8 incompressible fluid. It was shown that this model yields the correct dependence of cell diameter on turgor pressure and that certain cell-wall constants can be determined using the relation. Burstrom gt_gl. (1967) used the resonance frequency method for determining Young's modulus of internodes of etiolated pea seedlings. They found the modulus to increase nearly proportional to turgor pressure and that at water saturation the modulus is more than fifty times higher than at plasmolysis. Studies by Meynhardt (1964) indicated it may be possible to predict the susceptibility to splitting of different grape cultivars by an anatomical investigation of the berry tissue. He concluded that it seems possible that the subepidermal cell dimension ratio (longitudinal to radial) and the number of subepidermal cell layers may contribute to the resistance or susceptibility of grape berry tissue to Splitting. Working with tomatoes, Cotner g£_al. (1969) found that fruit with flattened epidermal cells were less susceptible to concentric cracks than those with rounded cells, but no such correlation existed for radial cracks. Clevenger and Hamann (1968) studied the mechanical properties of apple skin. They determined material proper- ties, including elastic modulus and Poisson's ratio, for three varieties of apples. All skins were found to be 9 anisotropic with the greatest strength in the longitudinal direction. Relaxation and creep experiments showed that apple skin tends to be viscoelastic in behavior. Four element (Burger) models were found to describe the action of the material very well. Akyurt (1969) and Akyurt gt_al. (1972) attempted to develop methods for studying the stress—strain relations in plant materials. With the cell wall idealized as a shell, the finite element method was proposed for the solution of the corresponding linear equilibrium problem. Akyurt showed that macrodisplacements as well as stresses and couple stresses acting on cellular bodies emerge as solutions of the field equations of the micropolar theory of Eringen (1962). The linear theory of viscoelasticity was also employed. Considine and Kriedeman (1972) devised a laboratory technique to measure the internal turgor pressure required for fruit rupture in order to assess resistance to splitting of grapes. Fruit of uniform maturity and known osmotic potential were immersed in a range of osmotics to create a known turgor pressure at equilibrium. "Critical turgor," the pressure which resulted in 50 percent of the grapes splitting, was approximately 15 atm in grape cultivars prone to splitting and 40 atm in resistant cultivars. They found splitting was not necessarily related to berry size or to the presence of seeds. No dominant relationship 10 was found to exist between the berry shape and the susceptibility to Splitting. The authors concluded that it is the epidermal and subepidermal layers which limit berry enlargement. Two uniform groups of berries were immersed in distilled water to emphasize this point. One group was intact, the other peeled. Within 30 minutes, the intact fruit were ruptured. Peeled fruit, on the other hand, absorbed twice as much water without suggestion of splitting. Raschke (1970) studied the transmission of changes in water potential in leaves. He found the epidermis of the leaf Zea mays transmits changes in water potential in the water supply of the leaf to the stomata within 0.1 second. Also, reduction in water supply can cause the subsidiary cells surrounding the stomata to collapse with- in 1.5 minutes, and the epidermis to shrink to one-third of its original thickness within 20 minutes. 2.2 Biot Theory of Elastic Porous Media Theories of deformation of a porous material contain- ing a viscous compressible fluid and the theory of flow of the fluid through the material have been developed and discussed in a series of papers by Biot and his co—workers (Biot, 1941, 1955, 1956b, 1962, 1963; and Biot and Clinger, 1941, 1942). The theories were first applied to consolid- ation and settlement of foundations for both isotropic and 11 anisotrOpic media. Later developments have been in dynamic problems (Biot, 1956a) and finite deformation (Biot, 1972). Paria (1957-58, 1958a, 1958b, 1966) applied Biot's theories to axisymmetric consolidation of isotropic material under static as well as impulsive loads, trans- verse isotropic semi—infinite mass under normal loads, deformation of viscoelastic body under pressure, spherical isotropic body under pressure, and flow of fluids through deformable bodies. Freudenthal and Spillers (1962) developed theoretical solutions, using Biot's theory, for the infinite layer and the half-space assuming a quasi-static consolidating elastic media. Only a limited number of experimental determinations of coefficients of the equations for Biot's theory have been reported. Biot and Willis (1957) measured the elastic coefficients for sandstone. Fatt (1957, 1959). reported compressibilities of petroleum-bearing sand— stones in the range of 0 to 15,000 psi. He also noted a useful model for sandstone can be developed by a sphere pack composed of a mixture of very hard and very soft spheres. Considerable interest in the use of Biot's theory combined with the finite element has been recently shown in the area of soil consolidation. Most researchers have 12 used variational principles equivalent to the governing equations in Biot's consolidation theory (Sanhu and Wilson, 1969; Yokoo gt_al., 1971; Hwang gt_al., 1971) to solve for pore pressure and settlement under various loading conditions. Hwang gt_al. (1972) used a formulation by the method of weighted residuals. III. DEVELOPMENT OF QUASI-STATIC THEORY OF A GAS-SOLID-LIQUID SYSTEM 3.1 Model Description Consider a medium which is the combination of a deformable solid material, a gas, and a liquid, compress- ible or incompressible. The solid forms the skeleton or the framework of the body and forms a division between two sets of small pores. One set of interconnected pores is filled with a liquid, the the other set contains a gas. The formulation of a mathematical theory of such a three—phase medium starts with the definition of certain relevant variables. Using the differential element pictured in Figure 3.1, we shall define: uZ as the displacements of the solid matrix parallel to the coordinate axes, Ux, Uy, Uz as the displacements of the gas, QX, Qy, Qz as the displacements of the liquid, f as the liquid porosity, defined by f = Vf/VB where Vf is the volume of the liquid and VB the volume of the bulk material within the element, g as the gaseous porosity, defined by g = Vg/VB’ 13 14 LIQUID REGIONS GAS REGIONS Figure 3.1 Differential Element of Medium 15 where Vg is the volume of the gas within the element. The volumes of liquid displaced through unit areas normal to the coordinate directions, X, Y, Z, would be fo, ny, and sz. Similar values for the gas are gUx, gUy, and gUZ. 3.2 Theoretical Development The total stress components of the bulk material, 1 Tij , can be expressed using components as Tij = Cij + Gij (Of + 0g) (3.1) where 0. results from forces applied to the solid part 13 of the body, 5 is the Kronecker delta, and of and 0 ij g result from forces applied to the liquid and gas, respect— ively. Denoting the liquid pressure, pf, and the gas pressure, pg, of = -fpf and 0g = gpg. (3.2) 1 The standard indicial system for a rectangular Cartesian reference frame is employed for the following section: Re- peating the subscripts i, j, or k implies summation, Kron- ecker's delta is denoted by 61-, differentiation with re- spect to space is indicated by subscripts preceded by comma. The subscripts f and g are not to be confused with the summation indices; they indicate liquid and gas, respect- ively. 16 Since the system is considered in equilibrium, pi and pg are assumed constant throughout their respective regions of the body. The strain energy of a porous elastic medium can be defined as the isothermal free energy of the gas-solid- liquid system. Let W denote the strain energy per unit volume. The variation of the strain energy for a volume V bounded by surface S is equal to the virtual work of the surface forces, i.e., If; 6de = £J(rx5ux + fyduy + fzduz + GxGUx (3.3) + Gycuy + czsuz + deQx + ryaoy + anoz) dS where f1, G and F1 are the tractions acting on the 1, solid, gas, and liquid regions of dS respectively. They can be expressed H: l i “ “13 ”J of 51; nj (3.4) 1 = cg 513 nJ '11 II where n is an outward normal to the surface. ,Forces can 3 be expressed in terms of r13, pf, and pg by using (3.1) and (3.2) to obtain 17 Gi = 6ij (-gpg) nj. Introducing the above expressions into (3.3) yields If! 6WdV = £!(fidui + Giaui + Fiin)dS = {!((Tij + 61j(gpg + fpf)nJ-)6ui + 61j(-gpg)5Uinj + dij(-fpf)6Qinj)dS. By defining W. = g(U- - u.) 1 1 1 (3.5) Vi = f6 ll f d1V (ui - Qi ) 20 These variables are now measures of the amount of each substance which has moved in or out of a given element attached to the solid frame. They also represent the increments of fluid and gas contents. The strain energy W must be a function of C. W, and the six strain components, 9 9 Y ) C: W). (3°13) W = W(ex, e y z y, z) Yx) Y 6W must be an exact differential; therefore, 9.11 8w r = r =““ xx 8ex YZ an 3W 3W T = -——'— T = "'—'" yx aey ZX aYy (3.14) ML .311. T = a T = BY zz ez XY z .31! §_W Pf = 3v pg = 8C' These relations lead directly to the formulation of the general stress-strain relations for a gas-solid- liquid medium. Biot (1962) points out several major aspects of this type of derivation where W is the iso- thermal free energy. The stress-strain relations in- clude phenomena which may depend on the physical chemistry of the gas-solid—liquid system; as well as, phenomena 21 which are expressible by means of thermodynamic variables such as interfacial and surface tension effects. 3.3 Linear Stress-Strain Relations Considering the case of an isotropic medium, the strain energy is a function of five variables, the three strain invariants, 11, 12, I3, and the fluid components C and W, i.e. w = w<11, 12, 13, r. W). (3.15) The strain energy is quadratic in form for a linear material, Love (1944). Only the first— and second-order variables are included in the energy expression.' The in- variant terms remaining are I1 = e + e + e 12 eyeZ ezey exey i(Yx Yy Yx ) It is more convenient to use the invariant '=_ ._.. 2+ 2+ 2_ _ _ 12 412 Yx yy yz 4ezey 4ezex 4e e . X y The quadratic form for W now becomes 22 av a He2 + u 15 - 2Ce§ + Mcz - 2Dew + Nv’ + pCW. Sokolnikoff (1956) shows that_the coefficients on any linear terms in W must be zero. Therefore = 2 + 2 + 2 + + 2W H(ex ey ez 2exey + 2exez 2eyez) - 4e e ) 2 2 2 _ _ + “(Yx + y + 7 4e eZ 4ezex x y y Z (3.16) -ZC(ex + ey + ez); —2D(ex + ey + e2)? + Mcz + NYZ + Pew, where H, u, C, D, M, N, P are coefficients which depend on the material properties. Substituting expression (3.16) in the general equations (3.14), we obtain = - + ._ .. W Txx He 2u(ey ez) C; D Tyy = He — 2u(eX + ez) - Cc - DY Tzz = He - 2u(ey + ex) - Cg - DY (3.17) T3’2 = uYX sz - “Yy Txy - “Yz pf = -De + NW + PC/Z p = -Ce + M; + PW/Z. Letting H = A C +'. 2n, C = 3M, 23 (3.17) can be given in the form XX YY ZZ yz Written in r' 1 T XX 1 YY T 22 r yz T ZX H We can define, using (3.18), the bulk modulus of F- A c A c O O 0 —GM 8N l.” Ace + Zuex - A c P Ac+2u e g + A c A+ C .. (1M -BN 2pc Y Zue z T ZX -BNe + -aMe + 211 matrix form —GM -BN aMC - BNW GM; - BNW GM; - BNW = UY Y NW + Pg/Z M; + PW/Z. 0000 O O 0 OT: 0 0 OT: jacketed compressibility as Kc + 2u/3 c and D = 00000 1: SN _aM -aM ’USOOO 2 "O O O O (3.18) (3.19) 24 for this "closed system" in which the fluid pores are sealed. Several relationships between coefficients in (3.16) can be determined by considering the non-negative nature of the strain energy W. Using (3.19), strain energy can be expressed as 2W = Kc e2 - Ce; - 2DeY + Mcz + NW2 + D§W + 2U/3[(ey-ez)2 2 2 (3.20) + (eZ - ex) + (ex - ey) ] lu(Yx Yy vz ) Letting e = C = W = 0 directly implies u 2 0. By putting Yx = yy = yz and ex = ey = ez the strain energy expression reduces to 2w = Kce2 - 2Cec - 2Dew (3.21) + Mcz + NW2 + PcY. By letting W = 0 in (3.19), we obtain the expression 2W = Kce2 - 2C§e + Mcz. (3.22) Since W must be equal to or greater than zero, it can be shown that Kc 2 0 and M 2 0. (3.23) 25 Using the quadratic form for solution with respect to e, we find the discriminate 402 c2 - 4KC(M2;2 — 2W) 2 o for real solutions. We find, for c f 0 This zero, when KCM — c2 2 2 2 o. (3.24) c Similarly, by letting C = 0 in (3.19), we obtain 2w = Kc e2 - 2Dv + NW2. (3.25) expression is never negative if Kc 2 0, N 2 0, and KCN — D2 2 0. (3.26) Again considering (3.20) when neither C or W is the discriminate yields (2c; + 2M)2 - 4KC(M;2 + NW2 + per - 2W) 2 0 solving for e. This relationship can be modified to (KCM - C2)c2 + (KCN — D2)W2 + (ch — 2cn)v; (3.27) — ZKCW 2 o. For real solutions for c, the discriminate of the solution must be non-negative, that is 26 (ch — 2013)“?2 — 4(KCM - C2)[(KCN — 132w2 — 2ch11 2 0. Therefore, considering 'l’ f 0, 2 2 2 2 2 2 (KCM - C )(KCN - D ) - Kc P — 4KCCDP + 4 C D (3.28) 2 0. For expression (3.28) never to be negative 2 2 -KC P 2 0. It was shown from (3.23) that KC is a real number which is greater than or equal to zero, thereby implying 2 P :0. (3.29) Therefore, P=0, for P to be a real constant as assumed. With proof that P = 0, the last two equations of (3.18) can be simplified to pf = -8Ne + NT (3.30) = — + pg dMe M; We can define, using (3.18), the bulk modulus for an "open system” where pg = 0 and pf = 0 as K = A+ 2n/3 (3.31) 27 which is the inverse of the "jacketed” compressibility. The open system would correspond to a jacketed compression test where the pore fluids are allowed to escape freely. Carrying Biot's (1956) analogy between a porous media and a thermoelastic solid, it can be concluded that Ac and A correspond to the adiabatic and isothermal Lamé coefficients for a nonporous medium. Solving equations (3.30) for c and W yields c = 0‘e + p /M g (3.32) W = 89 + pf/N Using C1. M B 2N A - AC — N - M (3.33) (3.18) can be transformed into the form + + = + Txx apg 8pf 2ueX Ae + a + = 2 e + Ae Tyy pg Bpf u y r + up + SD = 2ue + Ae 22 g f Z (3.34) yz uvx zx uvy xy UYZ Z;=OLe+p/M g Y + 8e + pf/N. Written in abbreviated form 28 r + 6 + = 2 e.. + 6.. Ac ij ij(apg Bpf) u 13 13 e + M C a pg/ (3.35) >6 I "' 89 + pf/N. The stress-strain relations can now be written in such a form as to yield the "effective stress," i.e. the total stress in excess of local fluid pressures i. = .. + 6.. + . T13 T1J 13 (Pf Pg) (3 36) OI‘ Tij - Gij ((1 - 00pg + (1 - B)pf) (3.37) = 2n eij + 613. Ae = + M C ae pg/ W = Be + pf/N. The "effective stress” for a fluid—solid system is commonly used in soil mechanics for the study of fluid saturated clays. The bulk modulus expressions for the two types of systems defined as the "closed system," (3.19), and the "open system," (3.31), can be combined with (3.33) to give K - K = A — A= + . (3.38) We obtain by combining (3.24) and (3.26) (Kc - BZN)N + (Kc - azM)M 2 0. Substituting from (3.38) into the above relation KC(N + M - 2MN) + 2KMN +a2M2 + BZNZ 2 0 (3.39) For expression (3.39) never to be negative, we conclude N + M - 2MN 2 0 (3.40) In summary, the non-negative nature of the strain energy yields the following limits to the coefficients of equations (3.18) and those derived from it, (3.41) M 2 O N 2 0 N + M - 2NM 2 0 (K — BZN)N 2 0 (KC — 32M)M 2 o 3.4 Governing Equations for Transient Phenomena The equations for the quasi—static theory of a gas- solid-liquid have been established. This theory shall 30 now be extended to cover the transient phenomena. The equations will describe the distribution of stress, fluid contents, and displacements as a function of time under given loads. It is important to note that the time variable t enters the theory through Darcy's law. There- fore, the transient problem in this case refers to a flow problem. Substitution of (3.35) into the equilibrium equation neglecting body forces as before yields = 2n e.. . — 6ij(apg + Bpf - Aekk),j (3.42) Substitution of the strain displacement relation eij = 5(uihj + uj’i) (3.43) into (3.42) produces .-.+ + ...—O. uu1.JJ (A u)uJ, p J]. g’i - Bpf’i = 0° (3'44) Equations (3.44) are a set of three equations with five unknowns, u, v, w, p pf. Two additional equations are g’ needed to complete the system. These equations can be obtained by introducing Darcy's Law governing the flow of each of the fluids. First consider the flow of the liquid, assuming it to be incompressible. Paria (1966) gives a modified form of Darcy's law 31 8 pf,i ‘ Cf 5? (Q1 “ Bi) (3.45) vffz where Cf =' k for an isotropic media, where vf is the f viscosity of the liquid, hr the permeability of the medium to liquid, and f the liquid porosity. If the solid displacements are zero, “i = 0, (3.45) reduces to the classical form of Darcy's law for an undeformed medium. Combining the divergence of (3.45) with the pf relation of (3.18) and (3.12) yields a -8Nuk,ki + Nf(uk - Qk),ki = Cfgf (Q1 - ui). (3.46) Darcy's equation for the flow of gas in the medium is 3 pg,i - C8 3? (Ui - ui) (3.47) Bag: with Cg== k for an isotropic medium, where vg is the g . viscosity of the gas, kg the gas permeability of the medium, and g the gas porosity. The density, p, for isothermal flow of gases is directly proportional to pg, hence 32%” = & age t g (3.48) a 2 xi Carman (1956). Equation (3.48) can be combined with the last two equations of (3.18) to produce (no sum on j) 2C 2) _ 2 =__E_ ( aMekk + Mt) ,jj g 3t (-OtMekk + Mg) (3.49) where as defined in (3.12) t = g(ui - Ui).J' (3.50) Combining these two equations gives (-GMuk,k + gM(ui - Ui),1)2 .33 2ngM 3 - a _ (no sum on j) Substituting (3.45) and (3.47) into the equilibrium equation (3.44) yields a set of nine differential equations and nine unknown displacements when combined with (3.46) and (3.51). Summarizing, the nine equations are a. ““i.Jj * (A + “) “3,31 ’ “Cg at (”1 ' “1’ .) = o 3 'BCf it. (Qi - ul (3.52) -8Nuk,ki + Nf (Iii - Q1),ij 3 = of??? (Qi " ui) 2 (-o¢Muk,k + gM(ui - Ui),i),qq 8 2C GM 3 = O - o g 2 GEM at ((ul U1).i) + g 3? uk,k. (no sum on q) 33 Written in conventional notation, these equations are uvzfi + (A + u)grad e - an Q% (U -'E) a .- _, (3.53) - BCf 5¥'(Q - u) = 0 -BNVe + va2(fi'-‘6) = cf 5% (5 - G) 2 ¢ 4'2 V (-aMe + gMV(u — U) ZCgaM 39 3 A A = 2gM 3? (V(u - U) + U? . The nine equations governing the transient problem are coupled with each other and hence would have to be solved simultaneously. As Paria (1966) points out for a liquid-solid media, this implies that the flow fields and the elastic field are not merely superposed, but that they react upon each other. This could also be concluded directly from the constitutive equations (3.18). 3.5 Determination of Elastic Coefficients The following series of hypothetical tests are intended to show the physical meaning of the elastic co- efficients d, 8, M, N, and Ac. The coefficient u will be assumed to be determined by standard means. A jacketed test refers to an experiment during which the sample is placed within an impermeable membrane. The pressure of the gas, pg, and the pressure of the liquid, 34 pf, are assumed controllable through tubes which penetrate the membrane and are connected to the appropriate region. See Figure 3.2. 3.5.1 Jacketed Compressibility Test - No Liquid Flow Consider a jacketed specimen subjected to a hydro- static pressure P' such that p = P'/a T = r = T = —P' W = 0. g xx yy zz Define the compressibility under the given test conditions as 6 = --. (3 54) Two relationships are obtained from (3.18) p' = A 5p' + 2n/3 op' — an; ' (3.55) P'/a = aM6P' + Me. (3.56) Assuming P' f 0, combining (3.55) and (3.56) gives 1 - a/a = (AC + 2u/3 - 92M)o. (3.57) From (3.24) Ac + 2u/3 -02M 2 0, therefore, G/a approaches unity as 6 approaches zero which corresponds to an incompressible media. 35 SAMPLE WITHIN AN IMPERMEABLE MEMBRANE Figure 3.2 Specimen for Hypothetical Tests 36 3.5.2 Jacketed Compressibility Test - No Gas Flow Consider that a jacketed specimen subjected to a hydrostatic pressure P' such that pf = P' Txx = Tyy = Tzz = -P' C = 0. Define the compressibility for these conditions as A = PI. (3.58) Two relations can be obtained from (3.18) -P' = -AcAP' - 2u/3 A P' - BNW (3.59) P' = BNAP' + NY. (3.60) Assuming P' f 0 implies l — 8 = (Ac + 2u/3 - BZN)A. (3.61) Equations (3.26) gives Ac + 2u/3 — BZN 2 0 therefore 8 equals one in the limit (incompressible case) and A = 0. Combining the results of the first two tests gives the relation + BZN = + azM. (3.62) 37 3.5.3 Jacketed Compressibility Test - No Liquid Flow, No Gas Pressure Consider a jacketed specimen subjected to a hydro- static pressure P' such that pg = 0 TXX = Tyy = Tzz = —p' w = 0' K = — 75: . (3.63) Equations (3.18) reduce to -P' = —K(Ac + 2u/3)P' - GMC (3.64) and 0 = aMKP' + Mg. (3.65) Assuming P' f O, the combining of (3.64) and (3.65) yields 1/K = Ac + 2u/3 - azM. (3.66) Combining the results of the first and third tests allows K to be expressed in terms of two measurable compressibilities a = a(l - g) (3.67) For a highly compressible gas, the compressibility 6 would be very near that of K implying a small value for a. 38 When pg = 0, (3.32) reduces to and a can be interpreted as the ratio of the change in gas volume to dilitation for the jacketed test. If the gas region within the specimen is connected with the atmosphere by a tube, c would be the amount of gas flow- ing through the tube. Another interpretation of a can be obtained from the first three equations of (3.34), where a is that portion of the gaseous pressure which produces strains. 3.5.4 Jacketed Compressibility Test - No Gas Flow, No Liquid Pressure Consider a jacketed specimen subjected to a hydro- static pressure such that pf=0 TXX=T =1: =_p' C20. Define the compressibility under these test conditions as e O = ———- . 3.68 p, < > Equations (3.18) again reduce to two equations p' = (AC + 2n/3) e p' + BNW (3.59) and O = BNOP' + NW. (3.70) 39 Considering P' f 0, and combining (3.69) and (3.70) produces the equation 1/9 = Ac + 2n/3 — BZN. (3.71) When the result (3.71) is combined with that of Test 3.5.2, a can be expressed in terms of two measured compressibilities m II .._. I I (3.72) The combination of (3.71) and the results of Test 3.5.3 yields the relation 1/9 + BZN = l/K + azM (3.73) which can be shown to be another form of (3.62). Considering the definition of Y for a uniform prorsity, (3.12), and (3.42) can produce the expression W = f(ui, — Q. .) = Be + pf/N i 1 (3.74) = f(eS - ef) where eS and ef are the dilitation of the solid and the liquid regions, respectively. If the liquid is in— compressible, i.e. ef = 0, then (D 'U f—h S —e_ - “r (3'75) It can be seen than pf/N 2 0 and eS/e s 1, therefore 4O 8 2 f. The value of A could be much smaller than the value of 6 for an incompressible liquid and a highly compressed gas, where the liquid is allowed to escape, leaving the solid and gas portions only to support the load. This implies that B is very close to unity for a soft solid material and air. 3.5.5 Jacketed Compressibility Test - No Liquid Flow, No Gas Flow Consider a jacketed Specimen subjected to a hydro- static pressure P' such that TXX = Tyy = Tzz = -P' C = W = 0. Define the compressibility under these test conditions as w = - e/P'. (3.76) Substituting these conditions into (3.18), —P' = — (A0 + 2u/3)wP'. (3.77) Considering P' f 0, implies 1 33 A = _ (3.7821) = A + 2p. (3.78b) 41 Combining (3.66) and (3.71) with (3.78b) gives 1. KIH <3h4 + BZN (3.79) + azM, M and N can now be solved for in terms of measured compressibilities; 'l 1 1 l 5 -2 M = (— - 35);: = (g - K)(1 - z) (3.80) _ 1 1 .1 _ _ 1 i ‘2 N - (— - 6)“? —( - 9M1 - K) (3.81) 3.5.6 Alternate Procedure for Determination of M Consider an unjacketed specimen subjected to a hydro- static pressure P' such that n = _——— . (3.82) Also define a second gas flow related parameter y1= c/P'. (3.83) Equations (3.18) can be used to obtain P' = oan' + Mylp'. For P' f 0, the above relation shows that M can be 42 expressed as 1 M = ;fi_:_VT , (3.84) Biot and Willis (1957) describe a possible experi- mental procedure for determination of y1. A unit volume of porous material is placed within a closed chamber. Gas is injected into the chamber under pressure and the volume of injected fluid is measured. The volume of gas injected per unit pressure will be the sum of the solid— liquid compressibility ¢, the volume of the gas which has entered the pores Y1. and a fixed quantity describing the elastic properties of the chamber and the gas. The differences between the volumes injected with and without the porous material in the chamber will be given by (3.85) AV = 6 + Y1 - Cg where Cg is the gas compressibility. Measurement of the unjacketed compressibility n, and knowledge of the gas compressibility and the gas porosity provide an alternative method for determination of the coefficient M. During the unjacketed test, gas will flow in and out such that the gas pore space and the solid- liquid matrix must undergo the same strain for linearly elastic media. Therefore, the porosity of the gas, g, will not undergo any strain. The dilitation of the gas can be given by 43 =...P'= . eg Cg Ui,i (3 86) where C8 is the compressibility of the gas. The dilitation of the solid-liquid region can be expressed as the sum of the solid dilitation plus the liquid dilitation. If the fluid is considered incompressible and the relative flow, W, is zero, the total dilitation of the solid-liquid region is made up of the dilitation of the solid alone, i.e. e = -nP' = u. S 1,i° Y can now be expressed in terms of compressibilities as 1 __E__g(uii-Uii) _ Y1 - p. p, g(Cg n). (3 87) 3.5.7 Alternate Procedure for Determination of N Consider a jacketed specimen subjected to a hydro— static pressure such that pf = P' T = r = T = -P' C = 0. As in second test, define the compressibility under these conditions as A = -‘———. (3.88) 44 Also define a second liquid flow parameter Y2 = W/P'. (3.89) Equations (3.18) can be used to obtain P' = BNAP' + yzNP'. (3.90) Considering P’ f 0, implies 1 N - BA + Y2 . (3.91) During the jacketed test with pf equal to the applied pressure, the liquid will flow in and out such that the fluid pore space and the solid-gas matrix must undergo the same strain for a linearly elastic media. Therefore, the porosity, f, will not undergo any strain. If the liquid is considered incompressible, e = 0, the sum of the f solid and gas dilitation must be zero, i.e. eS = -eg = —Cg(Apg). (3.92) Obtaining from (3.18) Apg for a change in P' when t = 0, eg can be expressed as e = 0C Me = -aC MAP'. (3.93) g g g Using definition (3.12) for Y, W fes = —— = ——— = -a fMA. 3.94 Y2 P' P' Cg ( ) 45 N can now be expressed in terms of the gas compressibility and liquid porosity 3.6 Closure Equations (3.35), Tij + 6ij(apg + Bpf) = 6ijAe + 2n eij g = ae + pg/M (3.35) Y — Be + pf/N, are similar to the conventional Hooke's law for a material made up of gas, solid, and liquid. The equations include the influence on the continuum of the two fluids which are not in the conventional Hooke's law formulation. Four additional material constants, c, B, M, and N, are necessary for complete use of the formulation. The pressure of the gas can be considered zero under most naturally occurring conditions; therefore, only one parameter, 8, has a significant effect on the stress values. Writing the set of equations in matrix form, 46 r " P s P -r Txx A+2u A A 0 (l 0 (x B ek Tyy A A+2u A 0 0 0 a B ey Tzz A A A+2u 0 0 0 a B eZ TX? 0 0 0 11 0 C) 0 0 nyI 1x2 == 0 0 0 0 11 0 0 0 sz 1&2 0 0 0 0 0 11 0 0 sz C a a a 0 O 0 —l/M 0 —pg :11 .1 “-8 B B 0 0 0 0 -l/MA _-pf.i allows one to see they can be used to replace the conventional relations for development of a finite element model. Associating relative flows with the stress vector and pressures with the strain vector is the most convenient formulation. It is important to note several of the advantages of the present approach. With this approach, it is not necessary to describe the structure of individual cells or the type of interaction which occurs between cells. The model assumes only homogeneity of action in a bulk material sense, not uniformity of cell shape and composition. Measurement of material properties is made on bulk materials not generalized from the action of 47 individual cells. In addition, neither the gas or liquid porosity needs to be determined using this formulation although the effects of both fluids are included in the model. The goVerning equations for the transient phenomena, i.e., the equations governing the medium with fluids in motion, were developed by introduction of Darcy's law of flows for the gas and the liquid (see Equation (3.53)). Additional material pr0perties are necessary for these equations including the fluid porosity terms and terms describing the prOperties of the fluids. These equations imply that the fluid flows and elastic action are not merely superimposed, but are coupled together. Assuming incompressibility of all three materials would uncouple the equations involved in describing the action of the body. IV. EXPERIMENTAL DETERMINATION OF COEFFICIENTS 3 AND M FOR APPLE FLESH A device for applying a hydrostatic pressure to a cylindrical specimen of apple flesh was used to deter— mine the three bulk compressibilities, 6, K, and w, as defined in (3.54), (3.63), and (3.76). The change in volume (volumetric strain) of a cylindrical specimen was related to the applied hydrostatic stress while using different values for the gas pressure, pg. The liquid flow was assumed to be zero in all cases. Equations (3.67) and (3.80) were used to calculate a and M after the three compressibilities were determined. Tests simulating the hypothetical cases discussed in Sections 3.5.2 and 3.5.4 were not attempted because of the diffi- culty involved in measuring and controlling turgor pressure, pf, within the tissue. 4.1 Equipment and Procedure A device develOped by Finney (1963) and used by Brusewitz (1969) was used in this work. The pressure chamber was modified to allow the hydrostatic pressure outside the specimen to vary independently of the air 48 49 pressure within the specimen. The equipment is shown in Figure 4.1 and schematically in Figure 4.2. A carefully cut cylindrical sample 3.81 cm (1.5 in) long with 2.54 cm (1.0 in) diameter was placed on the sample stand and then covered with a highly flexible membrane. The membrane was sealed around the edge of the stand, thereby prohibiting water from entering the sample. Air pressure and air flow for the sample were controlled through the sample stand. Water was used to apply the external pressure of 137.9 kPa (20 psi). The change in water level within a precision bored glass column was used to measure the change in volume within the system. The column had an inside diameter of 5.00 t 0.01 mm (0.1968 : 0.0004 in). Volume change could be measured with a sensitivity of 1 0.00982 cm3 (0.006 in3). All water was boiled to remove the air and then allowed to return to room temperature before each test. It was necessary to remove the water from the system between the calibration run and the actual test to insert the test sample. A calibration curve was determined before each test to account for system characteristics. A steel cylinder was used as the sample during calibration. The system was pressurized to 171.9 kPa (25 psi) before each test to check for leaks. 50 Figure 4.1. Equipment for Measurement of Bulk Compressibility 51 TO AIR SUPPLY PRESSURE REGULATOR I :®+— PRESSURE GUAGE I 4— GRADUATED GLASS COLUMN FILL VALVE :1: +—'— CHAMBER We. SAMPLE A SAMPLE STAND :Q')<——— PRESSURE GUAGE 2 +— PRESSURE REGULATOR 2 4— TO AIR SUPPLY Figure 4.2. Bulk Compression Apparatus 52 Pressure regulator two (Figure 4.2) was removed for determination of K making the air pressure within the sample equal to zero gauge pressure. A plug was in- serted in the sample stand for determination of w. The plug prohibited air from flowing out of the sample. Pressure regulator two (Figure 4.2) was used to vary the internal pressure of the sample for determination of 6. The internal pressure was applied after the external hydrostatic pressure had been applied. Readings were taken for internal pressures of 34.5, 51.7, 68.8, 86.2, and 103.4 kPa (5, 7.5, 10, 12.5, and 15 psi). The two coefficients, a and M, were calculated once the three compressibilities had been determined. Alpha was plotted as a function of the gas pressure using (3.67) a = a(1 - 6/K) where a is the ratio of the gas to the applied hydro- static pressure. The value of a at pg = 0 was deter— mined by projecting the regression 1ine through the vertical axis. M was calculated using (3.81) l 1 M - (w — K)a2 Using a at pg = O. 53 4.2 Results of Experimentation The equipment was found to be unsatisfactory for the determination of the coefficients a and M. The three major problem areas were repeatability of calibration, resolution, and control of applied pressures. The use of the sample stand and membrane to enclose the sample appears to be a satisfactory technique, but the construction of the pressure chamber has some undesirable characteristics. Consistent measurable differences between the two compressibilities K and m were found for the several varieties of apples tested. In all cases, K was found to be larger than w. Unsatisfactory repeatability of calibration made determination of the absolute magnitude of the difference unreliable. The coefficient a appeared to be a nearly linear increasing function of the gas pressure, pg. Calibration of the system was necessary due to the expansion of the equipment and compression of the fluid. The calibration indicated that the system deformation was of the same order of magnitude as that to be measured in the deformed sample. It would be necessary to reduce the volume change of the apparatus and reduce the variability of the volume change for reliable results. 54 4.3 Discussion of Apparatus Limitations The unsatisfactory calibration repeatability of the present system stems from several sources. The large expansion of the apparatus was a major source of error. This expansion was mainly due to the method of attaching the top portion of the apparatus. A rubber gasket, 0.318 cm (0.125 in) thick, was compressed between the upper and lower portions by tightening four bolts. The inability to maintain constant bolt tension created variation in calibration which could not be corrected. The need to remove the water to change samples was a source of error. This disturbance may have created significant variability in the compressibility of the water. Readings for the fluid level in the glass column could be made to i 0.05 mm (0.002 in). That is a volume change of i 0.00932 cm3 (0.0006 ina). A resolution of i 0.005 cm3(0.0003 ina) or less would be necessary to produce satisfactory results. This resolution is equivalent to a diameter of 3 mm (0.118 in) or less. This size of column would be impractical unless the variation of the apparatus itself is reduced. The system was noted to be very sensitive to changes in room temperature during operation. A change of 10C could produce a water volume change of approximately 0.2557 cm3 (0.00867 ina) for the amount of water in the present system. 55 Standard gas pressure regulators were used to control pressures. Difficulty was experienced in holding stable readings with these regulators. Variation in pressure with time at one regulator setting also contributed to the difficulties encountered in determing a and M. Suggestions for an improved apparatus have been in— cluded in Section IV, Suggestions for Further Study. V. FINITE ELEMENT FORMULATION The finite element method is a numerical procedure for solving differential equations and can be used in conjunction with the constitutive equations developed for a gas-solid-liquid medium to calculate stresses in a body of arbitrary shape. The element equations can be derived by minimizing the potential energy of the system. The assemblage of the element equations yields a system of algebraic equations which are solved for the desired quantities. A detailed discussion of the general theory of the finite element method is given in Zienkiewicz (1971), Oden (1972), Desai and Abel (1972), and Martin and Carey (1973). ,The region under consideration is divided into small elements connected at node points for the finite element formulation. The unknowns, stresses, gas pressure, and liquid pressure, are approximated over each element by polynomials. The polynomials for the simplex triangular axisymmetric element (Figure 5.1) are 0t+0Lr+0LZ 1 2 3 C'. II v=3+3r+Bz (5.1) l 2 3 56 57 Figure 5.1. Simplex Triangular Axisymmetric Element 58 ¢1+ ¢2r + 932 6 + 6 r + 6 z. 1 2 3 (5.1) Solving the equations for the coefficients using the nodal values of u, v,-pg,and-pf establishes the inter- polating polynomials for the region. Interpolating polynomials are expressed as a product of shape functions, Ni's, with where and the nodal values. N u + N u + N u 11 22 33 lv1 + sz2 + N3v3 — + _ + - In the above case (5.2) 33) N1(-pf1> + N2(-pf2> + N3(-pf3) functions r + clz)/2A (a1 + b1 (a + bzr + c2z)/2A 2 (a3 + b3r + c32)/2A r223 - r322 r3 - r2 59 etc., in cyclic order, and A is the area of the element. More complicated shape functions can be derived for poly— nomials involving higher-order terms. 5.1 Development of the Variational Equations Variational equations for the present theory were developed following the procedure as outlined by Segerlind (1974) for solid mechanics. Example finite element equations in the following development are for the axi- symmetric case. The column vector of nodal values, {U}2, consists of the nodal displacements, and pressure terms -pf and —pg for the gas-solid-liquid theory. This column vector is T {U} = [ul 2V ,-pgl)pf1! ' ° ' 2 u ,Vna—pgnr—pfn] 1 where n is the number of nodes and u and v the nodal displacements parallel to the R and Z axes, respectively. The total potential energy of the system can be expressed as H = SE - WL (5.4) 2 Standard matrix notation common to finite element work is used in the following development: { }'s denote a column vector,[ ] a matrix, superscript T denotes the transpose of the matrix. 60 where SE is the strain energy and WL the work done by the internal and applied loads (Segerlind, 1974). For a body subdivided into a number of elements, the expression for the total potential can be written in summation form E = 2 (SE - WL ) (5.5) where E is the number of elements and SE(e) and WL(e) are the respective components in each element. The strain energy for an arbitrary differential volume is d(SE ) = é{e}T{T}dV. (5.6) (e) Equation (3.10) simplifies to the above form by letting T {E} II r-H (D (I) (O 2 - 9 - } Pg pf T (5.7) {I} r , §,Y}. II c-a 9.! Assuming that each force and the corresponding displace- ment are linear, the total strain energy is obtained by integrating over the volume of the body giving Ii{e}T{r}dv V(e) SE (e) (5.8) T Jim} [B] TED] [B](UldV. V (e) 61 The material properties matrix,[ D], is defined in {T} = UI]{E}. This matrix is R _jiy _19_ 0 a B l-u l-u R E: 0 a 8 l-u [D] = R O a B (5.9) R(1-2u) symmetric 2(1—u) 0 0 -l/M 0 -1flq with R = E (1 '11) (1 +'IJ)(1 - 2U) The B matrix is defined by the relation {6} = [B]{U} and contains shape functions as well as their derivatives. The general form of[ B] is 62 3N1 8N2 3N T —— —n 0 32 0 32. 0 0 0 32 0 0 a a _N1.0 00.3.1530 00... Nno 00 gr Br Br N N1 0 0 0 PIE 0 0 0 .3 0 0 0 ? I‘ I‘ [B]= (5.10) ’BN 1 3N a a N 3N 7.1003331200...3_En00 2 3r 32 3r 32 '—5f 0 0 N10 0 0 N2 0 0 0 ND 0 0 0 0 N1 0 0 0 N2 0 0 0 N n where n is the number of nodes. The work done by the applied loads can be separated into four components: 1. Work due to concentrated loads, 2. Work resulting from the stress components acting on the outside surface, 3. Work due to movement of liquid across the out- side boundary, and 4. Work due to movement of gas across the outside boundary. The work due to body forces will be neglected in this analysis. 63 The work done by a concentrated force is the value of the force multiplied by the distance through which it acts. Denoting the nodal forces as T {P} = {pr1, p21, 0, 0, ...., Prn, pzn, 0, 0} (5.11) the product becomes T T wc = {U} {P} = {P} {U}. (5.12) The zeros placed in the {P} matrix make the relation compatible with the displacement matrix as defined in relation (5.3). Work done by the stress components onthe surface is described by = + . Wq é(uqr qu)dSl (5 13) l where u and v are the displacements and q1. and q2 are the stress components parallel to the coordinate axes. Equation (5.13) can also be written as T T w = 2K f{U} [M ] {qr} dL . (5.14) A T . . The[ Mi] '8 are used to denote the shape function matrlces. The work done by the flow of the liquid across the boundary can be expressed using (3.6) and (3.5) wf = I v pdeZ = f[ M2]{UlvdSZ. (5.15) 32 32 64 A similar expression for the gas flow is wg = £3 ngd83 = éBEM3]{U}wdS (5.16) 3. The total potential energy is E T T H = z I (MW [8] [DJEBJ{U})dV e=l V _ (e) (e) T T q T ]T -2R g {U} [Mg {qIZ‘}rdL1 - 2n £:U} [M2 {v(e)}rdL2 1 T T T -21r {\EU} [M3] {w(e)}rdL3 - {U} {p}. (5.17) Minimization of the energy expression produces E z ([k 1{U}—{f }) = {p} (5.18) 9:1 (e) (e) where the element stiffness matrix is [k 1 = f [BJTEDJEBJdV (e) V (e) (e) and the element force matrix is _ Tq {f } - 2nd [M11 {q21dL (e) A1 1 T + £042] {v(e)}dL2 2 m3) . T + £043] {w(e) 3 _ 65 5.2 Implementation of the Method An existing axisymmetric finite element computer program and associated subroutines were modified for use with the present theory. Modifications included in- corporating the new constitutive equations and allowing for four unknowns at each node. This program utilized isoparametric elements. Figure 5.2 outlines the structure of the computer solution. The negative values on the main diagonal of the material properties matrix, [0], produced negative terms on the diagonal of the stiffness matrix not found in the conventional formulation. A sub— routine, SIGN, was used to change the sign of a negative terms on the main diagonal of the stiffness matrix before solving the system of equations. A spherical body with a radius of 3.33 cm (1.31 in) was selected for study. The symmetry of the sphere and a load about the vertical axis made it necessary to consider only one quadrant. The body was subdivided into 22 elements. Figure 5.3, using quadratic quadrilaterals (8 nodes/element) and quadratic triangles (6 nodes/element) to fit the curved surfaces of the body and approximate the unknowns within each element. A layer of uniformly thick elements, 0.1524 cm (0.06 in), was placed at the surface to represent the skin of a fruit. Smaller elements 66 STAET READ - Characteristics of Grid + - Set Size of Global Stiffness Matrix - Initialize Global Stiffness Matrix INITIALIZATION + CALL SUBROUTINE L - to calculate element stiffness matrix gBESMAXS .—L + INSERT ELEMENT STIFFNESS INTO GLOBAL STIFFNESS MATRIX + CALL SUBROUTINES Figure 5.2. - to find displacements :BDYVAL - Insert . 811's + SHBJ 4. *IumPB +- .+ + SUWI) + OMIJSUEMIHINE — to calculate stresses BESTAXS + [(XHPUTI + END Outline of Finite Element Computer Program 67 No. of elements = 22 No. of Nodes = 83 Degrees of Freedom = 332 % Width of Band Matrix = 68 )3 3 1; II pp Scale 1 Figure 5.3. Grid for Finite Element Application 68 were used at the top to increase the accuracy since the load was applied parallel to the z-axis. Determination of the element matrix for an iso- parametric element requires numerical integration. Nine integration points were used for each quadrilateral and seven for each triangle. Location of the integration points and weighting factors can be found in Zienkiewicz (1971) or Segerlind (1974). Material pr0perties (Mohsenin, 1970) were chosen to simulate the properties of an apple. An elastic modulus of 13790 kPa (2000 psi) and a Poisson's ratio of 0.30 were used for skin. An elastic modulus of 5171 kPa (750 psi) was used for the parenchyma (flesh). These values for the skin and parenchyma were used in all the numerical problems. Poisson's ratio of the parenchyma was varied to determine its influence on the stress distribution. The following values were used for the other material properties: Parenchyma a = 0.1 M = 6895 kPa (1000 psi) 8 = 0.95 N = 3.10 kPa (0.45 psi) Skin a = 0.0 M = 6895 kPa (1000 psi) 8 = 0.0 N = 3.10 kPa (0.45 psi). 69 The values of the gas pressure, -pg, and the liquid pressure, -pf, were assumed constant throughout the whole body. A pressure of 0 kPa was used for the gas pressure in all cases. 5.3 Results of Finite Element Analysis Three types of numerical problems were run using the quarter sphere grid described in Section 5.2. 5.3.1 Unrestrained Sphere The first problem simulated a body with no skin. One set of material properties was used for the entire body, which was subjected to a turgor pressure and allowed to expand freely. It was determined that a homogeneous body expands due to turgor pressure, -pf, but no stresses are developed within the body. The 3.33 cm (1.31 in) radius sphere with the following material properties E = 5171 kPa (750 psi) u = 0.25 a = 0.10 M = 6895 kPa (1000 psi) 8 = 0.95 N = 3.10 kPa (0.45 psi) enlarged by 6.34 percent in the radial direction due to a turgor pressure of 689.5 kPa (100 psi). 70 5.3.2 Restrained Sphere The second problem simulated the action of an apple with no external restraints. A thin layer with a higher elastic modulus was used to represent the skin. The turgor pressure in this case created stress in both the parenchyma of the apple and the skin. The stress in the parenchyma was compressive and uniform throughout the body. This compressive stress was directly proportional to the turgor pressure for any given Poisson's ratio. The magnitudes of all stress components in both regions, however, were highly dependent on Poisson's ratio. Figure 5.4 shows the hydrostatic stress in the center of the apple for various Poisson's ratios. The turgor pressure was 2758 kPa (400 psi) in each case. The turgor pressure created tension and shear stresses in the skin. Figure 5.5 shows the relationship between maximum principal stress and maximum shear stress in the skin as related to the Poisson's ratio of the parenchyma. Turgor pressure was 2758 kPa (400 psi) in each case. 5.3.3 Flat Plate Contact The third numerical study determined the stress distribution within the spherical body due to compression by a flat plate perpendicular to the vertical axis of 600 500 400 300 STRESS (kPa) 200 100 Figure 5.4. I I I I J I I I I 71 pf = 2758 kPa 0.1 0.2 0.3 0.4 0.5 POISSON'S RATIO Hydrostatic Stress in Parenchyma of Body SKIN STRESSES (kPa) 6000 - 5000 — 4000 -- P 3000 F" p 2000 "‘ 1000 "‘ . ....Je Figure 5.5. 72 p = 2758 kPa Maximum Principal Maximum Shear I 1 I 1 I I I 1 0.1 0.2 0.3 0.4 POISSON'S RATIO Maximum Principal Stress and Shear Stress in the Skin 73 .symmetry. A two-step procedure was necessary to solve the problem. The expansion due to turgor pressure was determined. The set deflection was then subtracted from the expanded equilibrium state and boundary conditions input for the second step. The body was restrained to meet the flat plate deflection conditions during the second run. The body shown in Figure 5.3 was given the follow- ing material properties: Skin Parenchyma E = 13790 kPa (2000 psi) E = 5171 kPa (750 psi) u = 0.30 u = 0.35 a = 0.0 a = 0.10 B = 0.0 B = 0.95 M = 6895 kPa (1000 psi) M = 6895 kPa (1000 psi) N = 3.10 kPa (0.45 psi) N = 3.10 kPa (0.45 psi) The body with no external restraints expanded 0.3138 cm (0.1236 in) for a turgor pressure, -pf, of -2068 kPa (300 psi). The parenchyma developed by a hydrostatic stress of —186.8 kPa (—27.1 psi). A maximum principal stress of 2537 kPa (368 psi) and a maximum shear stress of 1324 kPa (192 psi) deve10ped in the center of the skin. Expansion of the body due to turgor pressure was 74 restrained for the second run to yield a deflection of 0.0762 cm (0.03 in) along the vertical axis of symmetry. The plane of contact to produce a flat surface at that level had a radius of slightly greater than 0.86 cm (0.34 in). The deflection developed sizable maximum shear stresses within the parenchyma, as shown in Figure 5.6. The maximum shear stress of 227.5 kPa (33 psi) appeared approximately 0.318 cm (0.125 in) below the surface along the axis of symmetry. Shear stress in the skin under the area of contact was decreased slightly on the surface and increased slightly on the inner side. Effect on the shear stresses in the skin damped out rapidly beyond the area of contact. The principal stress of the largest absolute mag- nitude was a compressive stress. The stress increased in magnitude from -186.8 kPa (-27.1 psi), as it would be due to turgor pressure alone, as shown in Figure 5.7. 5.4 Closure Fruit splitting generally does not involve external loading but rather variation of the internal conditions. The gas-solid-liquid model developed, used in conjunction with the finite element method, has been shown to be capable of describing the stress distribution within a fruit created by variation of internal conditions such as turgor pressure. Good agreement exists between the 7.) 227.5 lflL4 1Hl3 55J3 216 13.7 6.9 Figure 5.6. Flat Plate Contact - Distribution Curves for T max (kPa) 76 —689 -552 —Zfl3 —241 ~207 —193 Figure 5.7. Flat Plate Contact - Distribution Curves for Principal Stress (kPa) 77 results of the numerical problems and the observations of Considine and Kriedeman (1972) (see page 9, Review of Literature). A fruit without skin expanded freely due to increased liquid (turgor) pressure in both cases. The restraint created by the skin produced stresses in the body with turgor pressure. VI. SUGGESTIONS FOR FURTHER STUDY This work has established a theory for a gas-solid- liquid medium similar to a biological material such as a fruit.‘ The finite element method was shown applicable for solution of numerical problems. Numerous numerical problems can now be approached using the theory in con- Junction with the finite element method. Several areas which should be given high priority are: 1. Study the effect of skin properties; i. e. thickness, elastic modulus, etc., on the stress within the skin and body. 2. Consider the effect of external distributed loads. 3. Consider the case where the body is subjected to large gas pressure, such as in the pressure bomb technique for measuring turgor pressure. 4. Extend the program capabilities to consider the viscoelastic material. 78 79 More accurate determination of material properties is needed before more specific use of the model can be made. Equipment needs to be developed which has satisfactory resolution for determining bulk compress- ibilities. A suggested design is shown in Figure 6.1 for a bulk compressibility apparatus. Pressure sources other than gas regulators should be considered. A column of heavy liquid such as mercury is one possible pressure source. The present theory should be extended to include viscoelastic materials. Biot (1961) has used the corre- spondence principle to deve10p the viscoelastic case for a fluid-solid medium. 80 mapdnaaa< mewafinfiwmohnsoo xfism vmuwmmmsm .H.® madman 52o A0555: 55:, 29:6 \Ifioimmmm hops; . __ k 3) in! H HE) I (11 —|> h‘d rIJ fi sees: a 0858 SE 01:26 \\\/ / x + L 808 1TI>h50pwE AH< L :I. .Alllllwnfib venom :onfloth hammzm MH< DH BIBLIOGRAPHY Akyurt, Akyurt, Biot, M. Biot, M. Biot, M. Biot, M. Biot, M. Biot, M. M. 1969. M., G. 1972. A. 1941. BIBLIOGRAPHY Constitutive Relations for Plant Materials. Ph.D. thesis, Purdue University. L. Zacharich and C. G. Haugh Constitutive Relations for Plant Materials. Trans. of the ASAE., 15(4):766-769. General Theory of Three-Dimensional Con- solidation. Journal of Applied Physics, 12:155-164. A. and F. M. Clinger 1941. Consolidation Settlement of a Soil with an Impervious Top Surface. Journal of Applied Physics, 12:578-581. A. and F. M. Clinger 1942. A. 1955. A. 1956a. 1956b. Bending Settlement of a Slab Resting on a Consolidating Foundation. Journal of Applied Physics, 13:35-40. Theory of Elasticity and Consolidation for a Porous Anisotropic Solid. Journal of Applied Physics, 26:182-185. Theory of Deformation of a Porous Visco- elastic Anisotropic Solid. Journal of Applied Physics, 27:459-467. General Solutions of the Equations of Elasticity and Consolidation for a Porous Material. Journal of Applied Physics, 23:91-96. 81 82 Biot, M. A. and D. G. Willis 1957. Biot, M. A. 1962. Biot, M. A. 1963. Biot, M. A. 1972. Broyer, T. C. 1952. Brusewitz, G. 1969. Burstrom, H. G. 1967. Carman, P. C. 1956. Cleland, R. E. 1971. The Elastic Coefficients of the Theory of Consolidation. Journal of Applied Mechanics, 24:594-601. Mechanics of Deformation and Acoustic Propagation in Porous Media. Journal of Applied Physics, 33:1482-1498. Theory of Stability and Consolidation of a Porous Media Under Initial Stress. Journal of Mathematics and Mechanics, 12:521-541. Theory of Finite Deformation of Porous Solids. Indiana University Mathematics Journal, 12:597-620. 0n Volume Enlargement and Work Expenditures by an Osmotic System in Plants. Physiologia Plantarium, 5:459-469. Consideration of Plant Materials as an Interacting Continuum. Ph.D. thesis, Michigan State University. , I. Uhrstrom and R. Wurscher Growth, Turgor, Water Potential, and Young's Modulus in Pea Internodes. Physiologia Plantarium, 20:213-231. Flow of Gases Through Porous Media. Academic Press, Inc., New York. Cell Wall Extension. Annual Review of Plant Physiology, 22:197-222. Clevenger, J. T. and D. D. Hamann 1968. The Behavior of Apple Skin Under Tensile Loading. Trans. of the ASAE., 11:34-37. Considine, J. A. and P. E. Kriedemann 1972. Fruit Splitting in Grapes: Determination of Critical Turgor Pressure. Australian Journal of Agricultural Research, 23:17-24. Cotner, S. S., 1969. 83 E. E. Burns and P. W. Leeper Pericarp Anatomy of Crack-Resistant and Susceptible Tomato Fruits. Journal of American Society for Horticultural Science, 94:1367. Desai, C. A. and A. F. Abel 1972. Eringen, A. C. 1962. Falk, S., C. H. 1958. Fatt, I. 1957. Fatt, I. 1959. Finney, E. E. 1963. Freudenthal, A. 1962. Frey-Wyssling, 1952. Haines, F. M. 1950. Introduction to the Finite Element Method. Van Nostrand Reinhold Co., New York. Nonlinear Theory of Continuous Media. McGraw-Hill, New York. Hertz and H. 1. Virgin On the relation between Turgor Pressure and Tissue Rigidity. I. Experiments on Resonance Frequency and Tissue Rigidity. Physiologia Plantarium, 11:802-817. Compressibility of Sandstones at Low to Moderate Pressures. Bulletin of the American Association of Petroleum Geolo- gists, 42:1924-1957. The Biot-Willis Elastic Coefficients for a Sandstone. Journal of Applied Mechanics, 26:296-297. The Viscoelastic Behavior of the Potato, Solanum Tuberosum, Under Quasi-Static Loading. Ph.D. thesis, Michigan State University. M. and W. R. Spillers Solutions for the Infinite Layer and the Half-Space for Quasi-State Consolidating Elastic and Viscoelastic Media. Journal of Applied Physics. 33:2661-2668. A. Deformation and Flow in Biological Systems. Interscience Publishers, Inc., New York. The Relation Between Cell Dimensions, Osmotic Pressure and Turgor Pressure. Annals of Botany, (N.S.) 14:385-394. 84 Hwang, C. T., N. R. Morgenstern and D. W. Murray 1971. On solutions of Plane Strain Consolidation Problems by Finite Element Methods. Canadian Geotechnical Journal, 8:109-118. Hwang, C. T., N. R. Morgenstern and D. W. Murray 1972. "Application of the Finite Element Method to Consolidation Problem" in Applications of the Finite Element Method in Geotechnical Engineering. Proceedings of a Symposium at Vicksburg, Miss. Love, A. E. H. 1944. A Treatise on the Mathematical Theory of Elasticity. Dover Publications, New York. Martin, H. C. and G. F. Carey _ 1973. Introduction to Finite Element Analysis. McGraw-Hill, New York. Matzke, E. B. and R. M. Duffy 1955. The Three-Dimensional Shape of Interphase Cells Within the Apical Meristem of Ana- charis Densa. American Journal of Botany, 42:937-945. Mela, M. J. 1967a. Elastic Theory of Cells and Mitochondria in Swelling Process. I. The Membrane Stresses and Modulus of Elasticity of the Egg Cell of Sea Urchin. Biophysical Journal, 7:95-110. Mela, M. J. 1967b. Elastic Theory of Cells and Mitochondria in Swelling. II. Effect of Temperature Upon Modulus of Elasticity of Sea Urchin and of Oyster. Biophysical Journal, 8:83-97. Meynhardt, J. T. 1964. A Histological Study of Berry-Splitting in Some Grape Cultivars. South African Journal of Agricultural Science, 7:707-716. Mohsennin, N. N. 1970. Physical Properties of Plant and Animal Materials. Vol. I. Gordon and Breach Science Publishers, New York. Oden, J. T. 1972. Nilsson, S. B. 1958. Paria, G. 85 Finite Elements of Nonlinear Continua. McGraw-Hill, New York. , C. H. Hertz and S. Falk On the Relation Between Turgor Pressure and Tissue Rigidity. II. Theoretical Cal- culations on Model Systems. Physiologia Plantarium, 11:818-837. 1957-58. Axisymmetric Consolidation for a Porous Paria, G. 1958a. Paria, G. 1958b. Paria, G. 1958c. Paria, G. 1966. Philip, J. R. 1958a. Philip, J. R. 1958b. Elastic Material Containing a Fluid. Journal of Mathematics and Physics, 36:338- 346. Deformation of a Porous Visco-ElaStic Body Containing a Fluid Under Steady Pressure. Bulletin of the Calcutta Mathe- matical Society, 50:71-76. Deformation of Porous Transversely Iso- tropic Elastic Material Containing a Fluid. Bulletin of the Calcutta Mathematical Society, 50:169-179. Deformation of Spherical Bodies of Fluid— Saturated Elastic Material - I. Bulletin of the Calcutta Mathematical Society, 50: 180-188. "Flow of Fluids Through Porous Deformable Solids." In Applied Mechanics Surveys. H. N. Abranson ed. Spartan Books, Washing- ton, D. C. Propagation of Turgor and Other Properties Through Cell Aggregations. Plant Physio- logy, 33:271-274. The Osmotic Cell, Solute Diffusibility, and The Plant Water Economy. Plant Physiology, 33:264-271. 86 Probine, M. C. and R. D. Preston 1962. Cell Growth and the Structure and Mechanical Properties of the Wall of Internodal Cells of Nitella Opaca. Journal of Experimental Botany, 13:111-127. Raschke, K. 1970. Leaf Hydraulic System: Rapid Epidermal and Stomatal Responses to Changes in Water Supply. Science, 167:189-191. Sanhu, R. S. and W. L. Wilson 1969. Finite Element Analysis of Seepage in Elastic Media. Journal Engineering Mechanics Division of ASCE, 95:641-652. Segerlind, L. J. . 1974. Introduction to the Finite Element Method. (Class notes for AE 990, Advanced Topics, Spring 1974), Michigan State University. Slayter, R. O. 1967. Plant-Water Relationships. Academic Press, New York. Sokolnikoff, I. S. 1956. Mathematical Theory of Elasticity. McGraw- Hill Book Company, New York. Tennes, B. R. ' 1973. Physical Properties of Fruit Related to Cracking. Ph.D. thesis, Michigan State University. Yokoo, Y., K. Yamagata and H. Nagaoka 1971. Finite Element Analysis of Consolidation Following Undrained Deformation. Soils and Foundations, 11:37-58. Zienkiewicz, O. C. 1971. The Finite Element Method in Engineering Science. McGraw-Hill, New York. STAT ”11111“ 1111111111 1293 LBR 111111111“ 12