# Analysis and numerical solution of nonlinear Volterra partial integrodifferential equations modeling swelling porous materials

we establish under a given set of assumptions for the initial-boundary value problem. Additionally, a special case of the VPIDE is reduced to an ordinary differential equation via a derived similarit)' variable and solved. In order to solve the full VPIDE we derive a novel approach to constructing pseudospectral dif ferentiation matrices in a polar geometr}' for computing the spatial derivatives. By construction, the norms of these matrices growr at the optimal rate of C(i ¥2), for N-by-N matrices, versus 0(N4) for conventional pseudospectral methods. This smaller norm offers an advantage over standard pseudospectral methods when solving time-dependent problems t hat require higher-resolution grids and, potentially, larger differentiation matrices. A method-of-lines approach is em ployed foi the time-stepping using an implicit, fifth-order Runge-Kutta solver. After we show how to set up the equation and numerically solve it using this method, we show and interpret results for a variety of diffusion coefficients, permeability models, and parameters in order to study the model's behavior. This abstract accurately represents the content of the candidate's thesis. I recommend its publication. signed o4fw6mmnn \y Lynn Bennethum iv

DEDICATION This thesis is dedicated to my wife Sharyl and my son Cole for the inspiration they give me everyday to be a better husband, a better father, and a better man.

ACKNOWLEDGMENT There are several people I need to acknowledge for the help they provided in the completion of this thesis. First, this thesis would not have been possible without the patience and insights of Lynn Schreyer-Bennethum. Lynn taught me the value of using physical insight to guide mathematical analysis. Kristian Sandberg introduced me to Pseudospectral methods and showed me how to think creatively as a numerical analyst. His trouble-shooting acumen saved me on many occasions and I cannot thank him enough for being so generous with his time despite running his own company and raising a family. I would also like extend my thanks to Bill Briggs for introducing me to Lynn and nudging me along as I came perilously close to leaving the graduate program. To Julien Langou, who helped me think more deeply about linear algebra. To Leo Franca for helping to bring me t o UC Denver. To Alexander Engau, for agreeing to be on my committee despite being in the unenviable position of being asked at the last minute. To Mat t hew Nabity who has been my compadre since the start of this journey and has become a life-long friend. To John Stinespring who, simply put, makes my life better. To Eric Sullivan and Kannanut Chamsri, my academic brother and sister, for patiently listening to my research talks and asking poignant questions. To Jinhai Chen for helping make my arguments more mathematically sound. To Jeff Barchers for believing in me enough to hire me at SAIC. To Richard Quanst rum who taught me, by his example, what it means to be a real man, I carry in my heart your pride in my accomplishments as a teacher, husband, and father - you will always be my coach. To Angela Beale and Lindsay Hiatt who kept me registered and up-to- date on my requirments, I would not have made it through a semester without them. To all the people whom I met over my years of education who helped me along the way by their kindness, advice, instruction, and encouragement. Finally, I must thank my family for their loyalty and unending support. My mother, Eleonore, who nurtured my creativity and encouraged me to "reach for

the stars," my brother Robbie who taught me the value of hard-work, my sisters Barbara and Cheryl who always stood by me growing up and who taught me to be a gentleman, my brother Lonnie who taught me that there is more to life than work (like music!), and my father Leonard who taught me to be self-reliant.

CONTENTS Figures xi Tables xvii Chapter 1. Introduction 1 1.1 Previous Work 5 1.2 Thesis Outline 7 2. Fundamental Continuity Equation 10 2.1 Brief Overview of Hybrid Mixture Theory 10 2.1.1 Macroscale Field Equations 16 2.2 Derivation of the Continuity Equation 18 2.3 Darcy's Law 20 2.4 Non-dimensionalization of the Model 25 2.5 Continuity Equation 27 2.6 Model Interpretation 28 2.6.1 Deborah Number 32 2.6.2 Integral Coefficient 35 2.7 Discussion 36 3. Analysis of the Volterra Partial Integrodifferential Equation 39 3.1 Formulation of the Boundary Value Problem 39 3.2 Existence of a Solution 46 viii

