Mathematics and Mechanics of Solids - Semantic Scholar

Article

Diffusion of chemically reacting fluids through nonlinear elastic solids: mixture model and stabilized methods

Mathematics and Mechanics of Solids 1–24 Ó The Author(s) 2014 Reprints and permissions: sagepub.co.uk/journalsPermissions.nav DOI: 10.1177/1081286514544852 mms.sagepub.com

R Hall Materials and Manufacturing Directorate, Air Force Research Laboratory, AFRL/RXBC Bldg 654, 2941 Hobson Way, WPAFB OH, USA H Gajendran and A Masud Department of Civil and Environmental Engineering, University of Illinois at UrbanaChampaign, 205 North Mathews Ave., Urbana, IL, USA Received 8 November 2013; accepted 12 March 2014 Abstract This paper presents a stabilized mixed finite element method for advection-diffusion-reaction phenomena that involve an anisotropic viscous fluid diffusing and chemically reacting with an anisotropic elastic solid. The reactive fluid–solid mixture theory of Hall and Rajagopal (Diffusion of a fluid through an anisotropically chemically reacting thermoelastic body within the context of mixture theory. Math Mech Solid 2012; 17: 131–164) is employed wherein energy and entropy production relations are captured via an equation describing the Lagrange multiplier that results from imposing the constraint of maximum rate of entropy production. The primary partial differential equations are thus reduced to the balance of mass and balance of linear momentum equations for the fluid and the solid, together with an equation for the Lagrange multiplier. Present implementation considers a simplification of the full system of governing equations in the context of isothermal problems, although anisothermal studies are being investigated. The method is applied to problems involving Fickian diffusion, oxidation of PMR-15 polyimide resin, and slurry infiltration, within a one-dimensional finite element context. Results of the oxidation modeling of Tandon et al. (Modeling of oxidative development in PMR-15 resin. Polym Degrad Stab 2006; 91: 1861–1869) are recovered by employing the reaction kinetics model and properties assumed there; the only additional assumed properties are two constants describing coupled chemomechanical and purely chemical dissipation, and standard values for viscosity of air and PMR-15 stiffness properties. The present model provides the individual constituent kinematic and kinetic behaviors, thus adding rich detail to the interpretation of the process in comparison to the original treatment. The last problem considered is slurry infiltration that demonstrates the applicability of the model to account for the imposed mass deposition process and consequent effects on the kinematic and kinetic behaviors of the constituents. Keywords Mixture theory, oxidation, slurry infiltration, stabilized method, variational multiscale method, PMR-15 resin

Dedicated to Kumbakonam R Rajagopal Corresponding author: A Masud, Professor of Mechanics and Structures. Department of Civil and Environmental Engineering, University of Illinois at Urbana-Champaign, 205 North Mathews Ave., Urbana, IL 61801-2352, USA. Email: [email protected]

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

2

Mathematics and Mechanics of Solids

1. Introduction In this paper, a stabilized mixed finite element method is presented for the diffusion of a chemically reacting fluid through a nonlinear elastic solid using a mixture theory based model. For a detailed introduction to mixture theory, interested readers are referred to comprehensive review articles by Atkin and Craine [1], Green and Naghdi [2,3], and the book by Rajagopal and Tao [4]. Mixture theory ideas have been used to model various phenomena, such as classical viscoelasticity [5], swelling of polymers [6], thermo-oxidative degradation of polymer composites [7,8], and growth of biological materials [9] and crystallization of polymers [10]. Malek and Rajagopal [11] proposed that processes for fluid mixtures are governed by the maximization of the rate of dissipation constraint. Karra [7] developed a mixture theory model and its constitutive relations based on this constraint for diffusion of a fluid through a viscoelastic solid. Karra and Rajagopal [8] also developed a mixture theory model for degradation of polyimides due to oxidation. A limitation of their model is that it cannot predict the oxidation layer thickness growth. Hall and Rajagopal [12] proposed a mixture theory model for diffusion of chemically reacting fluid through an anisotropic solid. Their model is based on the maximization of the rate of entropy production constraint, considering anisotropic effective reaction rates and the limits of diffusion-dominated (diffusion of the reactants is far more rapid than the reaction) and reaction-dominated (the reaction is far more rapid than the diffusion of the reactants) processes. This model in general can be applied to a variety of processes that involve directionality of flow, directionality of the reaction process, and that of the solid medium. Such processes arise in, as examples, the curing of composites using vacuum-assisted resin transfer molding (VARTM), in the prediction of oxidation layer growth in composites, and in slurry infiltration (SI) in the manufacturing of composites. The mixture theory model combines the composite constituent behaviors in an effective medium sense, reducing the computational cost of modeling chemically reacting multi-constituent mixtures, while retaining information involving the kinematic and kinetic responses of the individual constituents. The effective medium and individual constituent behaviors are each constrained to mutually satisfy the balance principles of mechanics. In this work we employ Variational Multiscale (VMS) framework [13–16] to develop a stabilized mixed finite element formulation involving the balance of mass equation for the fluid that is written in an Arbitrary Lagrangian Eulerian (ALE) form. The underlying idea in VMS is an additive decomposition of the solution field into coarse- and fine-scale components, which in the present context leads to a multiscale decomposition of the fluid density field that results in two coupled nonlinear problems termed as the coarse-scale and the fine-scale sub-problems. The space of coarse scales is identified with the standard finite element functional spaces, while the space of fine scales can contain various finitedimensional approximations that are subject to the condition that the spaces are linearly independent. The fine-scale solution is extracted from the nonlinear fine-scale sub-problem, and it is then variationally embedded in the coarse-scale equation, leading to a formulation that is expressed entirely in terms of the coarse scales. Consequently, the resulting stabilized form does not have an explicit appearance of the fine-scale density field of the fluid. Rather, the effects of fine scales are represented via the additional residual-based terms that in fact add to the stability of the numerical method. One of the applications of interest in the present work is the oxidation of polymer matrix composites (PMCs). Tandon et al. [17] conducted experiments to study oxidation processes in a high-temperature polyimide resin used in aerospace composites, and developed an oxidation reaction rate model that conforms to the observed experimental data. In this work, we implement this oxidation model in the context of mixture theory. Schoeppner et al. [18] and Varghese et al. [19] and Varghese and Whitcomb [20] developed finite element algorithms for the diffusion-reaction equation to model the oxidation in PMR15 resin and PMCs. In their work, the fibers and matrix were modeled in a discrete sense and thus their algorithm was computationally intensive. Varghese et al. [19] proposed an adaptive mesh strategy and decoupled subdomain strategy to reduce the computational cost of their algorithm. Their adaptive mesh strategy requires a prior knowledge of oxidation layer growth to constrain the unoxidized region, thus reducing the number of unknowns in the problem. The outline of the paper is as follows. In Section 2, we present the governing equations and the constitutive relations derived from the mixture theory for a chemically reacting fluid diffusing through a nonlinear elastic solid in a general three-dimensional context. The modeling assumptions and the onedimensional form of the general mixture theory are presented in Section 3. In Section 4, we present the

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

Hall et al.

3

weak form of the mixture theory governing equations and develop the VMS-based stabilized method. Section 5 presents the finite element results of the mixture theory for Fick’s diffusion problem, in the context of matching an analytical solution for demonstration of accuracy and stability of the numerical approach; oxidation of PMR-15 resin; and SI in PMCs. Conclusions are drawn in Section 6.

