A computational study of fluctuating viscoelastic forces on trapped interfaces in porousmedia

In immiscible liquid–liquid Newtonian flow through a porous medium, one phase often becomes trapped in corners or narrow regions by capillary forces as blobs (e.g. oil ganglia) deep within the matrix, whilst the second flow exhibits a steady and laminar flow (e.g. of water) that has negligible influence on the trapped liquid–liquid interfaces. However, recent microfluidic experiments have shown the situation radically changes when using a viscoelastic liquid, which is capable of exhibiting pore-scale unsteady flow that can deform such interfaces. Here, a computational model is developed which allows us to capture the forces that cause this behaviour and provide a framework for future investigations of this system. In this paper, the forces on trapped interfaces are investigated for the first time. Notably, when the viscoelastic flow becomes unsteady the forces on the trapped interfaces not only fluctuate but also become amplified, thus supporting experimental findings showing they can be used to free such interfaces. At Weissenberg values of 1.5 and above the fluctuations become important and the mean values of the forces on the interfaces decrease as the fluctuations grow. © 2020 The Authors. Published by ElsevierMasson SAS. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
During two phase immiscible flow within a porous media, one liquid can become trapped due to capillary forces. This phenomenon is particularly relevant in the oil recovery process, where in the secondary phase [1] a Newtonian fluid (typically water) is pumped in to remove oil that is left after primary recovery. However, after this stage a considerable amount (on average 60%-80% [2]) of oil remains trapped in the form of blobs ('ganglia' [3]) or is bypassed by the Newtonian fluid. Further pumping of the Newtonian fluid has no influence, as a steady laminar path is formed and the flow remains within channels that bypass the ganglia. Therefore, in the tertiary phase, 'Enhanced Oil Recovery' techniques are required that exploit additional physical mechanisms.
The use of polymeric (non-Newtonian) fluids during enhanced oil recovery is common to increase the viscosity of the displacing fluid and exploit elastic properties of polymeric fluids in order to sweep the bypassed regions. Different mechanisms for enhanced oil recovery have been proposed, but there is little quantitative evidence for the particular forces at play [4]. The forces that are responsible are challenging to measure experimentally due to the difficulties of resolving phenomena 'inside' solids on small spatial and temporal scales. In contrast, such measurements can easily be obtained from computational models of the process. * Corresponding author.
One proposed mechanism for enhance oil recovery is shear thickening, where the viscoelastic fluid experiences elongation and contraction while passing through a porous medium, so a high apparent viscosity is formed as the velocity is large enough to keep the polymer molecules stretched. Recently, the shear thickening mechanism has been observed experimentally in a single phase flow through a porous medium as an increase in pressure gradient, or apparent viscosity of the fluid on the Darcy scale. On the pore scale, the increase in apparent viscosity was shown to coincide with the onset of flow fluctuations [5,6]. The unsteady flow behaviour is known as 'elastic turbulence' as the observed behaviour is crudely similar to observations of inertial turbulence, however the unsteady behaviour is caused by the elastic properties of the fluid rather than inertial forces (which can be considered negligible at small length scales). Polymer solutions that display shear thickening behaviour can displace more oil from a rock core sample than solutions that do not, however the increase in oil recovered cannot be accounted for by the increase in pressure gradient alone [6]. Similar effects have been observed experimentally within pillared microchannels by De et al. [7]. The mechanism behind this additional oil displacement from rock samples may be caused by the flow fluctuations but it is challenging to identify due to the difficulties of visualising the pore space at the spatio-temporal scale of interest.
To overcome the challenges that arise from using rock samples, quasi-2D microfluidic analogues have been used to investigate the properties of the flow in porous media [5,6,8]. These  microfluidic geometries make it easier to visualise the flow behaviour experimentally, while being able to replicate some of the key features of displacement from a real rock sample. Notably, oil becomes trapped as ganglia in the pores of the matrix, see Fig. 1, as it would in the pores of a limestone rock sample, for example [6]. The trapped oil ganglia cannot be displaced by a Newtonian flow. Two-phase flow experiments with microfluidic networks have shown polymer solutions that exhibit a global increase in apparent viscosity correspond to cases when the pore scale flow exhibits an unsteady chaotic flow and oil droplet interfaces fluctuate, i.e. elastic turbulence is responsible for the increases in apparent viscosity [6,8]. The onset of these fluctuations in the flow corresponds to a jump in the amount of oil removed from a rock core sample [6].
Due to the small length scales in porous media and the relatively slow flows usually encountered, the Reynolds number is small, and so these time dependent effects can only be caused by the elastic properties of the fluid. This hypothesis is supported by the fact that the displacement amplitude of the oil droplets increases with larger flow rates that correspond to higher elastic forces [8]. However, a firm understanding of the mechanism driving the fluctuations and their influence on displacing trapped interfaces is still lacking. Computational modelling can improve our understanding of such mechanisms as it enables analysis of forces on the meniscus that cannot be obtained experimentally. Here, we will combine computer modelling with two flow configurations to investigate the effects of viscoelastic fluctuations. Studying the applied forces computationally will clarify which mechanisms contribute most to enhanced oil recovery. This will create a framework and preliminary results for future work, where we can extend the model to two phase flow with free surfaces, optimising the properties and driving forces of the displacing fluid for oil recovery and relating the pore scale dynamics to Darcy-scale dynamics.
Viscoelastic flow instabilities have been shown to occur in various geometries [9][10][11], but here we focus on three configurations of particular interest. Firstly, a cross slot geometry, where four rectangular channels meet at a 'pore' with the two inlets/outlets opposite each other in order to generate extensional flow, see Fig. 2A. This is a classical configuration which will be used for benchmarking as it has been studied in relevant conditions both experimentally [12][13][14] and computationally [15][16][17]. Secondly, the start up time of flow in a channel used to validate the time dependent component of the flow with an analytic solution. Thirdly, we consider a single symmetric pore based on the microfluidic network used by Clarke et al. [6], Fig. 1, with the two inlets/outlets adjacent to one another, see Fig. 2B.
In the cross slot geometry, two transitions in flow behaviour have been observed; firstly, from steady symmetric to steady asymmetric flow and then secondly from steady to time dependent flow [12]. Cruz et al. [17] propose that the first transition to steady asymmetric flow can occur when Wi is sufficiently large because there is streamline curvature and a free stagnation point -so the strain rates at the stagnation point are non-zero. These flow behaviours have been shown to be dependent on three dimensionless parameters: • Weissenberg number, Wi, the ratio of elastic to viscous forces [12,18]; notably, Sousa et al. [19] found that the time dependent flow behaviour was periodic at smaller Wi values (e.g. Wi = 6.6) and the flow became less regular as Wi increased (e.g. Wi = 282). • Viscosity ratio, β, which is the ratio of the solvent to total (solvent and polymer) viscosity [16]; as β increases so does the value of Wi required before a transition is observed.
• Aspect ratio of channels, ε = H/W , where H is depth and W is width; the smaller ε is, the larger Wi must be before a transition is observed [14]. Indeed, for shallow geometries (ε < 1) Cruz et al. [20] showed computational results using Upper Convected Maxwell and simplified linear Phan-Thien and Tanner (sPTT) models that predict no asymmetric state and the transition to time dependence occurs at Wi close to 0. Steady asymmetric flow was only observed for ε > 1. For computational 2D cross slot geometries, where ε → ∞, it is expected that both transitions will be observed and will occur at smaller Wi values than measured in quasi-2D or 3D experiments.
The channel geometry is used to validate the time dependent behaviour resulting from the model set up, following the approach of Zhang et al. [21]. An analytical solution exists for the start up time of Poiseuille flow for the Oldroyd-B model [22] and this is compared to the results from the numerical model. An instantaneous spatially uniform pressure gradient is applied at t = 0 to drive the flow.
The third geometry is chosen as oil can become trapped in corners of the central pore as shown in Fig. 1. It is a representative pore from the microfluidic network shown in Fig. 1, used to investigate fluctuating interfaces in Clarke et al. [6] with two adjacent inlets, see Fig. 2B. This has a square central pore with two inlets and two outlets, one centred on each side of a square pore with the inlets/outlets adjacent to one another (i.e. not opposite, as in the cross-slot). Only one transition, from steady to unsteady flow, is expected in the microfluidic geometry, because there is no free stagnation point for the flow inside it. The aspect ratio of the central pore in the microfluidic geometry in [6] is ε = 0.5 and the aspect ratios of the channels is between ε = 1 and 2.
These are not considered shallow, therefore, as a starting point we will use a 2D geometry. As demonstrated by Cruz et al. [20], the effect of this assumption is that the critical value of Wi where the fluctuations begin will be higher in the 2D case than in the more accurate 3D case. Within this geometry, the forces acting on the trapped drops will be considered in order to understand better the mechanism for drop displacement due to viscoelastic fluctuations.