3.3 Uniqueness of a Solution 49 3.4 Analytic Solutions and Reductions of the IBVP 56 3.4.1 Similarity Reduction of the IBVP 56 3.4.2 Viscous Case 59 3.4.3 Flory-Huggins Model 62 3.5 Discussion 65 4. Eigen-decomposition Pseudospectral Method 67 4.1 Pseudospectral Methods 71 4.2 Eigen-decomposition Pseudospectral Method 76 4.2.1 Derivation 76 4.2.1.1 Truncating the Sum 77 4.2.1.2 Approximating the Integral with a Quadrature 78 4.2.2 Error Analysis 80 4.2.3 Rank Completion 84 4.3 EPS Construction in Polar Coordinates 89 4.3.1 Construction 90 4.3.2 Numerical Examples 93 4.4 Discussion 95 5. Numerical Solution of the Partial Integrodifferential Equation .... 102 5.1 Numerical Method for Solving the Swelling Equation 102 5.1.1 Spatial Discretization 104 5.1.2 Time-stepping Methods 107 5.1.2.1 Semi-analytic Integration Rule Formulation 109 5.1.2.2 Method-of-Lines Formulation 110 ix

5.1.2.3 Pouzet Volterra Runge-Kutta Formulation I l l 5.1.2.4 Examples I l l 5.2 Numerical Solution of the IBVP 114 5.2.1 Two-dimensional Example 115 5.2.2 One-dimensional Examples 115 5.2.2.1 Flory-Huggins Model 116 5.2.2.2 Diffusion Coefficient Comparison 118 5.3 Discussion 119 6. Model Sensitivity Analysis 126 6.1 Diffusion and Permeability Models 127 6.2 Parameter Sensitivity 130 6.2.1 Moisture Content Curves 131 6.2.2 Viscoelastic Stress Curves 134 6.3 Discussion 135 7. Conclusion and Future Work 144 7.1 Model Analysis, Validation, and Extension 145 7.2 Generalizing the Applicability of the EPS Method 146 7.3 Generalizing the Existence and Uniqueness Proof for the IBVP . . 147 7.4 Extending the Applicability of the Drug Delivery Model 147 Appendix A. Derivation of Darcy's Law 149 B. Derivation of a Pouzet Volterra Runge-Kutta Method 156 References 159 x

FI GURES Figure 1.1 Porous polymer matrix. Photo taken from plc.cwru.edu 4 2.1 Averaging, Local Coordinates 13 2.2 Stress (a) versus strain (e - not to be confused with the volume fraction el) curve for linear elastic (left) and linear viscoelastic (right) materials 30 2.3 Stress versus time where the loading takes place over the time interval t < ti and the unloading takes place over the interval t>t\ 30 2.4 Temperature dependence of rate between transition states on tem perature 34 3.1 The model geometry is a cylinder but we assume angular and az- imuthal symmetry thus reducing the domain to a rectangle Q with boundary r = Tx U T2 40 3.2 Integral of temporal transform, B(uS) 48 3.3 Similarity solution to the IBVP (3.55) 61 3.4 Solution to the Flory-Huggins model (3.58) for a variety of times. . 64 XI

4.1 Trial function comparison for spectral, finite-difference, and finite- element methods. Spectral methods are characterized by one high- order polynomial for the whole domain, finite-difference methods are characterized by multiple overlapping low-order polynomials, and fi nite element methods are characterized by non-overlapping polyno mials with compact support - one per subdomain 69 4.2 Runge's example for interpolating f(x) = r5 L+i with evenly spaced nodes. As the number of evenly-spaced nodes (degree of the interpo lating polynomial) increases, the interpolating polynomial oscillates near the boundaries 73 4.3 Interpolating with Chebyshev nodes ameliorates Runge's phenomenon. As the number of Chebyshev nodes increases, the interpolation error decreases 74 4.4 L2 error comparison of the Chebyshev PS method versus second and fourth order finite difference methods applied to the second derivative of cos(7rx) over [—1,1] 75 4.5 Sparsity pat t ern for orthogonality test of eigenfunctions in a Carte sian geometry with Neumann boundary conditions. The rank of the differentiation matrix is Nc = N = 256 85 4.6 Sparsity pat t ern for orthogonality test of eigenfunctions in a Carte sian geometry with Neumann boundary conditions. The rank of the differentiation matrix is Nc = 139 (roughly 0.54./V) 86 xn