2. General mixture theory In this section, we first present the mixture theory-based model for diffusion of an anisotropic nonlinear viscoelastic fluid through an anisotropic elastic solid with mutual chemical reaction, as proposed by Hall and Rajagopal [12]. A basic assumption in the mixture theory is that the constituents of the mixture cooccupy the domain and as the mixture deforms, these co-existing continua deform with respect to each other. A set of appropriate constitutive relations that are based on the constraint of maximum rate of entropy production are also presented in [12]. In the present work, we consider the constitutive relations associated with unconstrained constituent volumes. Detailed derivation is available in [12]. The equations of mass and linear momentum balance for the diffusion of a chemically reacting fluid through a finitely deforming thermoelastic solid are given as follows [12]: Balance of mass :

Da r a ∂ra + ra divva = + div(ra va ) = ma dt ∂t

Balance of linear momentum: ra

Da v a = div(Ta )T + ra b + Ia dt

ð1Þ ð2Þ

where ra is the mass concentration and ma is the rate of mass transferred by chemical reaction, to constituent a, per unit mixture volume; va is the velocity of constituent a and Ta is its partial Cauchy stress, while Ia and b are the interactive force per unit mixture volume on constituent a and the body force per unit mass. The balance of energy and assumption of maximized rate of entropy production, together with Newton’s third law, lead to the following relations for the partial stresses on the solid and fluid, T s and T f ; the interactive force I f on the fluid, the constituent entropy ha , and the rate of fluid mass conversion, m f , all per unit mixture volume; and the heat flux q, per unit mixture area:  s f  rs rf rr I f = g f rr f  g s rrs  r (c f  cs ) r r r ð3Þ s f rr v f f s f f s s (ru) (h  h )  m (v  v )  mA (v  v ) r T s = rF s (

∂c T rf f s s (c  cs ))I s )  r (g + ∂F r

ð4Þ

rs s (c  c f ))I + mAL  D f r

ð5Þ