Problem formulation
The viscoelastic model is coupled with the incompressible Stokes flow within a domain, Ω, to give dimensionless equations, where τ is the extra stress tensor for the viscoelastic component of the flow, u is the velocity of the fluid and p is pressure. The viscosity ratio is β = η s /η 0 , where η 0 = η s + η p is the total viscosity from adding the solvent and polymer viscosities. We non-dimensionalise p and τ by the same scale, η 0 U/L where U is the characteristic velocity -the average flow velocity of steady Newtonian flow in a channel, and L is the characteristic lengththe width of the inlet channel. The extra stress, τ, is related to the conformation tensor, C , (the ensemble average of QQ , where Q is the end-to end vector of the dumbbell), by where the function F (C) depends on the constitutive model used.
The equation for C , is given by, where is the Oldroyd derivative. There are numerous models of varying complexity that can be used for F (C ) to capture viscoelastic behaviour [23]; here we focus on the simplest, well-known Oldroyd-B model with just two parameters, a constant shear viscosity and relaxation time. The Oldroyd-B model can be derived from a molecular model of Hookean dumbbells in a Newtonian solvent, assuming that the solvent velocity field is homogeneous and one can neglect hydrodynamic interactions and external forces, as stated in Owens and Philips [23] (pgs. 33-37), so that: where Wi = λU/L is the Weissenberg number with λ the relaxation time. Wi is the ratio of elastic to viscous forces. I is the unit tensor. The boundary conditions are no-slip, on the solid walls of the domain whilst the inlet/outlet conditions depend on the particular flow configuration considered.