4.7 Sparsity pattern for orthogonality test of eigenfunctions in a Carte sian geometry with Neumann boundary conditions for the completed operator. The rank of the differentiation matrix is now TV — 256 (TVC = 139, roughly 0.547V) 88 4.8 Comparison of the L°° norm relative error resulting from comput ing the second derivative of the functions {s i n( 5 | 5:( x+ 1))} _ using the EPS method with Gauss-Legendre quadrature nodes and weights versus the Standard Construction using Cheybshev-Lobatto nodes. Rank-completion versus a reduced rank construction the EPS method is also compared 97 4.9 Given TV, use the EPS method to construct the polar Laplacian in creasing Nc from 1, 2,..., TV and computing the error | | o^ — o£||oo for each set of TVC eigenvalues, {a\,a^}nc=l. The colorbar is log-scaled where the dark regions indicate smaller error and the light regions indicate larger error 98 4.10 Eigenvalues of Ar constructed via Chebyshev collocation versus the exact eigenvalues, \an\2. The eigenvalues of Ar are the zeros of the zero-order Bessel function of the first kind, Jo(ctn) — 0 for all n = 1,2,3, 99 4.11 Chebyshev expansion coefficients versus Fourier-Bessel expansion co efficients for the Gaussian pulse centered at r = 1/2 99 xm

4.12 The Poisson equation example comparing the completed EPS con struction to the standard construction using Chebyshev collocation: (a) The residual error, | |/ — Aru||oo. (b) The error \\u — A"1/^, and (c) L°° condition number 100 4.13 The maximum stable time step At for solving the radial part of the heat equation on a disk using the explict Runge-Kutta 4 solver. Here N denotes the size of the problem 101 5.1 The EPS differentiation matrix for the cylindrical Laplacian (5.8) using Nr = Nz = 64 nodes 108 5.2 Initial liquid volume fraction el(r, z, 0) with elmm = 0.1 and e'max = 0.9.116 5.3 Liquid volume fraction, el(r,z,t), plots over the cylindrical cross- section Q with a Kozeny-Carman permeability, K(el) = J^-J, linear diffusion coefficient, D(el) = el, and model parameters \x = 0.01 and r = 1. The grid size is Nr x Nz = 64 x 64 with Ncr = Ncz = 42(0.65A^r) and a conventional (explicit) RK4 time-stepper was used with constant time-step At = 10~6. Solutions are shown at (a) t = 0, (b) t = 0.2, and (c) t = 0.4 123 5.4 Liquid volume fraction el found by solving (5.36) using an MOL approach with the EPS discretization and MATLAB's odel5s for time-stepping 124 5.5 Relative L°° error comparing Chebyshev collocation to the EPS con struction over t € [0,10] for the Flory-Huggins model (5.36) 124 xiv

5.6 Liquid volume fraction, £l(r,t), plots over the radial grid with a Kozeny-Carman permeability, K(el) = ^ r and model parameters H = 0.1 and r = 1. The grid size is Nr = 750 with Ncr = 450(0.607Vr) and a variable step-size, 5th order implicit time-stepper was used. Solutions are shown for t

Normalized moisture content curves M/M^: (a) Fix KS = 0.1 and compare the moisture content for De = {0,0.01,100}, (b) Fix KS = —0.1 and compare the moisture content for De = {0,0.01,100} . . 140 Normalized moisture content curves M/M^: (a) Fix De = 0.01 and compare the moisture content for K,S = {—0.01,-0.1}, (b) Fix De = 100 and compare the moisture content for KS = {—0.01, —0.1}, (c) Fix KS = —0.1 and compare t he cases De — {0.01,1,100}, (d) Fix KS = —0.01 and compare the cases De = {0.01,1,100} 141 Viscoelastic stress 142 Viscoelastic stress 143 xvi

TABLES Table 2.1 Units for the dimensional parameters of the Darcy model (2.50). . . 27 2.2 Characterization of the Deborah number, De 35 2.3 Characterization of KS 36 5.1 Values of the log | | en(i )| | 00 error at times t\ = 0.25i/, t-2 = 0.50£/, t3 — 0.75tf, and t\ — tj comparing RK4 and ODE45 versus SAI for Example 1 120 5.2 Values of the log | | en(i )| | 00 error at times t\ = 0.25£/, t-2 = 0.50?;/, ^3 = 0.75i/, and t± = tj comparing RK4 and ODE45 versus SAI for Example 2 121 5.3 Values of the log He^^Hoo error at times t\ = 0.25/;/, t2 = 0.50i/, t3 = 0.75t/, and t\ = tj comparing RK4 versus PVRK4 versus SAI for Example 2 122 6.1 Characterization of «;s 130 6.2 Characterization of the Deborah number, De 131 B.l Butcher diagram 158 xvn

