www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE A Particular Solution for a Two-Phase Model with a Sharp Interface David A. Ekrut, Nicholas G. Cogan Department of Biological Mathematics Florida State University Tallahassee, Florida, United States ekrut@math.fsu.edu, cogan@math.fsu.edu Received: 23 October 2014, accepted: 8 March 2015, published: 28 April 2015 Abstract—Two-phase models can be used to de- scribe the dynamics of mixed materials and can be applied to many physical and biological phenomena. For example, these types of models have been used to describe the dynamics of cancer, biofilms, cytoplasm, and hydrogels. Frequently the physical domain separates into a region of mixed material immersed in a region of pure fluid solvent. Previous works have found a perturbation solution to capture the front velocity at the initial time of contact between the polymer network and pure solvent, then approximated the solution to the sharp-interface at other points in time. The primary purpose of this work is to use a symmetry transformation to capture an exact solution to this two-phase problem with a sharp-interface. This solution is useful for a variety of reasons. First, the exact solution replicates the numeric results, but it also captures the dynamics of the volume profile at the boundary between phases for arbitrary time scales. Also, the solution accounts for dispersion of the network further away from the boundary. Further, our findings suggest that an infinite number of exact solutions of various classes exist for the two-phase system, which may give further insights into the behaviors of the general two-phase model. Keywords-Multi-phase modeling; Two-phase mod- elling; Free boundary problems; Gel Dynamics; Analytic solutions; Exact solutions. I. INTRODUCTION Two-phase models are useful for capturing the interactions between fluids and/or viscoelastic ma- terial. Each phase is averaged over a control volume, where the volume-averaged phases are incompressible. There is no inertial component to the system, and the phases are immiscible. Each phase is governed by conservation equations. These models have been successful at describing how emergent structures develop though the inter- actions of the two phases. There are several known applications. Breward et. al. [1] developed a two-phase model to understand the role of viscosity and drag- friction in avascular tumor growth. An asymptotic solution solved explicitly for the volume fraction revealed that in the absence of viscosity and friction, tumor growth was regulated by oxygen tension. Numerical simulations showed increases in either the drag coefficient or viscosity param- eter reduces the speed tumor growth. This leads credence to the notion that the invasiveness of tumor cells is related to the viscosity of the cells. Well-differentiated cells are known to grow more slowly and considered more viscous due to over- lapping filopodia. Whereas, poorly differentiated Citation: David A. Ekrutl, Nicholas G. Cogan, A Particular Solution for a Two-Phase Model with a Sharp Interface, Biomath 1 (2015), 1503081, http://dx.doi.org/10.11145/j.biomath.2015.03.081 Page 1 of 14 http://www.biomathforum.org/biomath/index.php/biomath http://dx.doi.org/10.11145/j.biomath.2015.03.081 D. A. Ekrutl, A Particular Solution for a Two-Phase Model ... (less viscous) cells repel one another, contributing to the spread of tumors. An extension of this model with an additional phase [2] contrasts the role of the expansive growth (passive response) and foreign body hypotheses (active response) in tumorigenesis. Numerical simulations showed capsule formation could not result from an active response. Another model [3] was used to describe avascular tumor as a two-phase system where tumor spheroids exist in two states, one solid and one liquid. Time independent solutions reveal tumor size increases at an optimal rate of cell proliferation under nutrient-rich stress-free condi- tions. Simulations also provided a critical region for which a necrotic core forms at the tumor’s center. Several forces are required to balance conser- vation of momentum. For two-phase models, the viscosity of the phases and interstitial friction must be accounted, but for biofilm morphology, in ad- dition to hydrostatic pressure, osmotic pressure is also needed. One such model [4] describes the role of a network comprised of an extra-cellular poly- meric substance (EPS) in structural development in biofilm. Numerical simulations indicate as EPS is produced by bacteria, a rise in osmotic pressure contributes to the expansion of the biofilm region. Two-phase dynamics have also been used to sim- ulate biofilm growth and cell motility [5], [6]. A mobile cell contains polymer network phase com- prised of actin filaments, intermediate filaments, and microtubules. This phase is the exoskeleton to a cytoplasmic phase. The network contracts to propel the cell forward. Numerical simulations of these models have shown to contain traveling wave solutions. Another biological model describes to formation of channels in biofilm [7]. Steady-state analysis suggests that there is an optimal range for the pressure gradient to drive the formation of a channel between two flat plates. When regions occupied by differing materials have free boundaries, numerical methods are use- ful to track the sharp interface. The location of the interface can be followed explicitly by inter- face tracking methods [8]. Alternatively, interface capturing can be used to implicitly solve the same equations throughout the domain by capturing the appropriate interface conditions [9]. One such interface capturing method given by Du et. al. [10] has analyzed the behavior of a free boundary problem of a two-phase viscous fluid mixture with a prevalent viscosity in a single phase. The solution found by Du et. al. is pertur- bation solution of the front velocity at time t = 0 for a vanishing solvent phase. This solution was built to explore how the velocity of the interface moves in a consistent manner to develop numerical methods to handle the free boundary problem. The velocity is then tracked numerically for various initial profiles with the interface capturing method developed by the group. In each instance, the numerical solution is compared to the asymptotic solution and found to be accurate. In part, the purpose of this paper is to explore the accuracy of the perturbation solution given by [10] in comparison to an exact solution, which was found using symmetry analysis, also called Lie’s classical method. In each model previously dis- cussed, numerical, perturbation, and semi-analytic methods were used to provide insights into the behaviors of interest. And though these methods have had some successes in assessing two-phase models, few attempts have been made to attain generalized behavior of these systems with exact solutions. Lie’s method produces symmetry transforma- tions which can reduce a system of Partial Dif- ferential Equations (PDEs) in one spacial dimen- sion to a system of ODEs. These symmetries are generated by introducing infinitesimal transforma- tions, which leave the original system invariant. For classical symmetry analysis, expansion of this infinitesimal transformation, produces a linear sys- tem of PDEs, called determining equations, whose solutions provide the forms for the symmetry transformations. Non-classical methods have also been developed which, in some cases, lead to additional symmetries. The infinitesimal transfor- mations give rise to highly non-linear determining equations and can be difficult or impossible to Biomath 1 (2015), 1503081, http://dx.doi.org/10.11145/j.biomath.2015.03.081 Page 2 of 14 http://dx.doi.org/10.11145/j.biomath.2015.03.081 D. A. Ekrutl, A Particular Solution for a Two-Phase Model ... solve. For this reason, the analysis in the paper only includes the classical method, as it recovers the solution given by [10] that we are seeking. Lie’s classical method for producing symmetries has been successful in generating exact solutions for a system of PDEs describing viscous flow through expanding channels [11]. In this work conservation laws and point symmetries provide reductions, some of which lead to exact solutions of the flow in deformable channels. For elliptic, hyperbolic and mixed-type PDEs for Ricci flow, Wang [12] found several solutions, including trav- eling wave solutions, to hyperbolic geometric flow of Riemann surfaces. The work by Cimpoiasu et. al [13] used Lie symmetries to produce classes of solutions for the 2D nonlinear heat equation. It has also been shown that Lie symmetries generate the similarity solution for a class of (2+1) nonlinear wave equation [14]. In this paper, we generate an exact solution for the two-phase model using a point symmetry. In the first section, we outline a derivation of a two-phase system that represents the simplest version of the model and can be adapted for a variety of physical situations. Next, we briefly discuss how to develop symmetries and find a time translation, scaling symmetry, and a general Galilei time group. In the third section, we use a symmetry transformation to reduce the system of PDEs to an invariant system of ODEs. We make parameter assumptions similar to Du et. al. [10] to recover the exact velocity for their asymptotic solution and compare the exact to the perturbation solution. It is shown that the approximated free boundary solution is a close approximation to the general solution for t = 0. In the fourth section, we vary which physical driving forces dominate the two-phase model and generate additional exact solutions to the system. In the final section of this work, we discuss potential uses of exact solutions for the two-phase model and future directions of this work. II. THE TWO-PHASE MODEL In this section, we derive the equations to de- scribe a two-phase model as seen in the kinetics of biological gels as described in [5], [10]. Gels swell and deswell due to ionic fluctuations and chemical triggers. An example of this occurs in crawling cells. Myosin converts chemical energy in the form of ATP into mechanical energy by causing actin filament to contract, propelling cells into motion. Neutrophils and macrophages, cells integral to the immune system of humans, respond in this manner. Chemical gradients are left by cells foreign to the immune system, leaving a chemotactic trail for the immunological cells to follow [15]. Like in [10], we assume the viscous terms are prominent forces and inertial terms are negligi- ble. Gels are composed of a polymeric network given by φ1 and a fluid solvent φ2. Both phases are treated as Newtonian fluids that are immisci- ble. When considering the redistribution of mass within a control volume, the flux of the network is given by ∇· (φ1u1), where the network moves with a velocity u1. A similar argument is made for the solvent to give the following equations to conserve mass. ∂ ∂t (φ1) + ∂ ∂x (uφ1) = 0, (1) ∂ ∂t (φ2) + ∂ ∂x (vφ2) = 0, (2) where the sum of the volumes saturate to a fixed control volume, φ1 + φ2 = 1. Several forces act upon the network. The first is the force due to the network stress tensor σ1, which includes the viscous stress tensor and mass production. σ1 = µ̂1(∇u1 + ∇uT1 ) + λ1∇·u1, where µ̂1 is the shear viscosity and λ1 is the bulk viscosity. In 1-D, this becomes σ1 = µ1 d dx u1, (3) where µ1 = 2µ̂1 + λ1. Another force that we include is the frictional force created by interstitial interactions between Biomath 1 (2015), 1503081, http://dx.doi.org/10.11145/j.biomath.2015.03.081 Page 3 of 14 http://dx.doi.org/10.11145/j.biomath.2015.03.081 D. A. Ekrutl, A Particular Solution for a Two-Phase Model ... phases. If both fluids move in unison or if either volume fraction becomes negligible, drag will van- ish. With a frictional coefficient given by ξ, this drag force is given by ξφ1φ2(φ1 −φ2). Next, we need to account for both the hydrostatic pressure and osmotic pressure caused by swelling. If P is the total hydrostatic pressure, then the total pressure P acting on the network is given by φ1∇P . Ionizing chemicals in the solvent can cause the gel to absorb or release the fluid solvent, causing an osmotic pressure gradient ∇ψ(φ1) acting on the network. For this reason, φ1 is considered the active phase. For the form of the osmotic pressure term, we follow Cogan et. al. [5] and the references therein, and assume that ψ(φ1) = k2φ21(φ1 −φ0). The constant k2 accounts for the effects of the ionic environment, polymeric structure, and sol- vent concentration that contribute to swelling and deswelling. The value of φ0 is a reference vol- ume fraction. This structure allows for osmotic pressure to vanish in the event of φ1 = 0 or at some reference fraction φ0 that can be determined experimentally for various physical applications. Assuming constant shear and bulk viscosity, the momentum of these moving fluids can be given by balancing the forces described above. µ1 ∂ ∂x ( φ1 ∂ ∂x u ) + φ1 ∂ ∂x P(φ1,φ2) (4) − ∂ ∂x ψ(φ1) − ξφ1φ2(u−v) = 0 Similar arguments can be made to derive the forces of momentum within the solvent. The solvent is a Newtonian fluid with only viscous stresses acting on it. Fluid pressure acts on the solvent, but osmosis does not create pressure on the fluid itself. The fluid is actively absorbed and released by the gel. The final force is the drag or frictional force created by interstitial interactions. Combining these gives the momentum for the solvent. µ2 ∂ ∂x ( φ2 ∂ ∂x v ) + φ2 ∂ ∂x P(φ1,φ2) (5) + ξφ1φ2(u−v) = 0, where µ2 is the viscosity of the solvent. Summing (4) and (5) gives the following equation. µ1 ∂ ∂x ( φ1 ∂ ∂x u ) + µ2 ∂ ∂x ( φ2 ∂ ∂x v ) + (φ1 + φ2) ∂ ∂x P(φ1,φ2) − ∂ ∂x ψ(φ1) = 0. Since φ1 + φ2 = 1, this becomes µ1 ∂ ∂x ( φ1 ∂ ∂x u ) + µ2 ∂ ∂x ( φ2 ∂ ∂x v ) (6) + Px − ∂ ∂x ψ(φ1) = 0, where Px = ∂∂xP(φ1,φ2). Solving for Px gives Px = ∂ ∂x ψ(φ1) −µ1 ∂ ∂x ( φ1 ∂ ∂x u ) (7) − µ2 ∂ ∂x ( φ2 ∂ ∂x v ) . Next, we substitute φ2 = 1−φ1 in the equations of mass (1) and (2), and the momentum equation (4) to find the following system for analysis. ∂ ∂t (φ1) + ∂ ∂x (uφ1) = 0, (8) − ∂ ∂t (φ1) + ∂ ∂x (v(1 −φ1)) = 0, (9) µ1 ∂ ∂x ( φ1 ∂ ∂x u ) − ∂ ∂x ψ(φ1) (10) +φ1Px − ξφ1(1 −φ1)(u−v) = 0. Together equations (8-10) can be reduced to a system of ODEs using the following transforma- tion. u = f(t−αx), v = g(t−αx), (11) φ1 = m(t−αx), Biomath 1 (2015), 1503081, http://dx.doi.org/10.11145/j.biomath.2015.03.081 Page 4 of 14 http://dx.doi.org/10.11145/j.biomath.2015.03.081 D. A. Ekrutl, A Particular Solution for a Two-Phase Model ... where f, g, and m are to be determined and α is an arbitrary constant describing wave speed. Traveling wave solutions have been shown to exist for the two phase system [6]. For this reason, if one were to guess an invariant transforma- tion to reduce this system, the general traveling wave solution (11) may seem like an obvious first choice. But, this specific transformation came from a more general transformation found using symmetry analysis. Before producing the general transformation, a brief explanation of symmetry analysis is given in the following section. III. SYMMETRY ANALYSIS In this section, we give a brief explanation of the method for generating the invariant transforma- tions that will be used to generate exact solutions in later sections. For systems of PDEs in 1-D, symmetry transformations reduce the PDEs to a system of ODEs. Derived by Sophus Lie [16], Symmetry Analysis is the mathematical method for finding transformations to a system of PDEs that leaves the set of equations invariant, or un- changed. More recently, there has been substantial literature regarding symmetry methods. For further details, we refer the reader to books by Hydon [17], Bluman and Kumei [18], and Olver [19]. The following coordinate change is called the infinitesimal transformations. These can be thought of as a local perturbation on the original coordinate system. φ̄1 = φ1 + Φ1(t,x,u,v)� + O(� 2), t̄ = t + T(t,x,u,v)� + O(�2), x̄ = x + X(t,x,u,v)� + O(�2), ū = u + U(t,x,u,v)� + O(�2), (12) v̄ = v + V (t,x,u,v)� + O(�2), where Φ1, T , X, U, and V are called the infinites- imals. In general, one seeks to find invariance of a system of differential equations of the form Fi(t,x,u,v,φ1,ut,vt,φ1t,ux,vx,φ1x, ...) = 0, (13) with i = 1, 2, . . . ,n, where u, v, φ1 are functions of t, x. In the specific case of our two-phase model, the system Fi is given by the equations (8- 10). Under (12), a set of differential equations is produced for the infinitesimals T , X, U, and V . These differential equations are called the deter- mining equations because they determine the form for the infinitesimals. Solving these determining equations produced by (12) provides invariant transformations for the differential equations given by (13). The following is called the invariant surface condition, so called because it leaves the solution surface invariant under the change of coordinates. Tut + Xux = U, (14) Tvt + Xvx = V, (15) Tφ1t + Xφ1x = Φ1. (16) When the infinitesimals are solved in conjunc- tion with the invariant surface condition given by (14-16), the solutions u, v, and φ1 provide a trans- formation which reduces the original PDE (13) to an ODE. In other words, by using Lie’s method to find an infinitesimal change of coordinates, a two variable PDE can be reduced to an equation of a single variable to become an ordinary differential equation (ODE). Taking the physical nature of the problem into account, these reductions can lead to exact solutions to the PDE. Applying the transformation given by (12) on (8-10) yields a large system of linear PDEs. The determining equations are solved interac- tively to give the forms of the infinitesimals. Φ1 = 0, T = α, X = δx + Γ(t), (17) U = δu + d dt Γ(t), V = δv + d dt Γ(t). Due to the size of the equations, details of the determining equations are omitted. For more Biomath 1 (2015), 1503081, http://dx.doi.org/10.11145/j.biomath.2015.03.081 Page 5 of 14 http://dx.doi.org/10.11145/j.biomath.2015.03.081 D. A. Ekrutl, A Particular Solution for a Two-Phase Model ... details on an example, see the details given in the Appendix (A). In order for the PDEs given by the determining equations to be satisfied, two cases arise. Either δ = 0 or δ 6= 0. If δ 6= 0, then the friction coefficient ξ given in the momentum equations vanishes. The transformation given by α is a time translation, δ is a scaling symmetry, and Γ(t) is a general time dependant Galilei group, as used in fluid mechanics [20]. These symmetries can be used to find invariant reductions in the original system. Notice that for δ = 0 and Γ(t) = 1 in (17) and solving for u, v, and φ1 in (14-16) gives the transformation u = 1 α + f̂(t−αx), v = 1 α + ĝ(t−αx), φ1 = m(t−αx), Letting f̂ = f(t−αx)− 1 α and ĝ = g(t−αx)− 1 α gives the transformation (11). It should be noted that for our purposes, we are only interested in pursuing a classical symmetry analysis to recover the solution presented by Du et. al. [10]. It is possible that more solutions will arise from other methods as well. Non-classical symmetries arise in many cases. In the work performed by Arrigo et. al. [21], a nonclassical symmetry is emitted by a class of Burgers’ system. The Steinbergs symmetry method has provided exact solutions and reductions to the Calogero- Bogoyavlenskii-Schiff equation [22]. The Gardner method can generate an infinite hierarchy of sym- metries, as was shown with the KdV equations, Camassa-Holm, and sine-Gordon equations [23]. Non-classical symmetries have also been gener- ated for the fourth-order thin film equation using non-classical methods [24]. Further analysis could include any of these methods, as well as a classification of parameters which has the potential to produce more symme- tries. The purpose of this work is not an exhaustive search for symmetries, but an introduction to using symmetry methods to recover a more general solution to the two-phase problem described above and partially recovered by Du et. al. [10]. IV. RECOVERING THE EXACT SOLUTION FOR A FREE BOUNDARY PROBLEM As discussed in [10], since the viscosity of the solvent is of a much higher magnitude than that of the fluid, we assume the solvent viscosity µ2 is zero. Since, φ1 + φ2 = 1, we have φ2 = 1 −φ1. Now, we replace φ2 in the equations of momentum (4-5) and find µ1 ∂ ∂x ( φ1 ∂ ∂x u ) + φ1Px − ∂ ∂x ψ(φ1) (18) −ξφ1(1 −φ1)(u−v) = 0, (1 −φ1)Px + ξφ1(1 −φ1)(u−v) = 0, (19) where u, v, and φ1 are all functions of t, x as previously discussed and Px is the pressure gradient. Next, we solve (19) for Px to find Px = −ξφ1(u−v). (20) We see the mass equations (1-2) have now become ∂ ∂t (φ1) + ∂ ∂x (uφ1) = 0, − ∂ ∂t (φ1) + ∂ ∂x (v(1 −φ1)) = 0, Summing these two equations of mass gives ∂ ∂x (uφ1 + v(1 −φ1)) = 0. Imposing the average velocity is zero, we have uφ1 + v(1 −φ1) = 0, which gives v = − φ1 1 −φ1 u. (21) To match the form of equations given by Du et. al. [10], we let the osmotic swelling term take the form ψ(φ1) = φ1Ψ(φ1). This, together with Biomath 1 (2015), 1503081, http://dx.doi.org/10.11145/j.biomath.2015.03.081 Page 6 of 14 http://dx.doi.org/10.11145/j.biomath.2015.03.081 D. A. Ekrutl, A Particular Solution for a Two-Phase Model ... Multiphase φ1,φ2>0 φ1=0 Γ Ω1 Ω2 Fluid Fig. 1. This shows the region Ω2 of pure solvent (φ1 = 0) separated at the boundary Γ from the region Ω1 containing the mixture of both phases. (20-21), reduces the equation (18) to the following equation. (µ1φ1(u)x)x−(φ1Ψ(φ1))x− ξφ1 1 −φ1 u = 0. (22) As in Figure (1), we assume the mixture occu- pies the interior region (Ω1), while pure solvent occupies the external region (Ω2). As the gel- mixture swells/deswells, the interface between the regions (Γ) moves. To specify the motion Du et. al. impose standard jump conditions: [µ1φ1(u)x −φ1Ψ(φ1)] = 0 [u] = 0. The solution found by Du et. al. [10] approx- imates the front velocity for the free boundary problem at time t = 0. The solution for a piece- wise constant profile is given by, φ1 = { φ− if x < 0 φ+ if x > 0 , and the following can be derived u = { Ceβ−x if x < 0 Ce−β+x if x > 0 , (23) where β± = √ ξ µ1(1 −φ1±) , and C = −φ+Ψ(φ+) + φ−Ψ(φ−) µ1(φ+β+ + φ−β−) . The solution (23) was derived by assuming φ1+ → 0 at t = 0. In biological gels, regions of gel separate from regions of pure solvent. So, it is reasonable to assume that the network phase vanishes in this region of pure solvent. To make a graph of the solution given by (23), we assign the following initial profile. φ1 = { φ− = 1 6 if x < 0 φ+ = 0 if x > 0 . (24) The parameters used to generate the graphs are taken from [10], but are repeated in (I) for convenience. The graph Figure (2) represents the velocity front for a swelling gel in contact with a fluid solvent. This perturbation solution is an approximation for the velocity front at t = 0. However, there exists an exact solution to this system that captures this behavior for all values of φ1 at any point in time. For the infinitesimals given by (17), let δ = 0 and Γ(t) = 1. Solving the invariant surface condi- tion for u, v, and φ1 will lead to (11) in terms of the variable r = t − αx. As with the case found with solving for (23) , we assume the viscosity of the second phase is negligible in comparison to that of the first phase, letting µ2 = 0. To make the analysis easier, we allow only for swelling in the active phase, making φ0 = 0. Applying (11) reduces (8-10) to a single ODE. µ1α 4(αf − 1)2f ′2 −µ1α3(αf − 1)3f ′′ (25) +3k2γ 2α3f ′ + ξ(αf − 1)4 = 0, where g = 1 α , (26) m = γ αf − 1 , (27) Biomath 1 (2015), 1503081, http://dx.doi.org/10.11145/j.biomath.2015.03.081 Page 7 of 14 http://dx.doi.org/10.11145/j.biomath.2015.03.081 D. A. Ekrutl, A Particular Solution for a Two-Phase Model ... Fig. 2. This is the perturbation solution given by (23) at time t = 0. This shows the velocity for a region of vanishing network (top) in the region x > 0 for an initial profile (24) (bottom). This would represent expectations of a velocity front for a swelling gel to contact a fluid solvent. It should be noted that we could have just as easily solved (8-10) for f and left m to be determined. We are choosing to leave f general to assess the behavior of the velocity, because we wish to show the exact solution approximated by Du et. al. [10] can be recovered. Multiple solutions exist for (25). First, we at- tempt to recover an exponential solution similar in form to (23). If we assume the viscosity of the network phase φ1 has a greater impact on the system than viscosity and interstitial friction, as was assumed by Du et al [10]. we can divide by µ1. This gives the following equation from (25). α4(αf − 1)2f ′2 −α3(αf − 1)3f ′′ (28) +3 k2 µ1 γ2α3f ′ + ξ µ1 (αf − 1)4 = 0. For µ1 of a much larger magnitude than k2 and ξ, this becomes the following: α(αf − 1)2f ′2 − (αf − 1)3f ′′ = 0, (29) whose solution is f(r) = eα(κr+λ) α + 1 α . (30) This makes the analytic solution for the original system (1-2) and (4-5) to be φ1 = γ̂ αf − 1 = γe−α(κ(t−αx)+λ), φ2 = 1 −γe−α(κ(t−αx)+λ), u = eα(κ(t−αx)+λ) α + 1 α , v = 1 α , with µ1 = 0, k2 = 0, ξ = 0. The parameters of this solution can be matched to the parameters of the solution given by (23). We can see that if k = − β α2 and λ = 1 α ln(αC − 1), then the solution found above becomes: φ1 = γ C e β α t−βx, φ2 = 1 − γ C e β α t−βx, u = Ce− β α t+βx + 1 α , (31) v = 1 α , The parameter α remaining in the velocity of (31) gives flexibility on scaling time and adjusting the orientation of the velocity. Notice, as α →∞, this solution is the same as (23). The velocity be- comes identical, and the volume fraction becomes constant, as in the perturbation solution provided by (23). So, in essence, we have recovered the time function that was missed by the perturbation Biomath 1 (2015), 1503081, http://dx.doi.org/10.11145/j.biomath.2015.03.081 Page 8 of 14 http://dx.doi.org/10.11145/j.biomath.2015.03.081 D. A. Ekrutl, A Particular Solution for a Two-Phase Model ... method used to find (23). Next, we match the numerical results of (23) for the parameters given by (I). To do this, we set t = 0 and separate the solution for the velocity as follows. u = { Ceβx + 1 α , if x ≤ 0 Ce−βx + 1 α if x > 0 , (32) In Figure (3), we can see the solutions given by (23) and (32) super-imposed on the same graph for parameter values given by (I). It is clear that the perturbation solution is a close approximation for the exact solution for large values of α. As we would expect from inspection of the solution (32), smaller values of α will adjust the exact solution away from the perturbation solution. The largest impact α has on the system is in regards to the time scale and solvent velocity. Large values of α require larger time steps for movement in the system, while decreasing the solvent velocity. β µ1 ξ α 1 0.0108 0.018 1000 10 0.0037 0.616 10000 100 0.000338 5.64 100000 TABLE I THE PARAMETERS GIVEN IN THE ROW BEGINNING WITH β = 1 GENERATES THE RESULTS IN (3) TOP. THE NEXT ROW FOR β = 10 GIVES (3) MIDDLE WITH THE FINAL ROW GENERATING (3) BOTTOM WITH β = 100. There are several benefits of finding the ex- act solution, instead of using numerical methods. First, numerical results have a difficult time captur- ing the behavior at the region of contact between the phases, while the analytical solution easily gives interface behaviour without computationally expensive coding, as can be seen in 4. Here we can see the region of network at t = 0 moving uniformly away from the initial contact region x = 0. Smaller values of β fail to capture the sharp interface. But as β increases to β = 100, we see the interface remains sharp as time increases. This is expected, as these results coincide with the numerical simulations found in [10] by a moving mesh. Fig. 3. These are the perturbation solutions given by (23) at time t = 0 graphed with the solution given by 32 with β = 1 and the corresponding values for µ1 and ξ described in (I) given by the top, β = 10 middle, and β = 100 on bottom. The perturbation solutions are a close approximations for the exact solutions near the region of separation. We can see that the shape of each solution is preserved for each set of parameters, though the scale is modified. Biomath 1 (2015), 1503081, http://dx.doi.org/10.11145/j.biomath.2015.03.081 Page 9 of 14 http://dx.doi.org/10.11145/j.biomath.2015.03.081 D. A. Ekrutl, A Particular Solution for a Two-Phase Model ... Fig. 4. For network (φ1) profiles for t = 0 (top) and t = 1000 (bottom). As β increases, the exact solution becomes more accurate at capturing the expected behavior at the sharp interface. As we have seen, the analytic solution recovers the perturbation solution as well as the numerical results given by Du et. al. [10]. However, this is just a single solution to the nonlinear equation given by (25). It is possible that the other solu- tions are extraneous, but more likely, additional solutions describe other physical or biological phe- nomenon yet to be determined. Further exploration will be required to fit these solutions, but we look at the others here. A. Other Solutions to (25) Being a non-linear system, the solution to (25) is not unique. Even though the transformation given by (11) will clearly give traveling wave solutions, the structure of the traveling wave for each solution can vary widely as can be seen with the next two examples. If the viscosity of the network is negligible µ1 = 0, the following solution to (25) is given by f(r) = − k 1 3 2 γ 2 3 ξ 1 3 α 1 3 (r + δ) 1 3 + 1 α , (33) where r = t−αx. The structure of this solution is different from (30) in several ways. When plotting at a single moment in time t = 0, it looks like a pulse as seen by the first curve in figure (5). When animated (30) can be seen as a traveling wave solution, given by the black curves which moves in the positive t direction. Fig. 5. The solution given by (33) plotted at t = (0, 10). The first curve is at t = 0. As seen by the black curves, the velocity front travels like a wave as time increases. Alternatively, if the osmotic pressure has less of an impact than viscosity and friction, then with k2 = 0 as seen in the absence of ionizing agents for gels, we find the following solution f(r) = e α(−κr+λ)+ ξ 2µ1 r2 α + 1 α . (34) Like (30) this solution is exponential, but as seen in (6) the quadratic term gives an unbounded traveling wave. The velocity at t = 0 is given by the first curve. As time increases, the front velocity Biomath 1 (2015), 1503081, http://dx.doi.org/10.11145/j.biomath.2015.03.081 Page 10 of 14 http://dx.doi.org/10.11145/j.biomath.2015.03.081 D. A. Ekrutl, A Particular Solution for a Two-Phase Model ... travels as a wave moving to the right. This does not seem to have any physical analogue since the velocities are unbounded. Fig. 6. The solution given by (34) plotted at t = [0, 10]. The first curve is at t = 0. The velocity front shifts to the right as time increases, moving as a traveling wave. V. ADDITIONAL SOLUTIONS TO THE TWO-PHASE MODEL This section also provides theoretical solutions, which may or may not have physical relevance. We explore them here to account for the multitude of solutions that are emitted by the system (8-10). There are other two-phase models from physics that might have solutions contained here. For ex- ample, one such model describes granular flows where air is considered a non-viscous (nondense) phase with the rocks, debris, and other materials considered as a second highly viscous (dense) phase [25], [26]. Within these works, numerical simulations describe the flow behaviors air has on granular flow. The results suggest that drag has more than a negligible effect on the flow of granular materials of finite mass. The traveling wave solutions provided by (11) are given by a simple choice for Γ(t) in (17). Here, we explore different choices for the transformation and follow the reduction of the PDEs to ODEs. Then, we derive solutions to the ODEs by consid- ering various changes in the physical nature of the problem. By adjusting which physical parameters are the dominating driving force in the problem, we can generate different solutions, which may prove useful in exploring the nature of physical and biological phenomenon. First, we let δ = 0 and Γ = 1 t in (17). Solved with (14-16) will give the following transforma- tion. u = 1 αt + f ( x− ln(t) α ) , v = 1 αt + g ( x− ln(t) α ) , φ1 = m ( x− ln(t) α ) . Neglecting solvent viscosity, µ2 = 0, reduces the original system (8-10) to the following ODE, µ1f 2(α−f)f ′2 −µ1f3f ′′ (35) −3k2α(2φ0f − 3α)(α−f)f ′ − ξf5 = 0. with m = α f , g = αf α−f . Again, if we assume the dominating force is the viscosity and set k2 and ξ2 to zero, this can be solved to give f = κeλ(x− ln(t) α ) = κ α √ tλ eλx. (36) The complete solution to (1-2) and (4-5) becomes u = 1 αt + κ α √ tλ eλx, v = 1 αt + α α− κα√ tλ eλx κ α √ tλ eλx, φ1 = κ α √ tλ eλx, φ2 = 1 − κ α √ tλ eλx. If we assume friction and pressure dominate and let µ1 = 0, the ODE yields no real solution without further assumptions on the constants of integration. Biomath 1 (2015), 1503081, http://dx.doi.org/10.11145/j.biomath.2015.03.081 Page 11 of 14 http://dx.doi.org/10.11145/j.biomath.2015.03.081 D. A. Ekrutl, A Particular Solution for a Two-Phase Model ... This may imply that viscosity is required for non- constant solutions. Another choice for (17) is to let δ remain arbitrary and to consider ξ = 0, a requirement for invariance to be satisfied. Again we consider cases where the viscosity of the solvent is negligible, µ2 = 0. Choosing Γ = 0 we find the following transformation u = f(xe− δ α t)e δ α t, v = g(xe− δ α t)e δ α t, φ1 = m(xe − δ α t), which reduces (8-10) to the following three ODEs. (αf −rδ)m′ + αmf ′ = 0, −(αg −rδ)m′ + α(1 −m)g′ = 0, (37) (1 −m)(µ1f ′ −k2m(2φ0 − 3m))m′ +µ1m(1 −m)f ′′ = 0, with r = xe− δ α t. If both osmotic pressure and viscosity are negligible such that k2 = 0 and µ1 = 0, then the following solution satisfies (37). f = γe−λr + δ λr − 1 λα , g = (λr − 1)δeλr λα(δeλr − 1) + κ λα(δeλr − 1) , m = δeλr. The complete solution to (1-2) and (4-5), is u = γe−λxe − δ α t + δ λxe− δ α t − 1 λα e δ α t, v = (λxe− δ α t − 1)δeλxe − δ α t λα(δeλxe − δ α t − 1) + κ λα(δeλxe − δ α t − 1) (38) φ1 = δe λxe− δ α t , (39) φ2 = 1 −φ1. It should be noted that if either viscosity is the dominating force with k2 = 0, or if osmotic pressure is the dominating force with µ2 = 0, then m, f, and g are constant. This implies that friction is required for non-constant solutions. This is different from before, where we found viscosity to be the driving force for the model. In summary, we have found that each solu- tion requires a dominating force to generate non- constant solutions. This gives flexibility in assess- ing the two-phase model and suggests that exact solutions may exist for many differing physical phenomenon of interest. For example, it is possible that the solution given by (39) can be matched to results consistent to granular flow, since friction as a necessary component for granular flow [25], [26]. VI. DISCUSSION In this work, we found an exact solution which accurately replicates the results from a previously found numerical results. It has been shown that for α → ∞, the analytic solution found here is exactly the perturbation solution found by Du et. al. [10]. The exact solution has the benefit of time dependence, which is useful for assessing behavior of the two-phase system without the im- plementation of numerical methods. Additionally, we showed that many traveling wave solutions arise from the two-phase problem. Due to the time dependent general Galilei group, we have an unlimited number of choices to adjust the speed of the wave through time. These solutions also require specific dominating forces to attain. It is possible that such solutions only arise in specific physical circumstances. Though some of these solutions may be extraneous, further investigation is warranted to determine their uses. Although asymptotic and numerical methods yield useful information concerning the behavior of multi-phase systems, these methods require substantial efforts. Exact solutions have the benefit of being computationally inexpensive to simulate, and with Lie symmetries, are relatively simple to generate. There are several directions for future analysis that arise from this work. First, exploring the be- havior of the additional solutions may give further Biomath 1 (2015), 1503081, http://dx.doi.org/10.11145/j.biomath.2015.03.081 Page 12 of 14 http://dx.doi.org/10.11145/j.biomath.2015.03.081 D. A. Ekrutl, A Particular Solution for a Two-Phase Model ... insights into the nature of dominating forces in the two phase system. This may give insight into specific physical phenomenon in which these addi- tional solutions may be esoterically relevant. Also, additional symmetries may exist, which could be found using non-classical methods. Solutions aris- ing from non-classical methods would then need to be assessed to determine relevant matching phys- ical or biological behavior. Additionally, biofilms typically include growth terms to account for the production of new network. It is possible that symmetry solutions can capture this behavior as well. APPENDIX Deriving the infinitesimals for the two-phase model generates a large system of linear PDEs. For this reason, we have provided details of the process by way of an example in this appendix. For further details, see [27]. Consider the following nonlinear first order PDE ut = u 2 x (40) Under the transformation t̄ = t + �T(t,x,u) + O(�2), x̄ = x + �X(t,x,u) + O(�2), ū = u + �U(t,x,u) + O(�2), to order �2 (40) becomes Ut + utUu −ut (Tt + utTu) −ux (Xt + utXu) − 2ux(Ux + uxUu −ut (Tx + uxTu) − ux (Xx + uxXu)) = 0. Using the original equation (40) to eliminate ut and grouping coefficients of ux, we have Ut + u 2 xUu −u 2 x ( Tt + u 2 xTu ) − ux ( Xt + u 2 xXu ) − 2ux(Ux + uxUu − u2x (Tx + uxTu) −ux (Xx + uxXu)) = Ut − (Xx + 2Ux) ux + (2Xx −Tt −Uu) u2x = (Xu + 2Tx) u 3 x + Tuu 4 x = 0. Invariance requires the coefficients of ux to be zero, providing us with the following system. U(t,x,u)t = 0, X(t,x,u)x + 2U(t,x,u)x = 0, 2X(t,x,u)x −T(t,x,u)t −U(t,x,u)u = 0, Xu + 2Tx = 0, Tu = 0. These are called the determining equations, be- cause they determine the forms of the infinites- imals. These are linear PDEs, which are easily solved with standard techniques of integration. So, we have the following form for the infinitesimals. T(t,x,u) = c1 + c2t + c3x + c4t 2 + c5tx + c6x 2, X(t,x,u) = c7 + c8 + c9x + c4tx + 1 2 k5x 2 − (2k3 + 2k5t− 4k6x) u, U(t,x,u) = k10 − 1 2 k8x− 1 4 k4x 2 + (2k9 −k2) u + k5xu− 4k6u2, where ci and ki are arbitrary constants of inte- gration. Together with the invariant surface con- dition given by Tut + Xux = U we can find a transformation to reduce (40) to an ODE. The form for the transformation will vary depending on choices for the constants ci and ki. REFERENCES [1] C. J. Breward, H. M. Byrne, C. E. Lewis The role of cell-cell interactions in a two-phase model for avascular tumor growth J. Math. Biol, 45 (2002) 125-152 [2] S. R. Lubkin, T Jackson Multiphase mechanics of capsule formation in tumors J. Biomech. Eng.124 (2002) 237-243 http://dx.doi.org/10.1115/1.1427925 [3] H. M. Byrne, L. Preziosi Modelling solid tumor growth using the theory of mixtures Math. Med. Biol. 20 (2003) 341-366 http://dx.doi.org/10.1093/imammb/20.4.341 [4] N. G. Cogan, J. P. Keener The role of the biofilm matrix in structural development Math. Med. Biol. 21 (2005) 147166 [5] N. G. Cogan, R. D. Guy Multiphase flow models of biogels from crawling cells to bacterial biofilms HSFP Journal4 (2009) 11-25 x Biomath 1 (2015), 1503081, http://dx.doi.org/10.11145/j.biomath.2015.03.081 Page 13 of 14 http://dx.doi.org/10.1115/1.1427925 http://dx.doi.org/10.1093/imammb/20.4.341 http://dx.doi.org/10.11145/j.biomath.2015.03.081 D. A. Ekrutl, A Particular Solution for a Two-Phase Model ... [6] J. M. Oliver, L. S. Kimpton, J. P. Whiteley, S. L. Waters, J. R. King Multiple travelling-wave solutions in a minimal model for cell motility Math. Med. Biol.30 (2013) 241272 [7] N. G. Cogan, J. P. Keener Channel formation in gels SIAM J. Appl. Math. 65 (2005) 1839-1854 http://dx.doi.org/10.1137/040605515 [8] R. D. Guy, A. L. Fogelson, J. P. Keener Fibrin Gel Formation in a Shear Flow Math. Med. and Bio. 24 (2007) 111-130 [9] W. J. Boettinger, J. A. Warren, C. Beckermann, A. Karma Phase-field simulation of solidification Annu. Rev. Mater. Res. 32 (2002) 163-194 http://dx.doi.org/10.1146/annurev.matsci.32.101901. 155803 [10] J. Du, R. D. Guy, A. L. Fogelson, G. B. Wright, J. P. Keener An interface-capturing regularization method for solving the equations for two-fluid mixtures Commun. Comput. Phys. 14 ( 2013) 1-25 [11] S. Asghar, M. Mushtaq, A. H. Kara Exact solutions using symmetry methods and conservation laws for the viscous flow through expanding-contracting channels App. Math. Model.32 (2008) 29362940 http://dx.doi.org/10.1016/j.apm.2007.10.006 [12] J. Wang Symmetries and solutions to geometrical flows Sci China Math 56 (2013) 16891704 http://dx.doi.org/10.1007/s11425-013-4635-8 [13] R. Cimpoiasu, R. Constantinescu Lie symmetries and invariants for a 2D nonlinear heat equation Nonlinear Analysis 68 (2007) 22612268 http://dx.doi.org/10.1016/j.na.2007.01.053 [14] M. Nadjafikhah, R. Bakhshandeh-Chamazkoti, A. Mahdipour-Shirayeh A symmetry classification for a class of (2+1)-nonlinear wave equation Nonlinear Analysis 71 (2009) 5164-5169 http://dx.doi.org/10.1016/j.na.2009.03.087 [15] S. Nagy, B. L. Ricca, M. F. Norstram, D. S. Courson, C. M. Brawley, P. A. Smithback , R. S. Rock A myosin motor that selects bundled actin for motility Proc. Natl. Acad. Sci.105 (2008) 96169620 http://dx.doi.org/10.1073/pnas.0802592105 [16] S. Lie Klassifikation und Integration von gewohnlichen Differentialgleichen zwischen x, y die eine Gruppe von Transformationen gestatten Math. Ann. 32 (1888) 213281http://dx.doi.org/10.1007/BF01444068 [17] [10.1017/CBO9780511623967] P. E. Hydon Symmetry Methods for Differential Equations: A Beginner’s Guide 6th edition, New York : Cambridge University Press (2000) [18] [10.1007/978-1-4757-4307-4] G. Bluman, S. Kumei Symmetries and Differential Equations Springer-Verlag, New York Inc. (1989) [19] [10.1007/978-1-4684-0274-2] P. J. Olver Applications of Lie Groups to Differential Equations Springer-Verlag, New York Inc. (1993) [20] T. Kambe Gauge principle and variational formulation for flows of an ideal fluid Acta Mechanica Sinica19 (2003) 437-452 [21] D. Arrigo, D. Ekrut, J. Fliss, L. Le Nonclassical symmetries of a class of Burgers systems J. Math. Anal. Appl. 371 (2010) 813-820 http://dx.doi.org/10.1016/j.jmaa.2010.06.026 [22] G. M. Moatimid, R. M. El-Shiekh, A. A. A. H. Al- Nowehy Exact solutions for CalogeroBogoyavlenskiiS- chiff equation using symmetry method Appl. Math. and Comput.220 (2013) 455-462 http://dx.doi.org/10.1016/j.amc.2013.06.034 [23] M. K. Srivastava, Y. Wang, X. G. Zhang, D.M.C. Nicholson, H. P. Cheng The Gardner method for symmetries Physical Review (2013) [24] K. Charalambous, C. Sophocleous Symmetry properties for a generalised thin film equation J Eng Math 46 (2013) 1-16 [25] L. T. Sheng, Y. C. Tai, C. Y. Kuo, S. S. Hsiau A two-phase model for dry density-varying granular flows Advanced Powder Technology 24 (2013) 132-142 http://dx.doi.org/10.1016/j.apt.2012.04.001 [26] S. P. Pudasaini A general two-phase debris flow model J. Geophys. Res. 117 (2011) 1-28 [27] D. Arrigo, Symmetry Analysis of Differential Equa- tions: An Introduction, Wiley (2015) 73-84 Biomath 1 (2015), 1503081, http://dx.doi.org/10.11145/j.biomath.2015.03.081 Page 14 of 14 http://dx.doi.org/10.1137/040605515 http://dx.doi.org/10.1146/annurev.matsci.32.101901.155803 http://dx.doi.org/10.1146/annurev.matsci.32.101901.155803 http://dx.doi.org/10.1016/j.apm.2007.10.006 http://dx.doi.org/10.1007/s11425-013-4635-8 http://dx.doi.org/10.1016/j.na.2007.01.053 http://dx.doi.org/10.1016/j.na.2009.03.087 http://dx.doi.org/10.1073/pnas.0802592105 http://dx.doi.org/10.1007/BF01444068 http://dx.doi.org/10.1016/j.jmaa.2010.06.026 http://dx.doi.org/10.1016/j.amc.2013.06.034 http://dx.doi.org/10.1016/j.apt.2012.04.001 http://dx.doi.org/10.11145/j.biomath.2015.03.081 Introduction The Two-Phase Model Symmetry Analysis Recovering the Exact Solution for a Free Boundary Problem Other Solutions to (25) Additional Solutions to the Two-Phase Model Discussion Appendix References