T f =  r f (g f +

ha = 

∂ca m a _  cu u r ∂u

q rs r f f =  mlru + (h  hs )(v f  vs ) u r   1 1 f f f s s f s (g  g )  (v  v )  (v  v ) m = mcm 2

ð6Þ ð7Þ ð8Þ

where the chemical potential g a of constituent a is defined through g a [r

∂c ∂ra

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

ð9Þ

4

Mathematics and Mechanics of Solids

r, c, and u are the mixture density, mixture Helmholtz energy, and temperature, respectively, while ca are the constituent Helmholtz energies; material parameters cau and cm are respectively associated with the constituent entropies, and with mass transfer, while l is the mixture thermal conductivity tensor; F s is the solid deformation gradient, Av and AL are drag and viscosity coefficient tensors, and D f is the fluid rate of deformation tensor. The rate of mass transfer to the fluid, m f , is determined in coordination with the orientation average _ Because of the presence of only two constituents, the mass balance proof the rate of reaction tensor G. vides that the rate of mass converted to the solid is the one lost from the fluid: ms = m f

ð10Þ

In the diffusion-dominated approximation (diffusion of the reactants is far more rapid than the reac_ tion), the operator Pfn, Xs , tg provides the directional solid mass conversion rate in the direction 2n, per unit mixture volume, such that: 1 m = 4p s

4p ð

_ Pfn½a, Xs , tgda

ð11Þ

a=0

where n is the outward unit normal, Xs is the reference coordinate of the solid, a is the solid angle, and a _ second-order representation is assumed for the operator Pfn, X s , tg: _ Pfn, Xs , tg’n  G_ ½Xs , t n

ð12Þ

with the tensor G ½Xs , t =

ðt

G_ ½Xs , t dt

ð13Þ

0

thus providing an anisotropic measure of the extent of reaction of the solid. Employing in the present work the Lagrangian solid strain measure Es and referring G to material coordinates, the Lagrange multiplier arising from the constraint of maximized rate of entropy production is given by, in the general case, cf. [12]:   1 _0 0 0 s 0 rs r f f s f s f f s f s  ru  G G G E K (h  h )(v  v )  m (v  v )  (v  v ) 1 4 IJ KL MN OP IJKLMNOP r ð14Þ m= + ( f ) 2 D  AL  D f + cu u_ 2 + (v f  vs )  Av (v f  vs ) + G_ 0  A0G  G_ 0 + ru  lru + cm (m f )2 1 s 0 + G_ 0AB G0CD G0EF EGH KABCDEFGH 2 0 where KIJKLMNOP is a tensor that couples the mechanical and chemically influenced attributes of the model, in a way that is compatible with the results of the maximization of the rate of entropy production as described in [12]. Because ha and m f depend on m, Equation (14) is a cubic equation in m. To obtain a single-valued relation for m, the following approximations are made:

1. 2. 3.

We assume that the attributes of the Helmholtz free energy functions of the constituents and the mixture can be represented in terms of suitably condensed forms, cs = c f = c, hs = hf = h; Slow diffusion permits neglect of the squared relative kinetic energy terms ((v f  vs )  (v f  vs ))2 , which are assumed also negligible relative to the drag force; We assume that the reaction is near enough to equilibrium to neglect the squared difference in the chemical potentials of the constituents, and the product of the chemical potential difference with the relative kinetic energy.

The Lagrange multiplier is thus reduced to the following single-valued function:

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

Hall et al.

5

m=

0 1 _0 0 s 0 1 4 GIJ GKL GMN EOP KIJKLMNOP  + f 2 D  AL  D f + cu u_ 2 + (v f  vs )  Av (v f  vs ) + s 0 KABCDEFGH G_ 0  A0G  G_ 0 + ru  lru + 12 G_ 0AB G0CD G0EF EGH

ð15Þ

It can be noted that the tensor K0 will have mostly zero-valued components. If reaction processes such as oxidation are considered, in which the reaction is several times faster in the fiber direction than the transverse directions, thus promoting a unidirectional reaction assumption, and assuming transversely isotropic coupling to the strains, the term involving K0 reduces to the following expression, involving four independent constants:   s 0 s s s s s s KIJKLMNOP = G_ 011 (G011 )2 K10 E11 + K20 (E22 + E33 ) + K30 (E12 + E31 ) + K40 E23 ð16Þ G_ 0IJ G0KL G0MN EOP In the present work, the influence of the energy and entropy production relations is retained through the presence of the Lagrange multiplier, which is obtained via invoking the constraint of maximized rate of entropy production. The equations explicitly retained are the constituent momentum balances and the mass balance equation, which can be considered most strongly enforced. In accordance with the present study being isothermal, the traditional heat capacity measures of the constituents are lost through the assumption above that the constituent entropy functions can be replaced by an overall entropy function. In general for anisothermal processes, the Helmholtz and entropy functions of each constituent would be retained. It is interesting to note, however, that the present system of equations incorporates the rate of temperature in combination with a non-traditional overall material property cu (the density average of the cau properties), which may provide a simplified approach to accounting for a class of homogenized anisothermal effects. The present paper, however, considers only isothermal conditions.

3. One-dimensional mixture theory Consider a one-dimensional mixture domain O of length L with boundary ∂O = fxjx 2 f0, Lgg. The governing equations for the one-dimensional case under isothermal conditions are as follows: Balance of mass :

∂ra ∂ra a a ∂va1 + v +r  ma = 0 ∂t ∂x 1 ∂x

Balance of linear momentum :

a ∂T11 Dva + ra b1 + I1a  ra 1 = 0 ∂x Dt

ð17Þ ð18Þ

The corresponding stresses and interactive force on the constituents can be written as follows: s T11 = rF11

∂c ∂c r f f (c  cs ))  rs (r s + ∂F11 ∂r r

f ∂c rs s L ∂v1 f +  c )) + mA (c ∂r f r ∂x   ∂c s ∂r f ∂c f ∂rs ∂ rs r f f f s r   (c  c )  I1 = f r ∂r ∂x ∂rs ∂x ∂x r ∂u rs r f f (h  hs )  m f (v1f  vs1 )  mAv (v1f  vs1 ) ∂x r f T11 =  r f (r

ð19Þ

ð20Þ

ð21Þ

We consider the following Helmholtz free energy function that corresponds to the one-dimensional representation of a transversely isotropic thermoelastic solid permeated by a chemically reacting Newtonian fluid:

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

6

Mathematics and Mechanics of Solids

 o cs1 u 1 n f s 2 s Rur + k2f r f c = A + (B + c )(u  u )  (u  u )  c2 u ln s + u rT 2   s r 1 1 s 1 s 2 l + msT + as + 2(msL  msT ) + bs (E11 + ) +L s 2 r rT 2 ð n o m  0G dG0  1 Es + 2G_ 0 A L= (G011 )2 K 11 11 11 2 s

s

s

s

ð22Þ

ð23Þ

where L describes the coupling between the solid strain and the extent of reaction, consistent with the developments of [12]; ls , as , msL , msT , bs are the transversely isotropic material constants, which in one dimension reduce to the elastic moduli of the solid. rsT , rT are the true solid density and the true mix is the ratio of the universal gas constant to the molecular weight of the fluid. ture density, respectively. R 0 0G  0G are defined for convenient manipulations involving L.  K1 =  rK1 and A =  rA Remark: For the case of the slurry deposition process that is presented in Section 5.3, G011 represents the extent of material deposition. For this case, the term L provides coupling between the solid strain and the extent of deposition of the suspended particles. We assume that this deposition function G011 is in fact a function of the volume fraction of particles, which is considered a process parameter.

The one-dimensional representation of the Lagrange multiplier m is given as 0 2  1 _0 s 1 4 G11 (G11 ) r K1 E11 m=  f 2 AL ( ∂v1 )2 + c u_ 2 + Av (v f  vs )2  rA  0G (G_ 0 )2 + I11 ( ∂u )2  1 G_ 0 (G0 )2 rK  1 Es u 11 1 11 1 ∂x ∂x 2 11 11

ð24Þ

Also, from mass balance law and Newton’s third law we see that the solid and the fluid interactive forces have the following relationship: I1s + ms vs1 =  (I1f + m f vs1 )

ð25Þ

3.1 Modeling assumptions and methodology In mixture theory where both solid and fluid co-occupy the domain and fluid moves relative to the deforming solid, it is natural to write the fluid balance laws in an ALE framework [21–23]. For the class of problems considered in this work, the inertial effects on the solid are assumed to be negligible. Based on these modeling assumptions, the balance laws Equations (17) and (18) can be rewritten as follows: ∂rs ∂rs s ∂vs + v1 + rs 1  ms = 0 ∂t ∂x ∂x

ð26Þ

s ∂T11 + rs b1 + I1s = 0 ∂x

ð27Þ

f ∂r f

∂r f f m f ∂v1 (v  mf = 0 +  v ) + r 1 ∂t Y ∂x 1 ∂x

f f

∂T11 ∂v ∂v1f

) + r f b1 + I1f  r f 1 r f (v1f  vm =0 1 ∂x ∂t

∂x

ð28Þ

ð29Þ

Y

∂(  )

where represents the time derivative in the ALE frame [22,23] and vm 1 is the fluid mesh velocity. ∂t Y It is important to note that as the solid domain deforms, the Lagrangian mesh that is tied to material s s points deforms together with it. Consequently, the mesh velocity vm 1 is set equal to v1 , where v1 is the velocity of the solid domain. Accordingly, the constitutive relations can be rewritten as

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

Hall et al.

7

s T11 = rF11

∂c ∂c  rs r s ∂F11 ∂r

f T11 =  rfr

ð30Þ

f ∂c L ∂v + mA ∂r f ∂x

ð31Þ

I1f =  m f (v1f  vs1 )  mAv (v1f  vs1 )

ð32Þ

m f = G_ 011

ð33Þ

Remark: In [12] an expression for the rate of mass conversion for fluid m f is derived via maximization of the rate of dissipation constraint. However, in the present work we prescribe an oxidation rate given in [17] that is developed based on physical measurements. Likewise, in the SI model we prescribe a rate of particle deposition as is given in [24]. Because of these postulated rates, the physics involved in the consistent derivation of mass conversion given in [12] is circumvented.

4. Weak form and development of the stabilized method The initial conditions for the density and velocity fields of the two constituents and the displacement field of the solid are ra (x, 0) = ra0 ;

va1 (x, 0) = va0 ;

us1 (x, 0) = us0

8x 2 O

ð34Þ

The boundary ∂O admits decomposition into ∂Og and ∂Oh , where ∂Og \ ∂Oh = f, and we denote the unit normal to the boundary ∂O by n1 . The boundary conditions for the problem are r f = rf0 v1f = vf0 us1 = us0 f T11 n1 = t0f s n1 = t0s T11

on on on on on

f

∂Org 3 0, T½ ∂Ogf 3 0, T½ ∂Osg 3 0, T½ ∂Ofh 3 0, T½ ∂Osh 3 0, T½

ð35Þ

where rf0 , vf0 are the prescribed fluid density and velocity, and us0 is the prescribed solid displacement. t0f and t0s represent the prescribed fluid and solid boundary tractions, respectively. Let g a and wa1 denote the weighting functions for the balance of mass and linear momentum for the corresponding constituent, respectively. The appropriate spaces for these weighting functions are n o a V = ga jg a 2 H 1 (O), ga = 0 on ∂Org ð36Þ n

Q = wa1 wa1 2 H 1 (O), wa1 = 0

on ∂Oag

o

ð37Þ

The corresponding trial solution spaces for the fluid and solid density, fluid velocity, and solid displacement are n o a a S r = ra (, t)jra (, t) 2 H 1 (O), ra (, t) = ra0 on ∂Org 3 0, T ½ ð38Þ

n o

S f = v1f (, t) v 1f (, t) 2 H 1 (O), v 1f (, t) = v 0f on ∂O fg 3 0, T ½

ð39Þ

n

S s = us1 (, t) us1 (, t) 2 H 1 (O), us1 (, t) = us0

ð40Þ

o on ∂Osg 3 0, T ½

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

8

Mathematics and Mechanics of Solids

The weak form of governing equations for the solid–fluid system can be stated as follows. For all constia tuents a 2 fs, f g, 8t 2 (0, T ), 8g a 2 V, and 8wa1 2 Q, solve ra 2 S r , v1f 2 S f , and us1 2 S s such that the following system holds. Weak form of equations for the fluid: !

  f f

f f f ∂r

f s ∂r f f ∂v1  (g f , m f ) = 0 g , ð41Þ + g , (v1  v1 ) + g ,r ∂t Y ∂x ∂x

! ! f

∂wf1 f ∂v

 (wf1 , r f b1 )  (w1f , I1f ) + w1f , r f 1 + ,T ∂x 11 ∂t

Y ! ð42Þ f f f f f f s ∂v1  (w1 , T11 n1 )∂Of = 0 w1 , r (v1  v1 ) h ∂x Weak form of equations for the solid:    s s s s ∂r s s ∂v1 s s ∂r g, + g ,r + g , v1  (g s , ms ) = 0 ∂t ∂x ∂x  s ∂w1 s s ,T n1 )∂Osh = 0  (ws1 , rs b1 )  (ws1 , I1s )  (ws1 , T11 ∂x 11 R where ð, Þ = ðÞdO is the L2 ðOÞ inner product.

ð43Þ ð44Þ

O

4.1 Fluid sub-system: residual-based stabilization Our objective is to model the diffusion of a chemically reacting fluid through a nonlinear elastic solid, a phenomenon that is observed in the process modeling of composites, oxidation of resin/composites, and SI in porous media, to name a few. In the modeling of these processes, fluid mass concentration is invariably specified at the inlet boundary. Since the strong form of mass balance of fluid given in Equation (17) is a first-order hyperbolic equation, any specified mass concentration boundary condition at the inlet that is different from the initial condition results in a discontinuous fluid concentration field. This discontinuity introduces spurious oscillations in the computed solution right at the beginning of the nonlinear iterative process that can lead to non-convergent and therefore non-physical solutions. To address this issue, we consider the weak form of the balance of mass equation for the fluid that is written in an ALE form. We employ VMS ideas [13–16] and develop a stabilized weak form for Equation (41). Underlying idea of VMS is an additive decomposition of the solution field into coarseand fine-scale components as given below: r f = r^ f + r~ f

ð45Þ

g f = g^ f + g~ f

ð46Þ

where r^f , r~f represents the coarse-scale and fine-scale components of the density field and g^ f , g~ f represents the coarse-scale and fine-scale counterpart of the weighting function, respectively. Various scale separations of rf are possible in Equation (45). However, they are subject to the restriction imposed by the stability of the formulation that requires the spaces for the coarse-scale and fine-scale functions to be linearly independent. In the development presented here, the spaces of coarse-scale weighting functions are identified with the standard finite element spaces, while the fine-scale weighting functions can contain various finite-dimensional approximations, for example bubble functions or p-refinements or higher order Non-Uniform Rational B-Spline (NURBS) functions. Substituting Equations (45) and (46) into Equation (41) and employing the linearity of the weighting function slot in Equation (41), we obtain the coarse-scale problem and the fine-scale problem as given in Equations (47) and (48), respectively:

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

Hall et al.

9

! 

∂(^ r f + r~ f )

r f + r~ f )v1f r f + r~ f ) f ∂(^ f s ∂(^  g^ , v1 = g^ ,  (^ gf , mf ) = 0

+ g^ , ∂t ∂x ∂x Y

ð47Þ

! 

f f f

f f ∂(^ r + r ~ ) ∂(^ r + r ~ )v r f + r~ f ) r f f f s ∂(^ 1

~ W = g~ ,  (~ gf , mf ) = 0  g~ , v1

+ g~ , ∂t ∂x ∂x Y

ð48Þ

^r f

W

f



f



It is important to note that both systems are nonlinear, and are also fully coupled in terms of the scales. The key idea at this point is to solve the fine-scale problem Equation (48) locally, using analytical methods or numerical methods, and extract the fine-scale component, r~f . This can then be substituted into the corresponding coarse-scale problem given in Equation (47), thereby eliminating the fine scales, yet modeling their effects. 4.1.1 Solution of the fine-scale problem. We segregate the terms into coarse-scale and fine-scale terms and group

all the terms containing a coarse-scale density field: ! 

 f

f f ∂~ r ∂~ r v rf f r f f f s ∂~ 1

+ g~ , ^ =0 ~  g W~ = g~ , , v + (~ g f , R) 1 ∂t Y ∂x ∂x

ð49Þ

^ is the residual of the Euler–Lagrange equations of the coarse scales over element interiors and is where R given as

f

∂^ r r f v1f rf s ∂^

+ ∂^ ^=  v  m f (^ R rf ) ð50Þ 1 ∂t Y ∂x ∂x In obtaining the above form of the fine-scale problem, we have assumed that the fluid mass conversion rate is a function of the coarse-scale fluid density field only, m f (^ r f , r~ f )’m f (^ r f ). To reduce the complexity of the fine-scale problem and also to reduce the computational cost for evaluating the fine-scale solution field, we assume that the fine-scale field vanishes at the element boundaries: g~ f = 0, r~f = 0

on Ge

ð51Þ

Remark: The assumption that fine scales vanish at the inter-element boundaries helps in keeping the presentation of the ideas simple and concise. Relaxing this assumption in fact leads to a more general framework. This, however, requires Lagrange multipliers to enforce the continuity of the fine-scale fields across inter-element boundaries. It is important to note that Lagrange multipliers can be accommodated in the present hierarchical framework as well.

Using the backward Euler time integration scheme and assuming that the fine-scale fluid density field at the beginning of a time step is zero, r~nf = 0, we can obtain the time discretized form of Equation (49) as given below: ! 

 f

f f r ~ ∂~ r v rf f r f f f s ∂~ 1

~ ^ =0 W = g~ , + g~ ,  g~ , v1 ð52Þ + (~ g f , R) Dt Y ∂x ∂x The fine-scale fields are represented by bubble functions within each element and are given as g~nf + 1 = be2 gnf + 1 ,

r~nf + 1 = be1 rnf + 1

ð53Þ

where be1 , be2 are bubble functions and gnf + 1 , rnf + 1 are the coefficients associated with the fine-scale fields over the element; examples are shown schematically in Figure 1. Substituting Equation (53) into Equation (52), and following along Masud and Khurram [15] and Calderer and Masud [23], we can obtain the fine-scale density field via solution of Equation (52) as follows:

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

10

Mathematics and Mechanics of Solids

Figure 1. Schematic representation of quadratic and linear bubbles.

^ r~ f = tR

ð54Þ

 ^ is the mean value of the residual of the Euler–Lagrange equations of the coarse scales for where R Equation (47). The stabilization parameter, t, is given as be1 (be2 , 1)Oe

t= (be2 ,

1 e Dt b1 )Oe

+ (be2 ,

∂v1f ∂x

be1 )Oe + (be2 , v1f

∂be1 ∂x )Oe

 (be2 , vs1

∂be1 ∂x )Oe

ð55Þ

We now substitute the fine-scale solution given in Equation (54) into the coarse-scale problem, Equation (47): ! 

 f

f f ∂^ r ∂^ r v rf f r f f f s ∂^ 1

^  g^ , v1 W = g^ , + g^ ,  (^ g f , mf ) ∂t Y ∂x ∂x ð56Þ   f  1 ∂vs1 ^ ∂^ g f  f s ^ =0 + (^ g , + tR)  , v1  v1 tR Dt ∂x ∂x Equation (56) represents the modified coarse-scale problem with the fine-scale effects embedded implicitly via the coarse-scale residual terms. The first four integral terms in Equation (56) correspond to the standard Galerkin method for the balance of mass for the fluid. The last two terms in Equation (56) have appeared because of the fine-scale density field. It is important to note that the fine-scale density does not explicitly appear in Equation (56); rather, the fine-scale effects are implicitly reflected in this form via the modeling terms. Equations (56), (42), (43), and (44) are linearized and solved simultaneously for the density and velocity fields of the fluid, and density and displacement fields of the solid using Newton–Raphson solution procedure. This coupled system of equations is discretized-in-time using the backward Euler scheme, while linear and quadratic Lagrange elements with equal-order fields are employed in the spatial dimension. The resulting stiffness matrix for the full system is non-symmetric.

5. Numerical results We present three test cases that investigate the stability and accuracy of the numerical method developed for the mixture theory model described in Section 4. In Section 5.1, we solve a reduced mixture model that is equivalent to the Fick’s diffusion problem. A system comprising a first-order hyperbolic equation and an algebraic equation is solved and the results are compared with the exact solution. Section 5.2 presents the oxidation problem of PMR-15 resin wherein a full system of mixture theory equations is solved and the results are compared with the experimental and numerical results reported by Tandon et al. [17]. Section 5.3 models the SI process that is involved in the manufacturing of composites, and a parametric study of the reduction in the porosity of the solid as a function of slurry particle fraction and initial solid porosity is presented. The present developments have been carried out in the context of the one-dimensional finite element method, and extension to the three-dimensional case will be pursued in a subsequent paper.

5.1 Fick’s diffusion problem In this section we employ a reduced mixture model to solve Fick’s diffusion problem. The transient Fick’s diffusion equation can be derived from the mixture theory balance laws, Equations (17) and (18),

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

Hall et al.

11

based on the following simplifications: (a) solid is assumed to be rigid; (b) fluid is assumed ideal; (c) fluid inertial effects are neglected; and (d) fluid is assumed non-reactive. The constitutive relations for an ideal fluid and the interactive force between the fluid and rigid solid can be given as [25] f  T11 =  r f Ru

ð57Þ

I1f =  Av r f v1f

ð58Þ

where Av is the drag coefficient. The governing Equations (17) and (18) can be reduced to the following system of equations:



∂r f ∂(r f v1f ) =0 + ∂x ∂t

ð59Þ

∂rf  Ru  Av rf v1f = 0 ∂x

ð60Þ

Since the coupled system of Equations (59) and (60) serves as a reduced-order model for the mixture theory, we solve this first-order system to investigate the underlying characteristics of the mixture model wherein the conservation of mass equation for the fluid is hyperbolic. The diffusivity of the solid can be written in terms of the drag coefficient of the solid as D=

 Ru Av

ð61Þ

The derivation of Equation (61) is provided in the Appendix. Remark: Solving for fluid velocity from Equation (60) and substituting back into Equation (59), one can obtain Fick’s diffusion equation. Since our full mixture model results in a first-order system, in this work we have opted to solve the reduced system also in its first-order form to help serve as a test case for evaluating our numerical method.

The unknown fields in this problem are the fluid concentration and fluid velocity that are solved with zero initial conditions. The one-dimensional domain of length 0.001 m is exposed to air at the left end of the domain where fluid concentration is assigned a value of 22.8863E-3 kg/m3 and fluid velocity is  and drag coefficient Av are constrained to be zero at the right end of the domain. The gas constant R 21 assigned values of 286.987 J/kg-K and 1.63E17 s , respectively. We employ the backward Euler scheme for time integration and run the problem for a total time of 30,000 seconds. A variable time step increment is used: the time step employed during the first second is Dt = 1E-4, and it is increased to Dt = 0:1 for the remaining steps. It should be noted that Equation (59) is a first-order hyperbolic equation for fluid concentration. For a non-zero fluid concentration boundary condition applied at the inflow, the standard Galerkin finite element method results in oscillations around the steep front thereby causing numerical instability. We employ the VMS method as described in Section 4 to stabilize the formulation, and provide a comparison between the stabilized numerical result with the exact solution. Figures 2(a) and (b) show the performance of the new method for h-refinement wherein we have used linear Lagrange interpolation functions. These plots show the spatial profiles of the fluid concentration and velocity fields at 1000, 10,000, and 30,000 seconds. It can be seen that as the number of elements is increased, the computed solution converges to the exact solution, which is a numerical validation of the consistency of the formulation. Likewise, Figure 3 shows the convergence of the fluid density field for quadratic elements. Figure 4 shows the L2 norm of the error in the fluid density field for linear and quadratic VMS elements. Convergence rates of 1.54 for linear and 1.88 for quadratic VMS elements are obtained for the nonlinear first-order problem. Figures 5(a) and (b) show that numerical results compare well with the exact solution at 1000, 10,000, and 30,000 seconds, wherein the domain is discretized with 400 elements. In Figures 6(a) and (b), we show the spatial distributions of fluid density and fluid velocity for three different values of the drag coefficient for a domain of length 1 m. It can be seen that for a lower drag

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

12

Mathematics and Mechanics of Solids

Figure 2. Mesh refinement study at various time levels. (a) Fluid density: linear Lagrange h-refinement. (b) Fluid velocity: linear Lagrange h-refinement.

Figure 3. Mesh refinement study using quadratic Lagrange elements.

Figure 4. Convergence plot of the L2 norm of fluid density.

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

Hall et al.

13

Figure 5. Comparison between exact and finite element solution. (a) Fluid density along the domain. (b) Fluid velocity along the domain.

Figure 6. Variation of (a) fluid concentration and (b) fluid velocity for three different drag coefficients.

coefficient that corresponds to higher diffusivity, fluid propagates further down in the porous solid as compared to the cases of higher drag coefficients.

5.2 Oxidation of PMR-15 resin Thermo-oxidative aging of PMCs in high-temperature applications influences the life and performance of these materials. In this section, we present numerical results for the oxidation behavior of polyimide PMR-15 resin based on the oxidation reaction model developed in the work of Tandon et al. [17]. For the sake of completeness, we provide a brief description of the oxidation process in polymer. However, for a detailed description of the oxidation process and the reaction kinetics model, refer to [17,18]. The oxidation front in polymer materials advances through a combination of diffusion and reaction mechanism. The exposed surface reacts with the diffusing air that depletes the amount of polymer available in that region. Once this region is fully oxidized, it acts as a medium through which air/oxygen diffuses and an active oxidation zone is formed ahead of the fully oxidized zone. Thus, at any instant of time, the oxidation process in polymers comprises a fully oxidized zone, an active oxidation zone, and a neat resin zone, as shown in Figure 7.

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

14

Mathematics and Mechanics of Solids

Figure 7. Schematic representation of the thermo-oxidation process.

The oxidation reaction rate implemented in this work is given in [17] as f  fox )R(r f ) G_ 011 = ( 1  fox  R0 2j(1  0:5j) f.fox f R(r ) = 0 f  fox j=

 f br  f 1 + br

ð62Þ ð63Þ

ð64Þ

where G_ 011 is the rate of reaction, f is the oxidation state variable that indicates the availability of polymer for oxidation, fox specifies the fully oxidized state of the material (Zone I), R0 is the saturated rate  is the inverse of the saturation air/oxygen concentration. The evolution equation for of reaction, and b the oxidation state variable f is given as df = am f dt

ð65Þ

where a is the constant of proportionality. f varies in the active oxidation region while it assumes a value of fox in the fully oxidized region and a value of 1.0 in the unoxidized region. In the numerical test presented below, we consider a one-dimensional domain of length 1 mm. The left end of the domain is exposed to air, and the simulation is run under isothermal conditions at a uniform temperature of 288°C. Material parameters used in the simulation are given as follows: (i) the true air density at 288°C, rfT = 0:6273 kg m3 ; (ii) viscosity of air, AL = 29:5E-6 kg=m s; (iii) gas constant,  = 286:987 J=kg  K; (iv) body force, bf = 0; (v) molecular weight of air is MWair = 0:02897 kg=mol; R

 = 32:4412 m3 kg; (vii) oxidation state, fox = 0:187; (viii) reaction rate, R0 = 1:69E-2 kg m3 s; (ix) (vi) b

s 3 true solid density, of solid, fs = 0:1; (xi) diffusivity of the

rT = 1320 kg m ; (x) porosity

solid, 2 s D = 8:933E-13 m s; (xii) Young’s modulus, E = 2:6 GPa; (xiii) k2f = 0; (xiv) a = 0:35 m3 kg; (xv)  0G =  0:25E12; and (xvi) K  =  1:0E9. It is noted that the only new parameters that are not conA 1  0G and K  1 . The remaining parastrained by direct measurements are the last two parameters, that is, A meters are either specified in the original work [17], or are standard reported values (limited to the viscosity of air AL = 29:5E-6 kg=m s and Young’s modulus of PMR-15 Es = 2:6 GPa). The one-dimensional domain is discretized with a graded mesh of linear Lagrange elements. The subset of the domain, [0,0.0001]m, is discretized with 100 elements and the rest of the domain also with 100 elements . The fluid and the solid constituents coexist over this domain. A fluid concentration of 22.8863E-3 kg/m3 is specified as the boundary condition for the fluid and a load of 1 atm is applied on the solid at the left end of the domain. The problem is run with time steps of 1E-5 seconds for 1000 steps, followed by a time step of 1E-3 seconds for 10,000 steps, and with a time step of 0.1 seconds for a total time of 100 hours. The drag coefficient Av for the oxidation problem is defined in terms of diffusiv v f f  ity of the solid, as A = r (r + r )Ru rT mD, where m is the Lagrange multiplier. For the derivation of this expression, refer to the Appendix.

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

Hall et al.

15

Figure 8. Oxidation layer growth with time for various values of the oxidation state variable.

Figure 9. Oxidation layer growth with time for various values of reaction rate.

Remark: In our model, the fluid properties and its initial/boundary conditions are defined in mass concentration units. Since fluid properties in Tandon et al. [17] are provided in molar concentration units, they have been converted to appropriate units for the present system of equations using the standard conversion relations.

The active oxidation zone that lies between the fully oxidized zone and the unoxidized core has a continuous variation of f from fox to 1, respectively. Figure 8 shows the positions of the actively oxidizing domain (Zone II, Figure 7) versus time, for four different values of f associated with a reaction rate of 1.69E-3 kg/m3-s and a solid diffusivity of 8.93E-13 m2/s. The upper curve f=1 depicts the position of the Zone III–II boundary (oxidation front) between the oxidized and unoxidized regions, and the lower curve f=fox depicts the position of the Zone II–I boundary of the actively oxidizing zone and the fully oxidized zone. The remaining values of f correspond to positions within the actively oxidizing domain. The oxidation layer growth results shown in Figures 9–11 are plotted for f = 0:3. A parametric study was done for the oxidation layer growth with time and results are shown in Figures 9–11. Figure 9 shows the variation in oxidation layer growth for different reaction rate parameters for a duration of 100 hours. The solid line shows the results from the mixture theory, where it can be seen that the reaction rate of 2.41E-4 kg/m3-s produces an oxidation layer growth of 66.9 mm as

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

16

Mathematics and Mechanics of Solids

Figure 10. Oxidation layer growth with time for various values of oxidation state.

Figure 11. Oxidation layer growth with time for various values of drag coefficient.

compared to 74.7 mm for the reaction rate of 1.69E-3 kg/m3-s at the end of 100 hours. The mixture theory results follow a similar trend in comparison with the Tandon et al. [17] numerical results. Figure 10 shows the growth of oxidation layer for 0.1 and 0.187 oxidation state values. Since f 2 ½fox , 1, the local value of f indicates the amount of polymer that is available for oxidation. An oxidation state value of 0.1 indicates the spatial location where almost 90% of the polymer is available for oxidation, as compared to a value of 0.187 that indicates that only 81.3% of the polymer can be oxidized. For a constant oxidation rate, a lower value of fox indicates that there are more sites available for oxidation and therefore the rate of growth of the oxidation layer will be slower, as can be seen in Figure 10. Figure 11 shows the influence of the diffusivity of the solid on oxidation layer growth in PMR-15 resin. It can be observed that a diffusivity value of 1.667E-12 m2/s advances the oxidation layer at a higher rate in comparison to the lower diffusivity values of 1.3E-12 and 8.933E-13 m2/s. The oxidation layer depths of 74.7, 90.1, and 100.1 mm are observed for solid diffusivity values of 8.933E-13, 1.3E-12, and 1.667E-12 m2/s at the end of 100 hours, respectively. Tandon et al. [17] studied the oxidation layer growth via the diffusion-reaction equation assuming an ideal fluid permeating through a rigid solid. Accordingly, in their model the deformation of the solid and viscous effects in the fluid are neglected. In the present work where we employ the mixture theory, a Newtonian fluid and an elastic solid are considered. Since the unknown fields in the mixture model are

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

Hall et al.

17

Figure 12. Fluid and solid kinematic and force quantities along the domain at the end of 100 hours. (a) Solid density along the domain. (b) Fluid density along the domain. (c) Fluid stress along the domain. (d) Interactive force along the domain.

fluid density, fluid velocity, solid displacement, and solid density, the kinematic and the force measures can be readily obtained from the simulations. Figure 12 shows the variation of the fluid and the solid kinematic and force quantities for solid diffusivity values of 8.93E-13, 1.30E-12, and 1.67E-12 m2/s. The plots shown are obtained for a saturated oxidation state value of fox = 0.187 and a reaction rate of 1.69E-3 kg/m3-s. Figures 12(a) and (b) show the variation of solid density and fluid density along the domain at the end of 100 hours. Full oxidation of all presumed available sites results in a fixed solid density, as indicated in Figure 12(a). Since there are only two constituents in the present model, loss of mass from one is the gain in mass for the other. Consequently, the density of the solid increases as shown in Figure 12(a) wherein the apparent solid density has a higher value as compared to the neat resin region. This is rather contradictory to the experimental observations, as the density of the PMR-15 resin is expected to decrease with increased levels of oxidation. (It does, however, correspond to initial weight gains in certain oxidizing material systems before substantial mass loss to the environment occurs. The transfer of mass out of the material system is not explicitly addressed here.) If the two-constituent mixture model is extended to the three-constituent model where the third constituent is allowed to evolve and also leave the domain, it can account for the experimentally observed weight loss in solid due to the oxidation process. Figure 12(c) shows that the variation in fluid stress is dominated by the hydrostatic pressure. Figure 12(d)

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

18

Mathematics and Mechanics of Solids

shows the distribution of interactive force between the diffusing fluid and deforming solid. It can be seen that the interactive force becomes zero in the neat resin region where the fluid has not reached yet.

5.3 Slurry infiltration problem SI is an important step in the processing of ceramic matrix composites (CMCs). In the SI process, a viscous fluid that is laden with particles of various sizes, composition, and volume fraction is injected into a fiber preform, wherein fluid primarily serves as a medium that carries the suspended particles to the preform. This cycle is repeated several times until the density of the preform increases and its porosity reduces to some desired design value. Once the SI process is complete, a second process called melt infiltration is carried out with a viscous fluid that can chemically react with the preform as well as the deposited particles to make a composite with the desired strength and density distribution [26]. In this section, we consider SI and employ properties of a porous PMC as a surrogate model for CMC material. We assume that water-based slurry has permeated the porous elastic solid and we model the process of deposition of suspension onto the fiber preform. Young’s modulus of the porous PMC is obtained via the rule of mixtures as given below: EL = Ef Vf + Em Vm

ð66Þ

where Ef , Em are the fiber and epoxy Young’s moduli, respectively, and are assigned values of 380 and 3.45 GPa. Vf , Vm are the volume fractions of the fiber and the matrix in the porous composite. For a 50% porous PMC, we evaluate the properties based on 40% fiber and 10% matrix composition. The carbon fiber density and the matrix density are 1950 and 1200 kg/m3, respectively. The water-based slurry is assumed to contain SiO2 particles of dimension 2–15 mm with 50% volume fraction. The viscosity of the slurry can be computed from Einstein’s equation [27] as follows: msl = mw (1 + 2:5fSiO2 )

ð67Þ

where msl is the viscosity of slurry, mw is the viscosity of water, and fSiO2 is the volume fraction of SiO2 particles in the slurry. Assuming 50% volume fraction of SiO2 particles, the slurry viscosity turns out to be 1.793E-3 kg/m-s. Given that the density of the SiO2 particles is 2650 kg/m3 and the density of water is 1000 kg/m3, slurry density can be computed as rsl = 0:5 3 rw + 0:5 3 rSiO2

= 1825 kg m3

ð68Þ

where rsl , rw , and rSiO2 are the density of the slurry, water, and the SiO2 particles, respectively. In the present model, it is assumed that the particle-laden fluid is uniformly present in the domain and the dependence of the rate of deposition on the flow velocity is ignored. Accordingly, the mass deposition rate of particles from the slurry onto the porous composite, as given in [24], is modified for the present case as follows: m_ f = kr f w

ð69Þ

where k is the filtration constant and w is the apparent mass fraction of the particles in the slurry. The filtration rate of the solid medium is assumed to be 83.8341E-3 s21. The initial apparent mass fraction of particles in the slurry can be computed as w0 = fs Vslp

rp rsl

ð70Þ

where fs is the solid porosity (fs = 1  rs rsT , rsT is the true density of the solid) and rp , Vslp are the density and volume fraction of particles in the slurry, respectively. The drag coefficient Av for the SI (perL meation)

problem is defined in terms of permeability of the solid K and the viscosity of the fluid A as v L A = A K. (For the derivation of this expression, see the Appendix). The permeability of the solid is

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

Hall et al.

19

Figure 13. Mixture constituents kinematic and stress measures along the domain at 30, 60, and 90 seconds. (a) Fluid density along the domain. b) Solid density along the domain. (c) Hydrostatic stress in the fluid along the domain. d) Axial stress in the solid along the domain.