Extensional flow in a cross slot
The model is validated in 2D using the cross slot geometry, shown in Fig. 2A. The inlet conditions for the velocity field are fully developed Poiseuille flow with mean velocity 1 in a channel from y = −0.5 to 0.5.
Then, assuming steady flow, the viscoelastic stress is calculated analytically, At the outlets, we assume the flow is parallel to the walls, that there is no normal stress generation and the pressure is zero: u u u · t t t = 0.
where n n n is the normal vector pointing out of the domain and t t t is the tangent vector. The domain for the cross slot geometry is tessellated with a square mesh, Fig. 3A, which is equivalent to the finite volume mesh used in Cruz et al. [17], i.e. it has the same number of elements of the same shape but an element is placed at the centre of the cross slot, rather than a volume. For future work with more complex geometries a triangular mesh would be more flexible and so the domain is also meshed with triangular elements, Fig. 3B, to see how this effects the results. For the triangular mesh finer elements were placed at the centre of the cross to ensure the behaviour around the stagnation point was captured and one layer of rectangular elements was placed along the walls to capture boundary effects. Whether the domain has square corners or filleted (the two edges either side of the corner are joined by a curve) corners has been shown to have little effect of the final solution [16,17]. This finding is confirmed here by using square corners with the square mesh and filleted corners with the triangular mesh.

Start up time for poiseuille flow
Miranda and Oliveira [22] present results showing the start up time of the viscoelastic flow in a channel. We used this method to validate the time dependent behaviour of the model. The boundary conditions are no-slip on the walls and periodic flow conditions applied on the inlet and outlet. Start up flow is initiated by applying a constant pressure gradient at t = 0.
For the Oldroyd-B model, the velocity was found analytically by Waters and King [24] as, where F n (y) = sin (Z n (1 + 2y)) , Z n = π 2 (2n − 1) , and for β 2 n = α 2 n − 16WiZ 2 n , α n = 1 + 4βWiZ 2 n , Note that the equations have been rescaled to match the dimensionless equations presented here. The channel for the numerical model has width L and length 10L and therefore a pressure gradient of 120 is applied to ensure an average velocity of 1 at steady state.

Microfluidic network geometry
The final model we will investigate is to apply the coupled Stokes and viscoelastic equations to a geometry based on one pore from the microfluidic network used by Clarke et al. [6]. A 2D geometry is used to reduce the computational requirements of the model and the geometry is simplified by using equal width inlets and outlets rather than the various widths used in the experiment. The geometry used is shown in Fig. 4. The effect of the flow on a fluid trapped in a corner of the domain is modelled by rounding one of the corners of the central pore, either the corner opposite both inlets, shown in Fig. 4A, the corner between the inlets, Fig. 4B, or the corner between one inlet and one outlet, Fig. 4C (due to symmetry, the results from the bottom left hand corner would be equivalent). As a first step, in this article we will not introduce a second liquid and associated free boundary, but rather consider the forces on a rigid surface with a no-slip boundary condition in the corner. The boundary condition on the walls for the viscoelastic fluid is no-slip. Periodic boundary conditions are applied at opposite boundaries for the velocity and viscoelastic stress. A pressure drop of 117 is applied, which ensures the average velocity is 1 in the channels for the Newtonian case.