1. Int roduct i on Generally speaking "drug delivery systems" are devices used for the pro cess of administering pharmaceutical compounds to achieve therapeutic effects. These systems can be divided into two major types: traditional delivery systems and controlled released systems. Traditional delivery systems are characterized by their immediate release of the drug which leaves absorption to be controlled by the body's ability to assimilate the drug concentration into different body tissues such as the blood. Drug concentrations from these systems typically undergo an abrupt increase followed by an abrupt decrease. Controlled release drug delivery systems, on t he other hand, are formulated t o modify t he drug release profile, absorption, distribution, and elimination for the benefit of im proving therapeutic efficacy, safety, patient convenience, and compliance. These systems are characterized by their maintenance of drug concentration to target tissues at a desired level for prolonged time periods. Depending on the release behavior, controlled release systems can be sub divided into three categories: passive pre-programmed, active pre-programmed, and active self-programmed. The models we will study focus on passive pre programmed release where the rate is predetermined and dependent upon the drug's release kinetics as it interacts with the body. We leave this restriction vague for the moment but will further restrict our focus as we derive the mathe matical model and impose simplifying assumptions. In lieu of definitions for t he remaining two categories, which we offer for the sake of curiosity and complete- 1

ness, we provide examples. An example of an active pre-programmed device is an insulin pump for type 1 diabetes patients where the insulin is delivered as a continuous flow to the body through a short t ube with a needle at the endpoint t hat is inserted under the skin, usually in the abdomen. The third cat egory, active self-programmed, is exemplified by devices t hat combine (provide dat a interfacing) insulin pumps with glucose monitors or, as a simple example, morphine drips. It should be evident by the definitions and descriptions given above t hat drug delivery has become more complex moving from uncontrolled release to sustained (controlled) release and programmable (controlled) release. It has also become more specific as we have witnessed a move from systemic delivery to organ (as in the case of cancer radiation therapy) and cellular targeting (as in the case of tumor-targeting delivery systems). As is the case in many fields, a feedback loop originates where increasing complexity in the technology drives the requirement for mathematical modeling and simulation and vice versa. Mathematical modeling of drug delivery is in its inchoate stages, intimated by the reviews given in [41, 60, 76], but has become an area of active and growing research. We will focus our attention on controlled drug delivery systems, also known as controlled release systems (CRS). Furthermore, we will limit the scope of our study to passive pre-programmed CRS and will subsequently drop this modifier. Therefore the CRS systems t hat fall under our scrutiny are those t hat are fabri cated by embedding a drug in a hydrophilic (water-loving) polymer matrix such as hydroxypropylmethylcellulose (HPMC) [76]. "Matrix" in this context refers 2

to the three-dimensional network containing the drug and other substances re quired for controlled release. The matrices can be prepared in several ways. Two examples include (1) mixing the powdered drug with a solvent, excipient (inactive ingredients), and pre-polymer before placing it in a polymerisation re actor or (2) preparing the matrix in advance and then putting it in contact with a highly concentrated drug solution able to swell the matrix whereby the solvent is removed afterwards (solvent swelling technique) [41]. A detailed exposition on polymeric matrices is beyond the scope of this thesis but there are several references available [41, 76]. However, we will briefly describe the properties of polymer matrices useful for controlled release systems. First, polymers are viscoelastic. By the name, it should be apparent t hat viscoelastic materials possess bot h viscous and elastic properties. Viscosity is the characteristic of a liquid t hat makes it resistent to flow. For example, honey is more viscous t han water. Elasticity is the charac teristic of a material t hat returns it to its initial state (instantaneously) after the external forces that deformed it cease. For example, a dry sponge is more elastic t han pumice. Hence, a viscoelastic material has the characteristic t hat it initiates returning to its original state after the external forces t hat deform it cease, but the elasticity is inhibited by the viscous behavior which resists the change in motion. Second, the polymer matrix is a porous material as seen by the polymeric matrix in Figure 1.1. Polymeric matrix systems can be classified according to their porosity (macroporous, microporous, and non-porous). We will focus our attention to the macro and micro porous systems where pore sizes range in 3