Figure 14. Reduction in solid porosity with time.

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

20

Mathematics and Mechanics of Solids

Figure 15. Reduction in solid porosity with time for 30%, 40%, and 50% SiO2 particles in the slurry.

Figure 16. Reduction in solid porosity with time for 40%, 50%, and 60% initial solid porosity.

taken to be 4.935E-17 m2. The chemical reaction and solid strain coupling parameters are assigned to be  0G =  0:25E3, K  =  1:0E0. A 1 In this problem, a one-dimensional domain of length 0.4 m is considered that contains both solid and fluid constituents that are uniformly present everywhere. We assume uniform material properties and temperature distribution. In addition, we assume that deposition of the suspended particles is occurring throughout the domain. The problem is run for 90 seconds with a time step of 5E-4 seconds. The solid displacement and fluid velocity are constrained at the left end of the domain. A load of 1E7 N is applied at the right end of the domain. Figures 13(a) and (b) show the reduction in apparent fluid density and increase in apparent solid density along the domain at the end of 30, 60, and 90 seconds. The initial apparent fluid density of 912.5 kg/m3 drops to 514.9 kg/m3 at the end of 30 seconds and further drops to 393.4 and 337.1 kg/m3 at the end of 60 and 90 seconds, respectively. This drop in fluid density is due to particle deposition on to the porous solid that results in an apparent solid density increase (see Figure 13(b)) from 900 kg/m3 to 1297.7, 1419.1, and 1475.5 kg/m3 at the end of 30, 60, and 90 seconds, respectively. In order to evaluate the evolution in the stress-carrying capacity of the solid, an external load is applied to the solid, which is held constant in time, that is, the solid is under constant compressive stress of 10 MPa throughout the process. Figure 13(c) shows the hydrostatic stress in the fluid, and Figure 13(d) shows the solid axial stress profiles along the domain. Since the rate of deposition of particles is constant along the

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

