Surgical implants, known as Tissue Engineered Nerve Guides ( TENG s) are often inserted to provide alignment and mechanical support to broken or damaged nerves. These implants are designed to be poroelastic and biodegradable, and the ideal rate of degradation would be that equal to the rate of regrowth of the recovering nerve. Inspired by the design of these TENG s, we develop a mixture theory based mathematical model to simulate the degradation of a poroelastic solid immersed in a uid bath. The temporal evolution of the solid's mechanical and transport properties is also modeled. The model comprises of the degrading solid, the degradation reaction products, and the uid in which the solid is immersed. The resultant governing equations are formulated in an Arbitrary Lagrangian{Eulerian ( ALE) framework. The weak formulation of the partial di erential equations ( PDE s) so derived is numerically implemented using a nite element method ( FEM ). The numerical model is studied for stability and convergence rates using the Method of Manufactured Solutions. Keywords: Fluid-Structure Interactions, Finite element method, Mixture theory, Poroelasticity 1 Introduction Peripheral nerves carry sensory-motor impulses from the brain to di erent parts of the body and vice versa. Injuries to peripheral nerves can have dire consequences, such as loss of sensitivity, motion, or both, in the injured organ. The inherent mechanism of the body to mend the damaged or broken nerves can have retargeting errors when rejoining broken nerve ends with large gaps. External support for realignment and mechanical strength of such nerve injuries involving large gaps is often provided in the form of surgical nerve autografts [ 1,2] (nerve segments obtained from a healthy donor site in the patient). Unfortunately, this method has several shortcomings, such as scarcity of viable donor cites, possibility of neuroma at the donor location and mismatch between the size of the graft and injured nerve. Tissue Engineered Nerve Guides ( TENG s) are ∗Corresponding author 1 an attractive alternative to nerve autographs. TENG s are sleeve-like porous sca olds, preferably biodegradable, that are sutured to the severed nerve ends [ 1,3].TENG s are expected to degrade and be bioresorped as the nerve regrows. Ideally, the rate of regeneration of the nerve matches theTENG 's degradation rate. Moreover, TENG s are expected to provide mechanical support and alignment to the broken nerve ends and are usually designed to be porous in order to enable e ective supply of drugs and nutrients to the recovering nerve [ 1]. Thus, TENG design requires careful study of their degradation chemistry and the corresponding evolution of their transport and mechanical properties. In this study, we have derived a formulation capable of simulating the complex chemistry, mechanics, and transport aspects of a typical TENG . An example of a typical TENG material is cross-linked urethane-doped polyester elastomers ( CUPE ) [4,5], and our eventual goal is to develop a model which can accurately estimate the degradation and mechanical response of this material when inserted at a target cite. Studies of degradation of biopolymers (e.g. [ 6,7]) show that the rate of hydrolytic degradation of elastomers is governed by the concentration of the polymer as well as the degradation reaction products. When a poroelastic polymer sca old is inserted in a dynamic environment such as in a limb, it will likely be subject to stresses along with di usive momentum transfer between the bodily uids and the solid. This necessitates the study of the mass and momentum balance equations of the polymer, the degradation products and the uid in which the sca old is inserted. Mixture theory provides an excellent framework for the study of all components in the mixture. Thus, we resorted to this continuum framework for the derivation of our governing equations. Reactive mixture theoretic formulations for biomechanical applications have been limited mostly to aneurysms (c.f. e.g [ 8]) and growth and remodeling of cartilage, bone and soft tissues (cf., e.g., [9{15]). Most of these works are based on volumetric or surface growth, through decomposition, without directly involving chemical kinetics. Moreover, the growth or remodeling is assumed to occur mostly through mechanical loading. Considerable focus of these studies is mainly on the evaluation of residual forces, tissue swelling or other mechanical phenomena. Constitutive models for the growth in biological tissues by applying the Clausius-Duhem inequality have been derived in [15], in which the authors derive coupled mass transport-mechanics based formulation for growth of tissues, by considering the viscous uid and solid to be a single vicoselastic entity. Balance equations are then obtained for this viscoelastic tissue. They obtain the formulation in terms of entities such as the mixing entropy, which may not be a readily available input data. Moreover, as mentioned earlier, a component-wise modeling of our problem is essential, as the chemical kinematics of the degradation of the poro-elastic solid may depend on the local concentrations of various components in the mixture. Thus, from a modeling point of view, degradation cannot be viewed simply as a reversal of tissue growth. Chemo-mechanical coupling has been modeled for non-biomechanical applications as well. For example, in [ 16] and [ 17] a model for the oxidative degradation of polyimides at high temperatures has been derived. The evolution of properties in a chemically reacting mixture is determined mainly by two processes - di usion and chemistry. Oxidative degradation of polyimides and growth of biological tissues are examples of reaction-dominated and di usion-dominated thermochemical processes respectively. An excellent derivation of these two asymptotic extremes in a di usive- reactive system consisting of a visco-elastic uid di using through a thermoelastic body have been derived in [ 18]. Di usion-dominated processes lead to occurrence of reactions all over the body, whereas reaction-dominated processes often lead to advancing reaction fronts in the body, thus the modeling approach for these two cases is di erent. A more general modeling approach is essential when the di usion-reaction rates are comparable, or when a more generalized formulation is sought. Hence, in our formulation, the chemistry as well as mechanical and transport properties of di erent components in the mixture are derived as a function 2 of the spatial location of material points in the mixture domain. Any heterogeneity in the mixture is naturally taken care of in this approach, which also allows for more exibility in determining the local behavior of the mixture, based on the responses of its components. We choose a classical continuum mechanics approach for mixtures (c.f. e.g. [ 19,20]) for the derivation of constitutive models for ows through poro-elastic degrading solids. Balance equations are derived for each component in the mixture and constitutive forms are then derived for stresses in the solid as well as forces of interaction between several components of the mixture. The mixture theoretic equations governing the motion of uid ow through poro-elastic solids have been derived in [ 21] and implemented using anFEM framework in [ 22]. In this work, we extend this theory to include bio-degradation of the poro-elastic polymer. This constitutes a reactive uid-structure interaction ( FSI) problem. We extend our previous non-reactive FSIformulation by including the relevant terms in the mass and momentum balance equations of all the mixture constituents and re-deriving the formulation. In this preliminary e ort, we assume that a single product of the degradation reaction is produced, which remains within the system. Thus, we model our thermochemical system as a three component mixture, comprising of the poro-elastic solid, the degradation product and the base uid in which the solid is inserted, assuming all the components to be incompressible. The constitutive equations are derived by invoking entropy maximization for the entire mixture (second law of thermodynamics for mixtures), which is a standard technique in mixture theory, [19, 23]. Notations used in the paper, along with kinematics and general mixture theoretic balance laws are given in Section 2. Constitutive equations for the speci c case of a degrading poroelastic solid immersed in a uid are derived in Section 3. Here, we also derive the Eulerian strong form of the governing equations, as well as their nal and simpli ed quasi-static ALE strong form. Corresponding weak formulation is presented in Section 4, along with the introduction of functional spaces for the trial and weighting functional spaces. The results pertaining to the numerical implementation of the FEM scheme are discussed in Section 5, along with a discussion on the convergence and stability of the numerical scheme. We conclude the paper with a summary and discussion section. 2 Kinematics and balance laws for mixtures 2.1 Notation and kinematics For a mixture of Nconstituents, we posit that each of these has its own reference con guration. Figure 1 shows a 3-component mixture, where aandXa(a= 1;:::; 3) denote the reference con guration and the referential particle position of component a, respectively. For our analysis, we assume that the reference con gurations for all components coincide with the initial con guration, and hence, Xaalso denotes the initial location of the particle. tdenotes the current con guration. aand tare subsets of End, thend-dimensional Euclidean point space, d= 2;3. The translation vector space of Endis denoted by Tnd. The particle at Xais mapped to a spatial point xin the current con guration via the di eomorphism a: a[0;T]! tgiven by, x= a(Xa;t): (1) Xais also used to denote the inverse map of the above mapping, that is, Xa= 1 a(x;t); (2) where 1 a: t! s. The properties of the particle are similarly mapped and hence evolve with time. As an example, if a constituent is generated in the mixture as a result of a chemical reaction, then the mass density 3 X1∂Ω1n1 ∂Ω2Ω1 χ1(X1,t)X2n2Ω2 χ2(X2,t) OxΩt∂Ωt ntX3 χ(,t)3X3∂Ω3 Ω3n3Figure 1: Con gurations: 1; 2; 3are the reference con gurations of species 1, 2 and 3 respectively, and X1;X2;X3are labels given to particles of the respective constituents. tdenotes the current con guration, and xis the common image of the di eomorphisms 1, 2, 3.na's denote the unit normals to the boundaries @ a's. Section 3 onwards, the subscript 1 will be replaced bysand will denote the solid component in a mixture consisting of N1 uids and a degrading poroelastic solid. Numerical subscripts will denote various uids in the mixture. of a `particle' of this constituent in the reference frame may increase with time. Also, using Eqs. (1) and(2), for any property  (Xa;t) in the reference frame (or, following the particle) currently at x, given by  (x;t) in the current frame, we can write, =  (Xa;t) =  (x;t) =   1 a(x;t) =   a(Xa;t): (3) Henceforth, any property will simply be denoted by , and its domain as a function will be decided based on the context. As is common in mixture theory, all components are assumed to co-exist at any given point xin the current con guration. Each constituent is characterized by a mass density distribution a=a(x;t), which represents the mass of component aper unit mixture volume. Then, the mixture mass density is de ned as =(x;t) =NX a=1a(x;t): (4) We assume that each constituent is incompressible. In this case, it is meaningful to talk about the intrinsic or true mass density of the constituent, which we denote by a, and which is the mass of constituent awhen said constituent is in pure form. The incompressibility assumption is embodied by the statement that a(a= 1;:::;N ) is a (known) constant. Relating ato ais done 4 by introducing the notion of volume fraction of constituent a, which we denoted by aand de ne through the following relation: a= a a: (5) We assume that the constituents of the mixture occupy the volume lled by the mixture. This assumption is expressed by the following relation, referred to as the saturation condition: NX a=1 a= 1: (6) Spatial gradients with respect to xandXa, respectively, are denoted as, rx=@ @xandrXa=@ @Xa: (7) Fora= 1;:::;N , We de ne the displacement, ua, deformation gradient, Faand its Jacobian Ja as ua(Xa;t) = a(Xa;t)Xaorua(x;t) =x 1 a(x;t); (8) Fa(Xa;t) =rXa a(Xa;t); (9) Ja(Xa;t) = det Fa(Xa;t) : (10) Also, the right Cauchy-Green tensor is given by Ca=FT aFa. The Eulerian description of the material velocity of constituent ais given by va(xa;t) =@ a(Xa;t) @t Xa= 1 a(x;t)=@ua @t Xa= 1 a(x;t): (11) Let a(x;t) be the Eulerian description of a quantity pertaining substance a. Then the material time derivative of ais denoted by  aand is de ned by,  a=@ a( a(Xa;t);t) @t=@ a(x;t) @t+rx a(x;t)va; (12) where, for generic vectors uandv,uvdenotes the inner product of uandv. The barycentric velocity of the mixture is de ned by v=1 NX a=1ava: (13) The time derivative of a quantity (x;t)following the barycentric velocity is denoted by _ and is given by, _ =@ (x;t) @t+rx v: (14) We will denote by ~vathe di usion velocity of substance a, which is given by: ~va=vav; (15) Another important concept for our analysis is that of ltration velocity of any uid b, which we de ne as, v lb= b(vbvs); (16) which is the relative velocity of the uid with respect to the solid, scaled by the uid volume fraction. Details of ltration velocity can be found in [ 21], and further discussion from the numerical FSIpoint of view can be found in [ 22]. In the text, we derive constitutive equations for a mixture consisting of one incompressible poroelastic solid and N1 incompressible uids, and later specialize the formulation for 2 uids. 5 2.2 Balance of mass The balance of mass for species ain an Eulerian framework is given by, @a @t x+rx(ava) = ^ca; (17) where ^cais the rate of production of species aper unit volume of the mixture. Dividing by a, the true mass density of the species a, and recalling that ahas been assumed to be a constant, we get @ a @t x+rx( ava) =^ca a= ^ a; (18) where the term ^ ais the rate of production of speci c volume (volume of a species per unit volume of the mixture) of constituent a. Rewriting the equation in terms of the ltration velocity, we get, @ a @t xr x( avs) +rxv la^ a= 0: (19) Pulling back Eq. (18) to the reference con guration of the solid, s, we have, @ a @t Xsr x avs+rx( ava)^ a= 0: (20) This remapping operation can be viewed as adopting an ALE framework (cf. [ 24]) in which the ALE map is the motion of the solid phase. Rewriting the ALE equation in terms of the ltration velocity, we obtain @ a @t Xs+ aFT s:rXsvs+FT s:rXsv la^ a= 0; (21) where, for any two (second order) tensors AandB,A:Bdenotes the (tensor) inner-product of A andB. Since v ls=0, in the ALE framework, the mass balance equation for the solid is @ s @t Xs+ sFT s:rXsvs= ^ s: (22) The mass balance equation for the mixture is postulated to be the result of adding all of the constituent-wise equations with the stipulation that: NX a=1^ca= 0: (23) Equation (23) is a constraint relation stating that no net mass is produced or consumed in the mixture. This said, it is important to observe that, in general,NP a=1^ a6= 0. Summing each contribution provided by Eq. (17), and separating the contribution of the solid phase from the other contributions, we get the mass balance equation for the mixture in the Eulerian framework, given by, rx vs+NX a=2v la! =NX a=1^ a: (24) Similarly, summing Eqs. (21) fora= 2;3;:::;N along with Eq. (22), enforcing the constraint in Eq. (6), the mass balance of the mixture in the ALE framework can be given the following form FT s:rXs NX a=2v la+vs! =NX a=1^ a: (25) 6 2.3 Balance of momentum For any component ain the mixture, the Eulerian form of momentum balance equation is given by a@va @t x+varxva |{z} ava=rxTa+aba+^pa; (26) where, for generic vector elds uandv, the notation urxvdenotes the action of the gradient of v ontou,∗and where Taandbaare the partial Cauchy stress tensor and external body force density, respectively. The term ^pais the momentum supply for constituent a. Sometimes, this quantity is also referred to as the local or internal force acting on constituent adue to its interaction with the other components in the mixture [ 19]. The constitutive relations for ^paandTadepend upon the type of mixture involved. Following [19], we de ne the stress tensor for the mixture as T=NX a=1(Taa~va ~va): (27) In addition, the momentum balance postulate for the mixture is taken to be @v @t x+vrxv =rxT+b; (28) where b=1 NP a=1aba. Equations (28) and(27), along with the balance of momentum equations for each constituent, yields the following constraint relation: NX a=1(^pa+ ^ca~va) =0,NX a=1(^pa+ ^cava) =0: (29) Another equivalent statement of the balance of momentum for the mixture is obtained by summing Eq.(26) over all components in the mixture, and applying the constraints given by Eqs. (23) and (29). So doing yields NX a=1a@va @t x+rxvava =rxTI+bNX a=1(^cava); (30) where TIis referred to as the internal Cauchy stress tensor of the mixture and is given by TI=NX a=1Ta: (31) 2.4 Energy balance and second axiom of thermodynamics Following [19], we choose the energy balance equation for the mixture as _= tr( TL)r xq+r+NX a=1a~vaba; (32) ∗This notation is typical in the uid mechanics literature (cf. [ 25]). In the continuum mechanics literature, the operation in question is typically denoted by (r xv)u(cf. [19, 23]). 7 along with the constraint, 0 =NX a=1[^a+~va^pa+ ^ca(a+1 2~v2 a)]; (33) whereaand^arepresent the internal energy and energy supply per unit volume of the component a, respectively. The quantities =NX a=1(a+1 2a~v2 a);q;andr denote the internal energy density, heat ux vector and external heat supply density for the mixture, respectively. The energy equation contributes to the governing equations of the system only when we are considering non-isothermal systems, in which temperature of the system varies. With this in mind, the energy equation also contributes to a framework for the discussion of the constitutive equations based on the entropy production inequality. We follow the treatment of this inequality in [19]. As such, we assume that all constituents of the mixture are at the same temperature . Furthermore, we denote by athe Helmholtz free energy of species aper unit mass of said species. Then, we also de ne the vector h=qNX a=1a(Ka+1 2~v2 aI)~va (34) to represent the mixture energy ux, where Kais the chemical potential tensor de ned by Ka= aITa=a: (35) In addition to handKa, we also introduce aas the speci c entropy of constituent a, as well as the following quantities: =1 NX a=1aaand a=a a=a(aa): (36) The quantities and adenote the speci c entropy of the mixture and the Helmholtz free energy of speciesaper unit volume of the mixture. With the above de nitions, we then have the general form of the second axiom of thermodynamics as [19, 21]: NX a=1 a_NX a=1tr (aKaLa)hrx NX a=1~va^pa+1 2^ca~v2 a 0; (37) where Lais the spatial gradient of va, i.e., La=rxva: (38) Using the de nition of ~va, the last term of inequality (37) can be adjusted so that the inequality can be written as, NX a=1 a_NX a=1tr (aKaLa)hrx NX a=1(^pa+1 2^ca(vavs))(vavs)0: (39) 8 Mass balance and saturation constraint If we add all the component-wise balance of mass equations (given by Eq. (17)), and use the saturation constraint given by Eq. (6), we get, NX a=1[varx a+ arxva^ a] = 0: (40) In terms of La, Eq. (40) can be re-written as, NX a=1[(vavs)rx a+ aI:La^ a] = 0: (41) 3 Constitutive relations We derive general constitutive relations for a mixture comprising of an elastic solid and N1 uids, and all the constituents are assumed to be incompressible. We assign mass production rates ( ^ca) to all species. The value of ^cais 0 in the case of species which are not chemcially reactive. We consider a mixture which allows for combined e ects of elasticity of the solid, heat conduction, di usion and buoyancy (or hydrodynamic di usive forces resulting from density gradients). The variables for which we would like to obtain constitutive forms are a; ; ^pa;Ta;hand^ a. We will rst derive the form of Ka, and then, with the help of Eq. (35), we will derive the form for Ta. We de ne a thermochemical process for a single temperature as a set of 7 N+ 4 functions whose values are given by x= a(Xa;t); ba=ba(Xa;t); ^pa=^pa(Xa;t); a= a(Xa;t); Ta=Ta(Xa;t);^ca= ^ca(Xa;t); a= a(Xa;t); q=q(x;t); r=r(x;t);  =(x;t);  =(x;t):(42) In writing Eqs. (42) we have made use of the relations in Eq. (3). We require that the constitutive equations be consistent with the second law of thermodynamics, along with relevant balance laws and the saturation constraint in Eq. (41). Enforcing this constraint along with the entropy inequality, inequality (39), we obtain a variational inequality written in terms of a Lagrange multiplier, , NX a=1 a_NX a=1tr ((aKa aI)La)hrx  NX a=1 (^parx a) +1 2^ca(vavs) (vavs) +^ a 0;(43) for all admissible thermochemical processes. Due to the indeterminacy introduced by the saturation constraint, constitutive equations need to be de ned for ( aKa aI) and ( ^par a) instead of Kaand^parespectively. Thus, for the assumptions for the mixture as stated in the beginning of this section, we assume general constitutive forms of several state variables as follows, ( a;;(^parx a);(aKa aI);h;^ a) =f(;rx;Fs;rxFs; b;rx b;va): (44) Here, b = 2, 3,. . . , N, theN1 uid components in the system. Equation (44) states that we view the quantities on the left of the equal sign to be functions of the kinematic variables listed as 9 arguments of the function f. In reality, it is well known that, as stated, the dependence implied by the left-hand side of Eq. (44) does not guarantee compatibility with the axiom of material frame-indi erence (cf. [19, 20, 26, 27]). We therefore specialize Eq. (44) as follows: f= ( a;;Fs^p a;FsK aFT s;Fsh;^ a) =f(;g;Cs;Gs; b;db;wb); (45) where, g=rXs=FT srx,Gs=rXsCs=FT srxCs,db=rXs b=FT srx bandwb=FT s(vbvs). The constitutive equations have exactly the same form as in the incompressible and non-reactive case in [ 21], except that we assume the dependence of aon the exact same parameters as the other state variables. In [ 14] (a mixture theoretic model of soft tissue growth and remodeling), the authors present arguments for a relatively simple constitutive form for a, given by a= a(Fs;a), a= 1;2;:::;N . We rather derive the most general constitutive forms using the assumptions in Eqs. (45), and introduce specializations later. Also, we de ne the internal free energy of the mixture as, I= I(;g;Cs;Gs; b;db;wa) =NX a=1 a: (46) Then, following the constitutive assumptions in Eqs. (44), we have,  a=@ a @@ @t+rxva +@ a @Cs:@Cs @t+rxCsva +@ a @g@g @t+rxgva +C@ a @Gs @Gs @t+rxGsva +NX b=2@ a @db@db @t+rxdbva +@ a @ b@ b @t+rx bva +@ a @wb@wb @t+rxwbva (47) Based on the assumptions in Eqs. (45), inequality (43) becomes, @ I @g_gC@ I @Gs _GsNX b=2@ I @wb_wb@ I @+ _+NX b=2 b@ I @dbrXstr(Lb) NX a=1 @ a @grxg+CNX a=1@ a @Gs r xGs! ~vatr" Fs K s+ 2@ I @Cs NX a=1~va @ a @wb! FT sLs# NX b=2tr" FsK bFT s b@ I @ bI@ I @dbdb I@ I @db FT sdb FT s+NX a=1~va @ a @wbFT s! Lb# mg NX b=2 ^ b +@ I @ b+ b 2(vbvs)2 +@ I @dbrXs^ b ^ s NX a=1 ^fa @ I @daNX b=2@ a @db! rxda! (vavs)0; (48) whereCdenotes the contraction operator as de ned in, for example, [ 27] and we have de ned m 10 and^faas, m=h+NX a=1@ a @~va; ^fa=Fs^p a@ I @ arx a+@ a @CsrxCs+NX b=2@ a @ brx b;(49) and, the following relations hold (as is evident from Eqs. (49) and (45)), m=Fsm(;g;Cs;Gs; b;db;wa); ^fa=Fs^f a(;g;Cs;Gs; b;db;wa):(50) The inequality (48) needs to hold for all admissible thermochemical processes, for all values of _g;dGs dt;dvb dt;La;rXstr(La);_andrxda, holding values of all other independent variables xed. Hence,8b= 2;:::;N1, we have @ I @g=0;@ I @Gs=0;@ I @vb=0;@ I @db b=0;@ I @=;@ I @daNX b=2@ a @db=0; (51a) NX a=1 @ a @grxg+CNX a=1@ a @Gs r xGs! ~va= 0;8r xg;8r xGs; (51b) K s=2@ I @Cs+F1 sNX a=1~va @ a @wb; (51c) FsK bFT s= b@ I @ bI+@ I @dbdb I+@ I @db FT sdbNX a=1~va @ a @wbFT s: (51d) Substituting these results in inequality (48), we have mrx NX a=1^fa(vavs)NX b=2 ^ b +@ I @ b+ b 2(vbvs)2 +@ I @dbrXs^ b ^ s0:(52) The rst four of Eqs. (51a) also show that I= I(;Cs; b;db); (53) so the dependence of Ion other independent variables is eliminated. Another important observation, stemming from inequality (52), is that the presence of volume fraction (or density) gradients can cause production of entropy, unlike the non-reactive case presented in [ 21]. This phenomenon is also observed in the case of reactive mixture of elastic materials, [26]. From Eqs. (51c) ,(51d) , and (35), it is easy to see that the partial stress tensors for the solid and uids, TsandTa,a= 2;:::;N , and the total internal stress tensor, TIare given by, Ts= ( s s)I+ 2Fs@ I @CsFT s; (54a) Tb= b b +@ I @ b @ I @dbdb I@ I @db FT sdbFT s+NX a=1~va @ a @wbFT s; (54b) 11 TI= ( I)I+ 2Fs@ I @CsFT sNX b=2 b@ I @ b@ I @dbdb I+@ I @db FT sda FT s ; (54c) where we have used the result:NP b=2NP a=1~va @ a @wbFT s=0, the proof of which can be found in, for example, [27]. 3.1 Simpli ed modelling assumptions The constitutive relations presented so far are the most general form of constitutive equations that can be derived under the assumptions in Eqs. (45) for ow through porous, chemically degrading media. In this study, we focus on a simpli ed form of constitutive equations, and the nite element based numerical implementation of the resultant model. As a rst step of simpli cation, we ignore the dependence of the partial free energy on gradients and higher order derivatives of the primary variables, and on the component velocities. That is, we adopt the simpli cation proposed by [ 21], so that we assume a= a(;Cs; b) (b= 2;3;:::;N ). Under these assumptions we see that Eqs. (54b) and (54c) simplify to, Tb= b b +@ I @ b I; (55) TI= ( I)I+ 2Fs@ I @CsFT sNX b=2 b@ I @ b; (56) while the constitutive form for the partial stress tensor for the solid is still given by Eq. (54a) . Thus, the constitutive forms of stress tensors in the reactive and non-reactive case are exactly the same. Following [ 21], we see that the internal stress tensor can be re-written in terms of the free energy per unit volume of the undeformed solid, WI(;Cs;Js b) =Js I, as TI=I+Te; (57) where the tensor Teis given by, Te=2 JsFs@WI @CsFT s: (58) Hence, inequality (52) reduces to mrx NX a=1^fa(vavs)NX b=2^ b +@ I @ b+ b 2(vbvs)2 ^ s0: (59) Pore pressure: Stress tensors in uids moving through porous media are typically written so as to include a component called pore pressure. A mixture theoretic perspective of pore pressure can be found in [21] and [28]. In the current context, the pore pressure for uid bcan be de ned as, pb=+@ I @ b: (60) This pressure consists of two components, a hydrostatic or base pressure, given by ph=and a component-speci c pressure, given by pbc=@ I @ b, so that for all uids b, pb=ph+pbc: (61) 12 Also, only the hydrostatic component of pressure acts on the solid, so that for the solid pressure, ps we can write, ps==ph: (62) Thus, we see that the hydrostatic component appears in the stress tensor of all the components in the mixture, including the solid, whereas pbcappears only in the stresses in uid components. In our current model, we will be assuming that all the uids are equally wetting and that the pores are perfectly saturated with the liquids. Hence, we assume that capillary e ects can be ignored, so that pbc0, which gives, pbp=ps; (63) forb= 2;3;:::;N . In other words, we assume that the contribution of hydrostatic pressure (which acts on the solid as well as all the uids) in the pore pressure is more signi cant than the contribution of the component-speci c part. This condition is also achieved when we are considering a mixture of miscible uids, for which interfacial uid- uid forces are absent. Substituting Eqs. (49), (55), (60), and (63) into Eq. (26), we get, b@vb @t+rxvbvb = brxp+@ b @rx+bbb+^fb; (64) for all uid components, b= 2;3;:::;N . Similarly, the balance of momentum equation for the solid and the mixture can be obtained from Eqs. (26), (30), (54a), (56), and (63). Equilibrium and linearization assumptions: With the de nition of pore pressure (Eq. (60)), the inequality (59) shows that the state given by rx=0; wb=0; b = 2;3;:::;N(65) and pb= 0; b = 2;3;:::;N ps= 0;) (with assumption in Eq. (63): p= 0), (66) corresponds to the state of thermodynamic equilibrium. In the absence of chemical reactions ( ^ b0), this system reduces to that described by Bowen, [ 21]. A state of equilibrium for a non-reactive porous medium ow exists simply when the conditions in Eqs. (65) are satis ed, and the additional constraint in Eq. (66) is essential only when chemical reactions are involved. Thus, the conditions in Eqs. (65) correspond to the state of thermo-mechanical equilibrium, and the addition of the constraint in Eq. (66) corresponds to the state of chemo-thermo-mechanical, or thermodynamic equilibrium of the system. At the state of thermodynamic equilibrium, we observe that the following conditions are satis ed, m(;0;Cs;Gs; b;db;0) =0; (67) ^f a(;0;Cs;Gs; b;db;0) =0: (68) We assume that departures from the equilibrium state are small, so that the linearized assumptions formand^fapresented in [ 21] hold. We observe that the balance of momentum equations in the reactive and non-reactive cases are exactly the same, except for the addition of the term NP a=1^cavato the momentum balance equation for the mixture, as seen from Eq. (30). Also, since the constitutive 13 forms of ^fa, and in turn of ^paas given by Eq. (49) have been evaluated only for a= 2;:::;N , the interaction force for the solid, ^pswill be evaluated using the last of Eqs. (29). We now introduce our nal modelling assumption, that all the processes (chemical reactions, deformations etc.) in our mixture are isothermal. Moreover, we also assume that the mixture is in thermal equilibrium, which means that all the components in our mixture are at a constant, uniform temperature at all times. Thus, we can disregard the energy equation and the constitutive restrictions imposed on m. 3.2 Balance of momentum equations Under the assumptions stated so far, the Eulerian form of momentum balance equations for the uids, solid and mixture are given, respectively by, b@vb @t+rxvbvb = brxp+bbb 2 bb s(vbvs); (69a) s@vs @t+rxvsvs = srxp+sbs+NX b=2 2 bb s(vbvs) +rxTeNX a=1^cava(69b) NX a=1a@va @t+rxvava =rxp+b+rxTeNX a=1^cava: (69c) Qausi-static approximation: We are interested in applications wherein the accelerations and thus, inertia e ects are negligible - which we call quasi-static processes. Hence, for such applications, the momentum conservation equations reduce to, 0= brxp+ b bbb bb sv lb; (70a) 0= srxp+ s sbs+NX b=2 bb sv lb+rxTeNX a=1^cava (70b) 0=rxp+rxTe+NX a=1( a aba^cava): (70c) We have used Eqs. (5)and(16) in writing Eqs. (70a) and(70b) . Equations (70a) {(70c) can be written in the ALE framework as, 0=Js  bFT srXsp+ b bbb bb sv lb ; (71a) 0=Js sFT srXsp+Js s sbs+JsNX b=2 bb sv lb+rXsPeJsNX a=1^ca av la (71b) 0=JsFT srXsp+rXsPe+JsNX a=1 a aba^ca av la ; (71c) where Pe=JsTeFT s, the rst Piola-Kirchho elastic stress tensor. For our analysis, we assume the free energy of the mixture to be given by, I= sWs(Cs;t); (72) 14 whereWsis the free energy of the pure solid in the current con guration. Assuming the solid to be hyper-elastic, TeandPeare then given by, Te= 2 sFs@Ws @CsFsT; (73) Pe= 2Js sFs@Ws @Cs: (74) We model our problem based on the degradation of biopolymers. In particular, we follow the hydrolytic degradation of polyesters as described in [ 6] to model our problem. The reduced hydrolytic degradation reaction as described in [ 6] involves three constituents: the polymeric solid, the product of degradation (monomer) and water, which is present in abundance in the system and acts as a base uid. Thus, we model our numerical problem as a three-component system comprising of these three components. We assume the degradation product to have a uid-like behavior, and will be referred to as the `monomeric uid' in the rest of the paper. Following the notation used so far, the properties pertaining to the solid, monomeric and base uid will be denoted by the subscriptss, 2 and 3 respectively. 4 Finite Element Implementation In this section we present a FEM based numerical framework, based on the reduced mixture theoretic model derived in the previous section. Also, due to the saturation condition, Eq. (6), the volume fraction of one of the components in the mixture can be evaluated if the volume fractions of all other components in the system are known. Hence, we eliminate the base uid volume fraction, 3, from the set of primary variables and for each time step, evaluate 3from Eq. (6)and solutions for sand 2. As is customary and for the ease of application of boundary conditions, we resort to the ALE framework. 4.1 Intial and boundary conditions and ALE strong form 4.1.1 Initial conditions. The initial conditions for the density and velocity elds are set to a(x;0) = 0 a(x); us(x;0) =u0 s(x); v lb(x;0) =v0 lb(x); (75) fora=s;2 andb= 2; 3,8x2 t. Also, v0 s=dug s dt x xed. 4.1.2 Boundary conditions. Referring to Fig. 1, the boundary @ tis partitioned into subsets (@ t)D dand(@ t)N d, such that (@ t)D d\(@ t)N d=;and(@ t)D d[(@ t)N d=@ t, whered= s, 2, and us. The superscripts D andNstand for `Dirichlet' and `Neumann,' respectively. For all t >0, the following boundary conditions are then admitted for our problem, d(x;t) = g d(x;t); forx2(@ t)D d; (76) us(x;t) = ug s(x;t);forx2(@ t)D us; (77) Ts(x;t)ns(x;t) =sg(x;t); forx2(@ t)N us; (78) 15 whered=s;2 and where g d,ug s, and sgare the prescribed volume fractions, solid displacement, and traction distribution for the solid. The boundary conditions involving the ltration velocity will be discussed in the next paragraph. Finally, we note that Ts:=pI+Te(79) represents the total stress in the mixture. As mentioned earlier, our nal formulation is written in terms of an ALE framework, for which it is essential to pull back all equations and boundary conditions to the solid reference con guration. In this context, we also partition the boundary @ sinto (@ s)D dand (@ s)N d, such that these partitions are simply obtained by pulling back the corresponding partitions (@ t)D dand(@ t)N dto the solid reference frame. The corresponding boundary values can then be assigned onto these boundary partitions on the basis of Eqs. (2) and (3). Boundary conditions for the ltration velocity: As is well-known (cf., e.g., [ 29,30]), equa- tions like those governing Darcy ow do not support the speci cation of boundary values for the ltration velocity. The only control on boundary values of the ltration velocity pertains to its normal component. Hence, the boundary condition involving the ltration velocity is as follows: for allt>0 v lbns=vn lb(x;t) for x2(@ s)Dnv lb; (80) where (@ s)Dnv lb[(@ s)N v lb=@ s, and where vn lbis a prescribed scalar function. If the boundary is impermeable, then vn lb= 0. With a slight abuse of terminology, we will still categorize the boundary condition in Eq. (80) as being a Dirichlet boundary condition limited to the normal component of the velocity. The weak implementation of Eq. (80) will be discussed in the following section. Finally, when it comes to the admissible traction boundary conditions pertaining to the ltration velocity equations, said traction can only act along the normal to the boundary. Pressure constraint: Similar to incompressible ow problems (cf. [ 22,29,30]), under pure Dirichlet velocity boundary conditions the pressure eld must be constrained in order to obtain a unique solution (the pressure solution is unique up to an additive constant). We therefore adopt the following (standard) constraint on the pressure: Z spd s= 0: (81) Thus, two types of boundary conditions are investigated in our work: 1.Dirichlet boundary condition for uson all boundaries along with the above pressure constraint. 2.Neumann boundary condition for the mixture in which the boundary traction shown in Eq.(78) is prescribed on a part of the domain boundary and a Dirichlet boundary condition is prescribed on the remaining boundary. For both the boundary conditions, Eq. (80) for the ltration velocity is enforced on the entire boundary (i.e. (@ s)Dnv lb@ s). 16 4.1.3 Strong form, Eulerian framework In an Eulerian framework, the strong form of the problem we consider is: Given the •body force elds, bs;b2;b3: s[0;T]!Td, pure component densities, s, 2, 32R+; •prescribed boundary values, g s:(@ t)D a(0;T]![0;1], g 2:(@ t)D 2(0;T]![0;1], ug s:(@ t)D vs(0;T]!Td,sg:(@ t)N vs(0;T]!Td,vn l2: (@ t)Dnv l2(0;T]!R, vn l3: (@ t)Dnv l3(0;T]!R, •prescribed initial conditions, 0 s(x) : 0![0;1], 0 2(x) : 0![0;1],u0 s(x) : 0!Td, v l2(x) : 0(x)!Td,v l3(x) : 0!Td, •and constitutive equation for Te, Eq. (73), nd s: t(0;T]![0;1], 2: t(0;T]!(0;T],u0 s: t(0;T]!Td,v0 s: t(0;T]!Td, v l1: t[0;T]!Td,v l2: t(0;T]!Tdandp: t(0;T]!R(where 0 tjt=0 t), such that,8Xs2 tand8t2[0;T], Eqs. (19) and(70a) are satis ed for a=s;2 andb= 2;3 along with Eqs. (24),(70c) ,(6), and the initial conditions (Eqs. (75)) as well as the prescribed boundary conditions: Eq. (76) and Eq. (80), with either •Eqs. (78) and (77), or, •Eq. (77) and the constraint in Eq. (81). 4.1.4 Strong form, ALE framework. Referring to Eqs. (1),(2), and (3), we de ne the map s: s! t, so that we have x= s(Xs;t), then the ALE strong form of the problem reads - given the •body force elds, bs;b2;b3: s[0;T]!Td, pure component densities, s, 2, 32R+; •prescribed boundary values, g s:(@ t)D a(0;T]![0;1], g 2:(@ t)D 2(0;T]![0;1], ug s:(@ t)D vs(0;T]!Td,sg:(@ t)N vs(0;T]!Td,vn l2: (@ t)Dnv l2(0;T]!R, vn l3: (@ t)Dnv l3(0;T]!R, •prescribed initial conditions, 0 s(x) : 0![0;1], 0 2(x) : 0![0;1],u0 s(x) : 0!Td, v l2(x) : 0(x)!Td,v l3(x) : 0!Td, •and constitutive equation for Pe, Eq. (74), nd s: s(0;T]![0;1], 2: s(0;T]!(0;T],u0 s: s(0;T]!Td,v0 s: s(0;T]!Td, v l1: s[0;T]!Td,v l2: s(0;T]!Tdandp: s(0;T]!R(where 0 tjt=0 s), such that,8Xs2 sand8t2[0;T], Eqs. (21) and(71a) are satis ed for a=s;2 andb= 2;3 along with Eqs. (25),(71c) ,(6), and the initial conditions (Eqs. (75)) as well as the prescribed boundary conditions in Eqs. (76) and (80), with either •Eqs. (78) and (77), or, •Eq. (77) and the constraint in Eq. (81). 17 4.2 Functional settings and weak form We de ne the time dependent spaces for the solid displacement, velocities, pressure and volume fractions as, Vus=fus us2L2( t)nd;rxus2L1( t)ndnd;us=ug son (@ t)D vsg; Vv lb=fv lb v lb2L2( t)nd;rxv lb2L2( t)ndndg; Vp=L2( t); V a=f a a2(0;1]\L1( t);rx a2L2( t); a= g aon (@ t)D ag;(82) forb= 2;3,a=s;1 andnd= 2;3 is the dimension of the Euclidean vector space. The corresponding Galerkin-weighting spaces for the variables are then given by, Vusw=fus us2L2( t)nd;rxus2L1( t)ndnd;us=0on (@ t)D vsg; Vv lbw=Vv lb; Vp w=Vp; V aw=f a a2(0;1]\L1( t);rx a2L2( t); a= 0 on (@ t)D ag;(83) Weak implementation of the boundary conditions for the ltration velocity. The boundary condition in Eq. (80) was implemented using a modi ed version of the method of Lagrange multipliers, as outlined by Babu ska in [ 31]. We introduce an auxiliary boundary variable, b2H1 2(@ t) and its associated test function, ~b2H1 2(@ t). The impermeability boundary condition for v lbusing the method of Lagrange multipliers is then given by, Z @ th ~b(v lbnvn lb) +b(~v lbn)i d = 0; (84) where nis the unit outward normal to the boundary @ t. In the ALE framework the condition is written as, Z @ sh Js~b(F1 sv lbnsvns lb) +b(F1 s~v lbn)i ds= 0; (85) where nsis the unit outward normal to the boundary @ s. The expected value of the boundary eldbwill be discussed after the weak formulation is presented. 4.2.1 Weak form in the Eulerian framework. Given the same data as in the strong formulation, nd us2Vus,v l12Vv l1,v l22Vv l2,p2Vp, 22V 2and s2V ssuch that8~us2Vusw,~v l22Vv l2w,~v l32Vv l3w,~p2Vp w,~ 22V 2wand ~ s2V sw, Z th ~ 2 @t 2jx+rxv l2+ 2rxvs^ 2 +~ s @t sjx+ srxvs^ si d t= 0 (86) Z t ~us^c2 2v l2NX a=1 aba (rx~us)p+rx~us:Te d t+Z (@ t)N us~us pnsg d = 0; (87) 18 Z t~p rx(vs+v l2+v l3)NX a=1^ a! d t= 0; (88) Z t ~v l22 sv l2b2 (rx~v l2)p d t +Z @ t (~v l2n)(p +2) +~2(v l2nvn l2) d = 0; (89) Z t ~v l33 sv l3b3 (rx~v l3)p d t +Z @ t (~v l3n)(p +3) +~3(v l3nvn l3) d = 0: (90) As can be seen in Eqs. (90) and(90), the Lagrange multipliers 2and3enforcing the boundary conditions for the ltration velocities have dimensions of pressure. This is not surprising as the pressure is the work conjugate of the divergence of said velocities. Following the arguments presented in [31], the expected solution for the 2and3is the trivial solution everywhere on @ t. 4.2.2 Weak form in the ALE framework. Given the same data as in the ALE strong problem, nd us2Vus,v l12Vv l1,v l22Vv l2, p2Vp, 22V 2and s2V ssuch that8~us2Vusw,~v l22Vv l2w,~v l32Vv l3w,~p2Vp w,~ 22V 2w and~ s2V sw, Z sJsh ~ 2(@t 2jXs+FT s:rXsv l2+ 2FT s:rXsvs^ 2) +~ s(@t sjXs+ sFT s:rXsvs^ s)i d s= 0; (91) +Z s" Js~us NX a=1 aba+^c2 2v l2! JsprXs~us:FT s+rXs~us:Pe# d s +Z (@ s)N vs JspF1 s~usnssg s ds= 0; (92) Z sJs~p FT s:rXs(vs+v l2+v l3)NX a=1^ a! d s= 0; (93) +Z sJs ~v l22 sv l2b2 r Xs~v l2:FT sp d s +Z @ sh ~3 JsF1 sv l2nsvns l2) +Js~v l2(p+2)FT snsi ds= 0;(94) 19 (a)kv l2k (b)kv l3k (c)kpk (d) s (e) 2 (f)kusk Figure 2: Manufactured solutions for di erent variables +Z sJs ~v l33 sv l3b3 r Xs~v l3:FT sp d s +Z @ sh ~3 JsF1 sv l2nsvns l3) +Js~v l3(p+3)FT snsi ds= 0:(95) Here, sg sis the prescribed traction eld on ( @ )s, which is related to sas follows ss=JsjjFT snsjjs: (96) Discrete Approximation: We partition the domain sinto a triangulation h sconsisting of cellsK, with diameter h. The cells Kare chosen to be non-overlapping quadrilaterals in the 2D case and hexahedrons in the 3D case, which span the spatial domain s. The decomposition of the boundary@ sinto Neumann and Dirichlet partitions is also followed by the boundary of h s. The interpolation functions over Kare chosen to be Lagrange polynomials. 5 Numerical Results The weak formulation derived in the previous section was implemented in a FEM framework, and the resultant numerical model was studied for accuracy. 5.1 Model numerical problem For the testing of our numerical scheme, we choose a 2-dimensional square of side length L= 1 m as the domain s. The elastic response of the solid is chosen as Neo-Hookean, with the solid energy 20 Table 1: Parameters Parameter Value ms 3103kg=ms2 2 1103kg=sm 3 1:5103kg=sm s 1103m2Table 2: MS constants Constant Value us 0:001 m v l1 0:1 m=s v l2 0:1 m=s p 1 PaTable 3: constants Case  2 2 s s 1 0.2 0.1 0.5 0.2 2 0.35 0.15 0.3 0.1 3 0.2 0.1 0.2 0.1 density given by, Ws(Cs) =ms 2(tr(Cs)3): (97) The value of msis chosen such that the constitutive response of the solid resembles that of bio- polymers. The viscosities of the uids were chosen to be of the order of the viscosity of water. The exact parameter values used in our analysis are given in Table 1. The numerical scheme is tested using the Method of Manufactured Solutions [ 32]. The functional forms of di erent variables (manufactured solutions) are chosen as, us= ussin2t t0 cos 2x+y L ^{+ sin 2xy L ^| ; (98) v l2= v l2cos2t t0 sin 2x2+y2 L2 ^{+ cos 2x2y2 L2 ^| ; (99) v l3= v l3cos2t t0 cos 2x2 L2 sin 2y2 L2 ^{+ sin 2x2y2 L2 ^| ; (100) p= psin2t t0 sin 2x+y L ; (101) 2= 2+ 2cos2t t0 cos 2x+y L ; (102) s= s+ ssin2t t0 sin 2x+y L ; (103) where ^{and ^|are the base vectors of the adopted Cartesian coordinate system. The values chosen for di erent constants are given in Table 2. The constants  2, 2, sand sare adjusted so as to cover di erent phases of degradation. The three combinations of the volume fraction constants used in our analysis are given in Table 3. Fig. 2 shows representative solution eld distributions for all the variables over the domain s. The solution is obtained using a mesh of size h=L=128, for the case 2 mentioned in Table 2, at time = 0 :7 s. The nature of the solution is qualitatively the same for all the cases, at all times, and for all re nements levels. 5.2 Solver and mesh speci cations The equations presented in this paper constitute a transient di erential-algebraic system, which was solved in COMSOL Multiphysics®. The weak formulation was studied for stability and convergence for the two sets of boundary conditions discussed in the previous section. The primary variables for the numerical study were us;v l2;v l3; s; 2andp. The solid velocity, vswas calculated in a post-processing step, using the kinematic relation: vs=@tus Xs. Implicit time stepping was implemented by the method of lines [ 33,34], using the variable-order variable time step backward 21 di erentiation formulas ( BDF ). The order of BDF was constrained between 2 to 5, whereas a maximum time step size of 1 103s was imposed. The resultant non-linear system of equations for every time step was solved using Newton's method (with no damping). PARDISO ([ 35{37]) was chosen as the linear solver for the resultant linear system of equations. The unit square domain was discretized uniformly into square elements of side length 1 =h, forh= 8;16;32;64;128 m. Piece-wise Lagrange interpolation polynomials were used for all elds. No additional stabilization technique has been implemented for the numerical formulation. 5.2.1 Choice of order of interpolation functions. It was observed that the solvability and accuracy of the unstabilized quasi-static formulation is dependent on the choice of interpolation order ( IO) of the interpolating functions. The numerical scheme did not provide satisfactory results when pressure was interpolated using an IOequal to or greater than that of the vector variables. This phenomenon was observed in the fully dynamic case (the formulation in which the time derivative and convective terms were included in the momentum balance equations) of non-reactive poro-elastic ow, [ 22] as well. However, unlike the quasi-static case presented in their paper, we observe that having equal IOfor pressure and vector elds does not yield convergent solutions to our numerical problem. This condition likely stems from a Brezzi-Babu ska (inf-sup) type condition. In order to get acceptable numerical results, the IOof pressure interpolation functions was set to be one less than that of the vector variables (velocities and solid displacement). We choose IOs of 2 and 3 for the vector variables, and corresponding order of 1 and 2 respectively for pressure. TheIOof the volume fractions was chosen to be either equal to, or less than that of pressure for each of the above choice of IOs for the rest of the variables. It was observed that having a higher order for volume fractions did not produce acceptable numerical results. Thus, we will be discussing three combinations of interpolation polynomials: •Case O-1: Quadratic polynomials for vector variables, linear for pressure and volume fractions, •Case O-2: Cubic polynomials for vector variables, quadratic for pressure whereas linear polynomials for volume fractions, •Case O-3: Cubic polynomials for vector variables and quadratic elements for pressure as well as volume fractions. TheIOof the auxiliary boundary variable bfor the implementation of the impermeability boundary condition for the ltration velocity was chosen to be equal to, or less than the order of v lb, forb= 2;3. Thus, Lagrange linear and quadratic IOs were chosen for bfor case O-1 and linear, quadratic and cubic polynomials for cases O-2 and O-3. We would also like to mention that acceptable results were not obtained when linear interpolation polynomials were used for bwhile using cubic polynomials for ltration velocities (cases O-2 and O-3). However, exact same convergence rates were obtained using either quadratic or cubic polynomials for bfor cases O-2 and O-3. Similarly, the convergence rates obtained using either linear or quadratic polynomials for case O-1 were exactly the same. Hence, we arbitrarily choose to report the results for linear binterpolation polynomials for case O-1 and quadratic for case O-2 and O-3. 5.3 Formulation results We now present results of our FEM implementation of the 2D Manufactured Solution problem discussed earlier. The simulations were run for the parameters and constants discussed in Tables 1 22 3.023.071.95 1.992.74 1.96 2.06 2 2.5 3 3.5 4 4.5 5-12-10-8-6-4-2024 2.502.50.92 0.97 1.01 1.121.01 2 2.5 3 3.5 4 4.5 5-9-7-5-3-1135 3.023.041.97 1.992.83 2.00 2.06 2 2.5 3 3.5 4 4.5 5-12-10-8-6-4-202 2.062.11.01 1.00 2.011.121.01 2 2.5 3 3.5 4 4.5 5-11-9-7-5-3-1135 2 2.5 3 3.5 4 4.5 5-18-16-14-12-10-8-6-4-202 2.382.921.01 1.00 2.022.56 2.00 2 2.5 3 3.5 4 4.5 5-13-11-9-7-5-3-1135Figure 3:L2andH1error norms for case 1 in Table 3, for cases O-1, O-2 and O-3 (from top to bottom respectively.), for the weak formulation along with Neumann boundary condition for us and 2 and the three cases described in Table 3. We observed that even without any stabilization, the numerical scheme did not produce any numerical instabilities for these parameter values. Similar convergence rates were obtained for both the sets of boundary conditions (Neumann boundary condition for us; and global constraint on pressure along with pure Dirichlet boundary 23 2.813.031.97 1.992.62 2.01 2.05 2 2.5 3 3.5 4 4.5 5-12-10-8-6-4-202 2.382.50.92 0.97 1.01 1.121.01 2 2.5 3 3.5 4 4.5 5-10-8-6-4-20246 2.813.031.97 1.99 2.62 2.01 2.05 2 2.5 3 3.5 4 4.5 5-12-10-8-6-4-202 2.052.11.01 1.00 2.001.131.01 2 2.5 3 3.5 4 4.5 5-11-9-7-5-3-1135 2 2.5 3 3.5 4 4.5 5-17-13-9-5-13 2 2.5 3 3.5 4 4.5 5-10-8-6-4-2024Figure 4:L2andH1error norms for case 2 in Table 3, for cases O-1, O-2 and O-3 (from top to bottom respectively.), for the weak formulation along with Neumann boundary condition for us conditions for us), with the impermeability condition implemented using the method of Lagrange multipliers. We rst present the results for the `Neumann boundary condition for us' case. Results for convergence rates are given in Figs. 3, 4 and 5 for the three cases mentioned in Table 3. Convergence results are presented in terms of the logarithms of both the L2norm and the H1 24 2.863.061.95 1.99 2.44 1.96 2.05 2 2.5 3 3.5 4 4.5 5-14-12-10-8-6-4-2024 2.402.50.92 0.97 1.01 1.121.01 2 2.5 3 3.5 4 4.5 5-10-8-6-4-20246 2.873.051.97 1.99 2.71 1.99 2.05 2 2.5 3 3.5 4 4.5 5-14-12-10-8-6-4-202 2.052.11.01 1.00 2.001.121.01 2 2.5 3 3.5 4 4.5 5-11-9-7-5-3-1135 2 2.5 3 3.5 4 4.5 5-18-14-10-6-22 2 2.5 3 3.5 4 4.5 5-10-8-6-4-2024Figure 5:L2andH1error norms for case 3 in Table 3, for cases O-1, O-2 and O-3 (from top to bottom respectively.), for the weak formulation along with Neumann boundary condition for us semi-norm of all the quantities determined in the solution. Such results are presented as `empirical evidence' collected in the numerical experiments. This said, by reporting convergence in the H1 norm we do not mean to suggest that said norm is controlled by the data of the problem for all of the quantities at hand. Proper expectations on the various elds comprising the problem have been 25 2.793.041.96 1.982.1 1.982.18 2 2.5 3 3.5 4 4.5-14-10-6-226 2.512.450.92 0.97 1.011.24 1.02 2 2.5 3 3.5 4 4.5-7-5-3-1135 2.813.051.97 1.99 2.952.032.18 2 2.5 3 3.5 4 4.5-14-10-6-22 2.12.241.01 1.00 2.011.34 1.02 2 2.5 3 3.5 4 4.5-11-9-7-5-3-1135Figure 6:L2andH1error norms for case 1 in Table 3, for cases O-1 and O-2 (top and bottom respectively.), for the weak formulation along with pure Dirichlet boundary condition on all boundaries for usand the global integration constraint on p. stated in the choice of their respective functional spaces. Hence, for example, the ltration velocities are only expected to converge in the L2norm andH(div) semi-norm. Clearly, if the numerically observed data show convergence in the H1semi-norm, H(div) convergence is implied. The convergence behavior of the numerical scheme is observed to be similar at di erent times. Hence, without loss of generality, the norms have been reported at t= 0:3 s. For all the three cases, we observe that near-constant convergence rate of about 2 was obtained for the L2error norms for both the ltration velocities, irrespective of the order of interpolation functions. Similarly, the L2 andH1error norms for pressure approach optimum values of 3 and 2 respectively for the cases O-1 and O-2. These results are comparable with the `quasi-static' formulation case in the non-reactive ow through poro-elastic solid study presented in [22]. TheL2error norm for solid displacement is seen to be above 2.5 in cases O-1 and O-2 for all ranges of volume fractions, whereas the H1error norm is above 2 in these cases. Solid velocity achieves optimum L2norm of about 3 in cases O-1 and O-2 as well. In case of O-3, we observe a relatively erratic behavior of the error norms of solid displacement at higher resolutions. The two additional scalar elds with respect to the non-degrading poro-elastic case, namely the volume fractions, were seen to have optimal L2convergence rates of about 2 for the combination of interpolation polynomials given by cases O-1 and O-2. However, in the case of O-3, the behaviour 26 of theL2error norm is unacceptable for higher resolutions. Thus, linear Lagrange polynomials seem to be the best choice for the volume fraction elds, along with the choice of shape functions discussed earlier for pressure and vector elds. A point worth mentioning is that the numerical solver failed to convergence when the solid volume fraction went above a value of about 0.85 or below 0.1. We would like to note, however, that the case when the solid volume fraction is extremely low (0.1 or less), would correspond to a physical scenario in which most of the solid has already degraded, thus leading to a loss of its initial porous structure. Under these circumstances, we expect that the constitutive equations for Peand Wsdescribed earlier will cease to describe the elastic behaviour of the solid accurately, since the `amount' of solid will be too low in the mixture. The failure of the numerical scheme to converge at very high volume fractions of the solid were also observed in the case of non-degrading mixture of solid and single uid in [22]. As mentioned earlier, similar tests were also conducted using the Method of Manufactured solutions for the second boundary condition, viz. Dirichlet boundary condition for usand the global constraint on pressure. The convergence rates are similar to those obtained in the case of Neumann boundary condition for us. Hence, to avoid repetition, we choose to report the results for just one representative case (case 1) for the `Dirichlet boundary condition for us' case. The results are shown in Fig. 6. We also conducted some numerical experiments to study the e ect of relative magnitude of the ltration velocities with respect to the solid velocity (or displacement) on the numerical outcome. The solid displacement amplitude, uswas held constant at 0 :01 m, while the amplitudes of both the ltration velocities v l2andv l3were identically varied. No changes were made to the rest of the variable amplitudes. The numerical scheme failed to converge for v l2andv l3values of 1 m=sand higher. For such high values of ltration velocities, it would be more appropriate to use the full dynamic formulation (cf. Eq. (26)) along with a suitable stabilization scheme. This is beyond the scope of this paper and would be the focus of our future work. 6 Conclusions A mixture theoretic formulation for the study of degradation of a poro-elastic solid immersed in a base uid is derived. Our transient 3-component system consists of the solid, the product of degradation known as the monomeric uid, and the base uid. The solid is embedded in the base uid, which also initiates the hydrolytic degradation reaction of the solid. The product of degradation reaction remains in the system, and is expected to a ect the rate of degradation as well. The constitutive equations for the three component system are derived in accordance with the principals of continuum mixture theory, so as to satisfy the second law of thermodynamics for the mixture. The model is easily extendable to include more number of uids, or reaction products as well. The equations are applicable for ow of multiple incompressible Newtonian uids owing through a degrading incompressible poro-elastic solid. The simpli ed, ltration velocity-based, quasi-static ALE formulation was implemented into an FEM based numerical model. The formulation was studied for stability and convergence rates using the method of manufactured solutions. The numerical model is stable over a wide range of practically relevant values of ltration and solid velocities, without any additional stabilization. Future work will focus on implementing the full analytic model, including the convection terms into a numerical model, with the possibility of including stabilization for the fully dynamic case. The eventual goal is to apply the model to studies pertaining to real bio-degradable sca old materials like CUPE when inserted at real target cites. 27 Acknowledgments This work was supported by the U.S. National Science Foundation (CMMI 1537008).