size from 0.1 — \\im and 50 — 200A, respectively. Drug release kinetics will be affected by polymer swelling, polymer erosion, drug dissolution/diffusion, drug distribution (inside the matrix), drug/polymer ratio, and system geometry (cylinder, sphere, etc.). The matrix systems under consideration will be dry and compressed without any fluid inside. We consider the gas filled pores to be part of the liquid phase (fluid and gas). Therefore in our models the initial conditions will contain a small percentage of liquid phase. Once the drug delivery device Figure 1.1: Porous polymer matrix. Photo taken from plc.cwru.edu. is immersed in biological fluid, the surrounding fluid pressure coupled with the hydrophilic properties of the polymer drives the fluid to penetrate the polymer matrix as described by Darcy's law for fluid flow (which we will describe in more detail in Section 1.1 and Chapter 2); the novel form of Darcy's law derive in [73] contains a constitutive equation modeling the viscoelastic properties of the polymer. The polymer swells as dictated by its viscoelastic properties and the pores in the matrix enlarge until they are of a size necessary for the drug to escape. The drug then diffuses into the surrounding fluid diminishing the drug concentration levels present in the delivery device. 4

1.1 Previ ous Work The equations considered in this thesis are nonlinear Volterra partial inte- grodifferential equations (VPIDE) derived by Weinstein and Bennethum [73, 74] and Singh et. al. [63]. These models are derived using Hybrid Mixture The ory (HMT) which involves upscahng field equations which govern the motion of materials and include the conservation of mass, conservation of linear and angular momentum balance, and conservation of energy from the microscale to a larger scale via volume averaging. Restrictions are then obtained on the form of the constitutive equations by using the second law of thermodynamics, also formulated as a field equation, at the large scale. A variable t hat results from upscahng via HMT is the volume fraction, a ratio of the volume of a particular phase over the sum of the volumes of all phases. For example, in a two-phase mixture with a liquid phase (indicated by I) and a solid phase (indicated by s), the liquid volume fraction, denoted e(, is the ratio of the volume of the liq uid phase (V1) to the sum of the volumes of bot h phases (Vs + V1). That is, e' = v*+v< • ^ i e n c l md volume fraction, e(, is the dependent variable considered in the models we study in this thesis. Note t hat if there are only two phases, s and /, we have el + e& = 1. In the case of the VPIDE considered, the mass conservation equation is upscaled and coupled to a novel form of Darcy's law [73]. Darcy's law is a constitutive equation relating the flux of a fluid field to a change in pressure [29]. As such it is often used to describe the flow of a fluid through a porous medium. In the case of the VPIDE derived by Weinstein [73], time rates of change of el are included as a constitutive variable. That is, Weinstein included the material 5

time derivatives of the liquid volume fraction, e , for m = 1, 2,... , p where p is determined by how long of a time history we want to include in our model. This temporal history accounts for the viscoelastic stress and the model for the stress assumes that the strain effects are cumulative, hence it contains an integral. That the integral term contains time derivatives of the volume fraction, not solid-phase strain, is justified in [51] where it is shown that at moderate to high fluid contents the normal components of the strain tensor are related to the volume fraction of the solid phase, es. Moreover, since the solid and liquid phases are all that comprise the material, we have that el = 1 — es. Hence, at the fluid content levels considered in these models, the time derivative of strain may be replaced by time derivative of liquid volume fraction. As a result of this novel form of Darcy's law, the VPIDE models liquid penetration into the drug- delivery device inhibited by the viscoelastic properties of the polymers. That is, the VPIDE can be viewed as the sum of a liquid penetration (nonlinear diffusion equation) term and a viscoelastic term (integral equation). An equation similar to the one derived in [73] was derived by Singh et. al. [63]. Singh also numerically solved the VPIDE using a finite-element method [64, 65]. However, the solver produced spurious oscillations that resulted in non-physical values for the volume fraction in [65]. Even though this work was completed in 2003, to our knowledge, it has never been revisited. Moreover, the equation in the drug delivery context [73] has not been solved. Models developed in the pharmaceutical literature that take into account the viscoelastic properties of the polymer have been phenomenological [12, 47, 60, 61, 62, 76] and lack the physics of the models derived by Weinstein [73]. In 6

addition, these models are typically linear partial differential equations, such as Fick's second law of diffusion in cylindrical coordinates [76], t hat can be solved either analytically using separation of variables or numerically using elementary methods [60]. There are diffusive models taking into account viscoelastic relax ation in polymers t hat have been developed. In particular Cohen and Whi t e [24] extended the work of Thomas and Windle [67] in modeling sharp fronts due to diffusion and viscoelastic relaxation in polymers. It should be noted t hat Cohen and White did not solve the full partial integrodifferential equation numerically but rather they reduced it, using a perturbation series, to a system of ordi nary differential equations which they solved with an Adams-Bashfort-Moulton numerical time-stepping scheme [24]. 1.2 Thesi s Outl i ne In this work we solve the model derived by Weinstein [73] using a novel for mulation for pseudospectral differentiation matrices, the Eigen-decomposition Pseudo-Spectral (EPS) method. The governing equation is a nonlinear Volterra Partial Integrodifferential Equation (VPIDE) of the second kind. These equa tions are particulary difficult to solve not only because they are nonlinear but also because the integral term poses some numerical challenges. Since this t erm is cumulative, over time, it collects round-off error and can pose stability prob lems for a numerical solver. Singh [64, 65] solved the problem using a finite- element method with a time-stepping technique derived by Patlashenko et. al. [53]. However, Singh's solution for the liquid volume fraction was unsatisfac tory in t hat it attained nonphysical values. We will show t hat one can obtain higher-accuracy t han the method suggested in [53] by using a method-of-lines 7

(MOL) approach. In this context we have extended the work of Weinstein and Singh as well as added mathematical rigor to their results. In Chapter 2 we provide a review of HMT relevant t o the derivation of the governing continuity equation and nondimensionalize the model arriving at two non-dimensional parameters, De (Deborah number [56]) and KS (ratio of viscoelasticity to diffusivity). In Chapter 3 we analyze the VPIDE by finding a sufficient solution space t hat provides a proof of existence and uniqueness on a cylindrical geometry, under a specific set of initial boundary conditions t hat match the physics of the drug-delivery application. We reduce a special case of the VPIDE, using a derived similarity variable, to an ordinary differen tial equation and solve the resulting boundary value problem using a shooting method. We find an analytic solution to the Flory-Huggins model, t hat was re-derived in [73], which is another special case of the VPIDE. These solutions provide us with an expectation for the behavior of the solutions to the model. In Chapter 4 we provide an overview of pseudospectral methods and introduce the EPS method. We use the same (spatial) regularity conditions imposed in the existence-uniqueness proof to derive an error formula for the EPS method. The numerical solver for the VPIDE is provided in Chapter 5. We use the solver on the Flory-Huggins model with bot h the EPS method and a conventional pseu dospectral method (Chebyshev collocation) comparing the numerical solution to the analytic solution from Chapter 4. This example allows us to compare the EPS method against Chebyshev collocation as well. In Chapter 6 we conduct several numerical experiments t o study the behavior of solutions to the VPIDE under a variety of conditions. We test different diffusion coefficients and per- 8

meability models as well as test the model's sensitivity to the non-dimensional parameters derived in Chapter 2. Chapter 7 contains ideas for further research. 9

2. Fundament al Conti nui ty Equat i on In this chapter we review some of the results from [73] relevant to the deriva tion of the transport equations being solved. This work relies heavily on Hybrid Mixture Theory (HMT) for which there are several excellent resources available containing a rigorous and detailed overview [8, 9, 73]. As such, we will not reca pitulate HMT but will give a brief overview of the relevant aspects and reference the listed resources as required. We also offer a physical interpretation for the derived equations couched in terms of a drug-delivery model. 2.1 Bri ef Overvi ew of Hybri d Mi xt ure Theory HMT involves upscaling field equations which govern the motion of materials from the microscale to the macroscale via volume averaging and then obtaining restrictions on the form of the constitutive equations by using the second law of thermodynamics, also formulated as a field equation, at the macroscale. The field equations often include the conservation of mass, conservation of linear and angular momentum balance, and conservation of energy. Recall t hat constitutive equations are specific t o the material being modeled - examples include Fourier's Law and Darcy's Law. There are several other upscaling techniques such as homogenization, however, since we base our analysis on equations developed by Weinstein [73] we focus on HMT. The averaging procedure used in HMT is for any multi-scale, multi-constituent, multi-phase material. We consider a two-scale (microscale and macroscale) two phase material - liquid (I) and solid (s), each composed of one constituent within 10