A Particular Solution for a Two-Phase Model with a Sharp Interface

Two-phase models can be used to describe 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 modelling; 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 material.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 interactions 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 dragfriction 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 parameter 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 overlapping filopodia.Whereas, poorly differentiated (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 conditions.Simulations also provided a critical region for which a necrotic core forms at the tumor's center.
Several forces are required to balance conservation of momentum.For two-phase models, the viscosity of the phases and interstitial friction must be accounted, but for biofilm morphology, in addition to hydrostatic pressure, osmotic pressure is also needed.One such model [4] describes the role of a network comprised of an extra-cellular polymeric 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 simulate biofilm growth and cell motility [5], [6].A mobile cell contains polymer network phase comprised 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 useful to track the sharp interface.The location of the interface can be followed explicitly by interface 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 perturbation 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 discussed, 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 transformations which can reduce a system of Partial Differential Equations (PDEs) in one spacial dimension to a system of ODEs.These symmetries are generated by introducing infinitesimal transformations, which leave the original system invariant.For classical symmetry analysis, expansion of this infinitesimal transformation, produces a linear system 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 transformations give rise to highly non-linear determining equations and can be difficult or impossible to 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 traveling 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 describe 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 negligible.Gels are composed of a polymeric network given by φ 1 and a fluid solvent φ 2 .Both phases are treated as Newtonian fluids that are immiscible.When considering the redistribution of mass within a control volume, the flux of the network is given by ∇ • (φ 1 u 1 ), where the network moves with a velocity u 1 .A similar argument is made for the solvent to give the following equations to conserve mass.
where the sum of the volumes saturate to a fixed control volume, 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.
where μ1 is the shear viscosity and λ 1 is the bulk viscosity.In 1-D, this becomes where Another force that we include is the frictional force created by interstitial interactions between phases.If both fluids move in unison or if either volume fraction becomes negligible, drag will vanish.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 ) = k 2 φ 2 1 (φ 1 − φ 0 ).The constant k 2 accounts for the effects of the ionic environment, polymeric structure, and solvent concentration that contribute to swelling and deswelling.The value of φ 0 is a reference volume 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.
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.
where µ 2 is the viscosity of the solvent.Summing (4) and ( 5) gives the following equation.
Together equations (8-10) can be reduced to a system of ODEs using the following transformation.
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 transformation 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 transformations 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 unchanged.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.
where Φ 1 , T , X, U , and V are called the infinitesimals.In general, one seeks to find invariance of a system of differential equations of the form F i (t, x, u, v, φ 1 , u t , v t , φ 1t , u x , v x , φ 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 F i is given by the equations (8)(9)(10).Under (12), a set of differential equations is produced for the infinitesimals T , X, U , and V .These differential equations are called the determining 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.
When the infinitesimals are solved in conjunction with the invariant surface condition given by (14-16), the solutions u, v, and φ 1 provide a transformation 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 interactively to give the forms of the infinitesimals.
Due to the size of the equations, details of the determining equations are omitted.For more 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 δ = 0.If δ = 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 (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 symmetries, as was shown with the KdV equations, Camassa-Holm, and sine-Gordon equations [23].Non-classical symmetries have also been generated 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 symmetries.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 where u, v, and φ 1 are all functions of t, x as previously discussed and P x is the pressure gradient.Next, we solve (19) for P x to find We see the mass equations (1-2) have now become Summing these two equations of mass gives Imposing the average velocity is zero, we have 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 Multiphase φ 1, φ 2 > 0 (20-21), reduces the equation ( 18) to the following equation.
As in Figure (1), we assume the mixture occupies the interior region (Ω 1 ), while pure solvent occupies the external region (Ω 2 ).As the gelmixture swells/deswells, the interface between the regions (Γ) moves.To specify the motion Du et.al. impose standard jump conditions: The solution found by Du et.al. [10] approximates the front velocity for the free boundary problem at time t = 0.The solution for a piecewise constant profile is given by, and the following can be derived where , and 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.
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 condition 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.
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) 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 attempt 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).
For µ 1 of a much larger magnitude than k 2 and ξ, this becomes the following: whose solution is This makes the analytic solution for the original system (1-2) and (4-5) to be 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: 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 becomes 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 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.
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.There are several benefits of finding the exact solution, instead of using numerical methods.First, numerical results have a difficult time capturing 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.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 solutions are extraneous, but more likely, additional solutions describe other physical or biological phenomenon 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 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 k 2 = 0 as seen in the absence of ionizing agents for gels, we find the following solution 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 travels as a wave moving to the right.This does not seem to have any physical analogue since the velocities are unbounded.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 example, 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 considering 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.
Neglecting solvent viscosity, µ 2 = 0, reduces the original system (8-10) to the following ODE, Again, if we assume the dominating force is the viscosity and set k 2 and ξ 2 to zero, this can be solved to give The complete solution to (1-2) and (4-5) becomes 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.This may imply that viscosity is required for nonconstant 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 which reduces (8-10) to the following three ODEs.It should be noted that if either viscosity is the dominating force with k 2 = 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 solution requires a dominating force to generate nonconstant solutions.This gives flexibility in assessing 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 implementation 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 behavior of the additional solutions may give further

Fig. 1 .
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.
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.

Fig. 4 .
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.

Fig. 6 .
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.