Hall et al.

21

domain, the fluid and solid stresses remain constant along the domain at a given (but otherwise arbitrary) time level. Figure 14 shows the decrease in solid porosity as a function of time. For a 50% initial solid porosity and with a 50% particle slurry, the maximum reduction in porosity is bounded by 0.25. As can be seen from Figure 14, the solid porosity asymptotes to 0.25 with time. Next, we present the results for the case where the porous solid is subjected to three infiltration cycles of 30 seconds each, for a total of 90 seconds. At the end of each cycle, the particle mass fraction w is reset to the initial particle mass fraction in the slurry w0 . Figure 15 shows the variation of the solid porosity with time for 50% porous solid and 30%, 40%, and 50% SiO2 particle volume fraction in the slurry. We see that as the particles get deposited, the porosity of the solid decreases. For all three different particle volume fractions in the slurry, this decrease in porosity is nonlinear, wherein the rate of reduction in porosity seems to be slowing down with time, which is indicated by the relatively flatter portion of the curve at the end of each cycle. From the perspective of the physics of the problem this means that while there is more relative reduction in porosity during early infiltration cycles, because of the closure of pores that happens due to the solid mass buildup, the relative reduction in porosity in subsequent cycles also slows down. Figure 16 shows a similar trend in reduction in porosity with time for three different initial solid porosities that are infiltrated with 50% particle slurry.