Numerical method
The equations are implemented in COMSOL Mulitphysics (v5.2a, Sweden), a commercial finite element software, using the Discrete Elastic-Viscous Split Stress formulation [25], to obtain a mixed finite element formulation, and the log-conformation method [26] to improve the convergence of the numerical method at higher Wi. The implementation of this method for viscoelastic fluid flow using finite elements has been described elsewhere (e.g. [27]) and we do not repeat it here. For all the calculations quadratic elements are used to solve the velocity components and linear elements are used for the pressure and extra-stress components. The Newton method is used with the backward differentiation formula for the time stepping with the 'free' option in COMSOL for adapting the size of the time steps taken. The solution is output every △t = 0.1 from t ∈ (0, 500). The

Extensional flow in a cross slot
The computational model is validated using the cross slot geometry, where the point of transition to asymmetric flow (from the symmetric Newtonian one, c.f. Fig. 5) and level of asymmetry have been previously measured by an asymmetry parameter DQ . The asymmetry parameter means that the transition of the flow to asymmetry can be quantitatively identified and compared to other models. The value of DQ is calculated as where Q 1 is the flux from inlet 1 that leaves through outlet 1 and Q 2 is the flux from the same inlet that leaves from outlet 2, see Fig. 2A Fig. 5D, with huge fluctuations observed as one jumps between the different flow configurations. Clearly then, with increasing Wi a highly nonlinear unsteady flow can be generated in our simulations, as observed experimentally [19].

Start up time for poiseuille flow
To validate the time dependent results of the model the results for channel flow are compared to the analytic solution [22]. In Fig. 6 the centre line velocity resulting from the analytic solution is compared to the centre line velocity, u 0 , from the model for Wi = 1 and Wi = 2 with β = 0.1 for both cases. The model results show good agreement with the analytical solution and so the model is valid for Wi ≤ 2, which is the range we will study in the next section.

Microfluidic network geometry results
Consider now the representative pore of the microfluidic network, where Fig. 7 shows that for steady cases (Wi = 0), the flow is diagonally symmetric. For Stokes flow and steady flow at low Wi the streamlines are mirrored along the diagonal with all the flow from the top inlet leaving through the right outlet and all the flow from the left inlet leaving through the lower outlet. As Wi is increased and the flow becomes time dependent, asymmetry is lost and some of the flow changes direction to leave through  The interfaces of trapped oil ganglia have been observed to fluctuate in the presence of a viscoelastic shear flow [6]. In the following analysis we consider forces on the oil droplet surface and infer how the interface may be deformed. If sufficient force is applied to the meniscus the capillary forces that hold it in place can be overcome and the oil drop can be released, this is beyond the current study and will be considered in future work. Two perpendicular forces are calculated on the drop trapped in a corner of the domain. The force out of the corner, F out , and the force across the corner, F across . Both forces are calculated by integrating the stress along the droplet surface, ∂Ω d , with F out = ∫ By plotting the mean and standard deviations of F out and F across it can be seen that as Wi increases so do the magnitude of the fluctuations, as quantified by the standard deviation, Fig. 10A. The mean force out of the corner is negative for the upper right corner so the drop will be pushed into the corner, however, for Wi ≤ 0.6 the across force has a larger magnitude. This implies that for Wi ≤ 0.6 there is a competition between the force across the drop, which is acting to destabilise it, and the force into the corner which would appear to be inhibiting drop detachment. The mean force out of the upper left corner is always positive, so the interface would be pulled out of the corner. The lower right corner is most interesting, as the mean force out of the corner is negative apart from in a region around 1.25 < Wi < 1.75 where it becomes positive, suggesting in this range the drop may be pulled out. Notably, this range also corresponds to maxima in the mean force out of the other corner, suggesting a potential sweet spot for displacing oil, if one is to believe the mean force out of the corner is the most important quantity. The standard deviation in Fig. 10 shows that the magnitude of the fluctuations increases as the Wi value increases. The magnitude of the mean force out of the corner on all the corner droplet surfaces appear to be a maximum for a value between Wi = 1.25 and 1.75. The mean magnitude of the force across the corner, Fig. 10B, stays around zero for the lower right and upper left corners, as a result of the symmetry of the flow and any deviation from zero is due to limitations on simulation time due to computational intensity. These results suggest that the mean force across the drop may not be the most important quantity, as it appears to show little sensitivity to Wi.
In Fig. 11, the force out of the corner is separated in to separate components: pressure shown by the standard deviation.