6. Conclusions We have presented a VMS-based finite element method [13,15,16] for the fluid–solid mixture theory model of Hall and Rajagopal [12] that is based on the constituent equations of motion and mass balance. The model addresses the energy and entropy production equations through an equation for the Lagrange multiplier that results from consideration of the full set of balance equations as a constraint during the process of maximization of entropy production. The present system of equations is applied to isothermal processes in the one-dimensional context. Employing VMS ideas, a multiscale decomposition of the fluid density field into coarse and fine scales and a priori unique decomposition of the admissible spaces of functions leads to two coupled nonlinear problems termed the coarse-scale and the finescale sub-problems. The fine-scale solution is extracted from the nonlinear fine-scale sub-problem, which is then variationally projected onto the coarse-scale space, leading to a formulation that is expressed entirely in terms of the coarse scales. Although the final formulation does not depend explicitly on the fine-scale density field for the fluid, the effects of fine scales are consistently represented via the additional residual-based terms, and they add to the stability of the numerical method. The resulting stabilized method for the mixture model is applied to hyperbolic propagation while recovering Fickian diffusion, anisotropic oxidation in composite materials recovering the data of Tandon et al. [17], and mass deposition. Results of the oxidation modeling of Tandon et al. [17] are recovered by employing the reaction kinetics model and properties assumed therein; the only additional assumed properties are two constants describing coupled chemomechanical and purely chemical dissipation. In all of these cases the mixture provides rich detail concerning the kinematic and kinetic behaviors of the constituents, in contrast to standard effective media approaches. The proposed solution scheme based on a single Helmholtz energy reveals the importance of an effective material property related to the temperature rate; further investigation in the three-dimensional context is needed to determine applicability to general anisotropic and anisothermal problems. Funding Partial support for this work was provided by AFRL under Contract No. FA8650-13-C-5214. This support is gratefully acknowledged.

Conflict of interest statement None declared.

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

22

Mathematics and Mechanics of Solids

References [1] Atkin, R, and Craine, R. Continuum theories of mixtures: basic theory and historical development. Q J Mech Appl Math 1976; 29: 209–244. [2] Green, AE, and Naghdi, PM. On basic equations for mixtures. Q J Mech Appl Math 1969; 22: 427–438. [3] Green, AE, and Naghdi, PM. A theory of mixtures. Arch Rational Mech Anal 1967; 24: 243–263. [4] Rajagopal, K, and Tao, L. Mechanics of mixtures. Singapore: World Scientific Publishing Company, Incorporated, 1995. [5] Rajagopal, K, and Srinivasa, A. On the thermomechanics of materials that have multiple natural configurations Part I: Viscoelasticity and classical plasticity. Zeitschrift fu¨r angewandte Mathematik und Physik ZAMP 2004; 55: 861–893. [6] Rajagopal, K. Diffusion through polymeric solids undergoing large deformations. Mater Sci Technol 2003; 19: 1175–1180. [7] Karra, S. Diffusion of a fluid through a viscoelastic solid. arXiv preprint arXiv:1010.3488, 2010. [8] Karra, S, and Rajagopal, K. A model for the thermo-oxidative degradation of polyimides. Mech Time Depend Mater 2012; 16: 329–342. [9] Humphrey, J, and Rajagopal, K. A constrained mixture model for growth and remodeling of soft tissues. Math Model Meth Appl Sci 2002; 12: 407–430. [10] Rao, I, and Rajagopal, K. A thermodynamic framework for the study of crystallization in polymers. Zeitschrift fu¨r angewandte Mathematik und Physik ZAMP 2002; 53: 365–406. [11] Malek, J, and Rajagopal, K. A thermodynamic framework for a mixture of two liquids. Nonlinear Anal Real World Appl 2008; 9: 1649–1660. [12] Hall, R, and Rajagopal, K. Diffusion of a fluid through an anisotropically chemically reacting thermoelastic body within the context of mixture theory. Math Mech Solid 2012; 17: 131–164. [13] Hughes, TJR. Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods. Comput Methods Appl Mech Eng 1995; 127: 387–401. [14] Masud, A, and Hughes, TJR, A space-time Galerkin/least-squares finite element formulation of the Navier-Stokes equations for moving domain problems. Comput Methods Appl Mech Eng 1997; 146: 91–126. [15] Masud, A, and Khurram, R. A multiscale/stabilized finite element method for the advection–diffusion equation. Comput Methods Appl Mech Eng 2004; 193: 1997–2018. [16] Masud, A, and Kwack, J. A stabilized mixed finite element method for the first-order form of advection–diffusion equation. Int J Numer Method Fluids 2008; 57: 1321–1348. [17] Tandon, G, Pochiraju, K, and Schoeppner, G. Modeling of oxidative development in PMR-15 resin. Polym Degrad Stab 2006; 91: 1861–1869. [18] Schoeppner, G, Tandon, G, and Pochiraju, K. Predicting thermooxidative degradation and performance of high-temperature polymer matrix composites. Multiscale modeling and simulation of composite materials and structures. Springer US: Springer, 2008, pp.359–462. [19] Varghese, J, Owens, BC, and Whitcomb, JD. Simulation of oxidation in textile composites. J Composite Mater 2011; 45: 1771–1782. [20] Varghese, J, and Whitcomb, J. Micromechanics of oxidation in composites with impermeable fibers. J Composite Mater 2009; 43: 2011–2043. [21] Hughes, TJ, Liu, WK, and Zimmermann, TK. Lagrangian-Eulerian finite element formulation for incompressible viscous flows. Comput Methods Appl Mech Eng 1981; 29: 329–349. [22] Khurram, RA, and Masud, A. A multiscale/stabilized formulation of the incompressible Navier–Stokes equations for moving boundary flows and fluid–structure interaction. Comput Mech 2006; 38: 403–416. [23] Calderer, R, and Masud, A. A multiscale stabilized ALE formulation for incompressible flows with moving boundaries. Comput Mech 2010; 46: 185–197. [24] Civan, F. Porous media transport phenomena. New Jersey: John Wiley & Sons, 2011. [25] Bowen, R. Porous elasticity: Lectures on the elasticity of porous materials as an application of the theory of mixtures, http://hdl.handle.net/1969.1/91297 (2014, accessed 22 January 2014). [26] Bansal, NP, and Boccaccini, AR. Ceramics and composites processing methods. New Jersey: Wiley, 2012. [27] Hosokawa, M, Nogi, K, Naito, M, et al. Nanoparticle technology handbook. Oxford: Elsevier, 2007.