Discussion
In this paper we have validated a finite element model for viscoelastic flow implemented in COMSOL Multiphysics and shown that the results compare well with literature. The computational model is used to gain insight into the forces generated on droplets that become trapped in the corners of a porous network and can be displaced by unsteady viscoelastic flow. At sufficiently high Weissenberg numbers, Wi > 1, the streamlines show that the flow from either inlet can change direction to exit through either outlet. This effect may be similar to that seen in De et al. [28] where increasing the value of Wi caused the flow to 'change lanes' in a periodic porous geometry. This may prevent preferential pathways from forming and therefore less oil is bypassed.
In Section 4, the forces that arise on the surface of droplets in three corners are compared. Here we will infer how these forces may deform the interface. In the lower right corner the mean force out of the corner pushes the drop into the corner for Wi < 1.25 and Wi > 1.75. The main contribution to the force is from the pressure term. In the upper left corner the mean force  out of the corner pulls the drop away from the corner towards the centre of the pore. In Fig. 1B, some of the drops show elongation in the corner that corresponds to the upper left corner of the modelled geometry. This are indicated by arrows. This may be an example of the so called 'pulling effect' [4]. The mechanism appears to be that the elastic components cause the asymmetry and this builds up large pressure gradients within the flow. In the upper right corner the mean force out of the corner pushes the drop into the corner and the across force pushes the surface of the drop towards the right hand outlet. The magnitude of the forces on the upper right drop are the largest of all three drops. Therefore, the model suggests this is the corner where the least oil is trapped. This appears to be the case on experimental images in Fig. 1B, as shown by the black arrows. However, this is also the case for the low Wi flow in Fig. 1A. The force out of the upper right hand drop is only due to the elastic forces as for Wi = 0 the mean force out of the corner is zero. This may be another example of the pulling effect. Fig. 10A shows for the upper left and upper right corner drops the magnitude of the mean force out of the corner is increased for Wi > 1 when compared to Wi = 0. This may be an example of the 'stripping effect' [4].
Currently, a larger value for the viscosity ratio is used in the model compared to the experiments. The effect of changing the viscosity ratio has been shown in Rocha et al. [16] where the transition to asymmetric flow appears at smaller Wi as the viscosity ratio is decreased. However there appears to be a lower limit of Wi ∼ 0.31 [20,29]. Different viscoelastic models also result in a transition to asymmetric or time dependent flow at different Wi numbers [17] and therefore it is unnecessary to be overly accurate about this value until a particular model can be shown to closely compare to experimental values.
The Oldroyd-B model is used here for its simplicity. A disadvantage of the Oldroyd-B model is that the Hookean dumbbells are assumed to be infinitely extensible which can cause the model to blow up at finite extensional rates [23]. Due to this assumption it can be difficult to model numerically as large stresses can occur. However, now that suitable meshes have been obtained for the hardest case of Oldroyd B simulation, it is simple to adapt the implementation to use similar viscoelastic models. It is not clear which viscoelastic fluid model is most appropriate for the extensional flow in a cross slot, as no model so far has been matched directly to an experiment. Experimentally, the transition to time dependent flow in a cross slot occurs at values of Wi ≥ 2.5 [12][13][14] where as computational studies show the transition to time dependent flow occurring for Wi < 2 [20]. Clearly, there is scope for further collaborations between experiment, modelling and computation for this complex flow. This paper has presented preliminary work towards investigating the effects of viscoelastic fluctuations on trapped drops. In the future this work will be extended to include free surfaces so that the deformation of the drops, due to the applied forces, can be visualised. Other viscoelastic models, such as the finitely extensible nonlinear elastic, sPTT or Rolie-Poly models, will be considered to compare how this effects the forces on the drop and discover whether the results are more realistic. The small scale behaviour will be used to generate Darcy scale models that will aim to predict how the viscoelastic fluctuations effect macroscale variables of interest such as flow rate, pressure gradients and oil saturation. The model may be extended to three dimensions and combined with 3D imaging data (e.g. X-ray computed tomography) so that the deformation of the trapped oil ganglia can be compared with experimental data in real rock core sample geometries, such as the data presented in Datta et al. [30,31]. Furthermore, the model may be used to investigate phenomena in heat transfer applications and CO 2 sequestration.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.