Appendix – Relation between solid diffusivity and drag coefficient In Section 5, we have presented a reduced-order mixture problem, oxidation of PMR-15 resin, and SI problem. A literature review reveals that the reduced-order mixture problem and the oxidation problem is in general modeled via a diffusion-reaction equation, while the SI problem is typically modeled via a Darcy equation. In the context of the mixture theory model, the fluid–solid interaction is accounted for via the interactive force field, which requires the specification of drag coefficient Av . The relation between the drag coefficient Av and the diffusivity of the solid D can be obtained by comparing the

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

Hall et al.

23

mixture theory equations and Fick’s diffusion-reaction equation. Similarly, the relation between the drag coefficient Av and the permeability of the solid K can be obtained by comparing the mixture theory equations and the Darcy equations for the SI problem. A1. Fick’s diffusion-reaction equation

Fick’s diffusion-reaction equation can be written as ∂r f ∂2 r f = D 2 + mf ∂t ∂x

ð71Þ

where D is the solid diffusivity. A2. Darcy equation

The fluid balance of the mass equation and Darcy’s law are given as follows: u1f = 

K ∂p AL ∂x

∂r f ∂(r f v1f ) = mf + ∂x ∂t

ð72Þ ð73Þ

where u1f is the filtration velocity, K is the permeability of the solid, and AL is the viscosity of the fluid.  Equations (72) and (73) can Assuming that the pressure of the fluid follows the ideal gas law, p = r f Ru, be combined as follows:  f  ∂ ∂r f K Ru f ∂r ð74Þ r  L s = mf A f ∂x ∂t ∂x where fs is the solid porosity. Equation (74) can be written in an expanded form as  f 2  f ∂2 r f K Ru  ∂r f K Rur ∂r  L s = mf  L s 2 A f ∂x ∂t A f ∂x

ð75Þ

A3. Mixture theory

The fluid balance of mass and linear momentum are given as ∂r f ∂(r f v1f ) = mf + ∂x ∂t

ð76Þ

f ∂T11 + I1f = 0 ∂x

ð77Þ

where the fluid body force and inertial effects are neglected.

A3.1 Reduced-order mixture problem Consider the constitutive relations for the fluid stress and interactive force as given in Equations (57) and (58). Substituting these constitutive relations into Equation (77), the fluid velocity can be written as v1f = 

 ∂rf Ru rf Av ∂x

ð78Þ

Equations (76) and (78) can be combined to give  ∂2 rf ∂r f Ru  v 2 = mf A ∂x ∂t

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014

ð79Þ

24

Mathematics and Mechanics of Solids

Comparing Equations (71) and (79), the drag coefficient can be written in terms of solid diffusivity as D=

 Ru Av

ð80Þ

A3.2 Oxidation and slurry infiltration problem Consider a simplified form of the constitutive relations for the fluid Equations (31) and (32), as given below: f T11 ’  rf r

∂c rf r  ’  Ru ∂r f rT

I1f ’  mAv v1f

ð81Þ ð82Þ

Substituting the above Equations (81) and (82) into Equation (77), the fluid velocity can be written as   ∂ rf r Ru f ð83Þ v1 ’  mAv ∂x rT Fluid velocity given in the above expression is substituted in the fluid balance of mass, Equation (76), and is written as follows:   f  ∂ Ru ∂r f r r f ∂ r  ð84Þ = mf mAv ∂x ∂x rT ∂t 2 f  f Rur ∂r f f ∂ r (r + r ) +    = mf  ∂t mAv rT ∂x2

ð85Þ

Comparing Equation (85) and Equation (71), we can obtain the following relation for solid diffusivity and drag coefficient for the oxidation problem as D=

 f Rur (r + r f ) mAv rT

ð86Þ

Comparing Equation (85) and Equation (75), we can obtain the relation between the drag coefficient and the permeability of the solid for the SI problem as Av =

AL fs (r + r f ) KmrT

ð87Þ

In the Section 5.3, we have presented numerical results for a simplified form of the above relation, AL Av ’ . K

Downloaded from mms.sagepub.com at UNIV OF ILLINOIS URBANA on August 22, 2014