Stationary dimpled drops under linear flow

Abstract The axially symmetric deformation of a drop in a viscous fluid, under the influence of an externally imposed flow having simultaneous rotating and compressional or extensional components, is addressed. In the previous studies, two families of stationary drop shapes were constructed by simulating the dynamics of drop deformation: stable singly connected shapes with respect to axisymmetric disturbances, and unstable toroidal shapes. These two branches coexist at the same flow conditions, but were not connected. In this study, we obtain a new family of branches of unstable highly deformed stationary drops connecting with the stable flattened shapes and the toroidal ones. We use a method based on classical control theory. The controller is designed for a two-state dynamic model of the system and is employed on a high-order nonlinear dynamic model of the drop deformation. Through this feedback-control-centred approach, an extended collection of unstable stationary solutions is constructed, which spans the range from the loss of stability to the dimpled shapes almost collapsed at the centre. In the latter region, which was never obtained in previous studies, a multiplicity of solutions is identified.


Introduction
In recent years, interest in non-trivial forms of fluid particles was stimulated by applications of non-spherical microparticles to have important potential as building blocks for self-assembled materials, including clustering of cells, imaging probes for therapy, drug carriers, and more (see e.g.Champion, Katare & Mitragotri 2007;Yan 2015;Zabarankin & Grechuk 2017).One of the advanced methods for producing microparticles of complex shapes is the solidification of drops deformed by the flow in microfluidic devices.The flow affecting the deformation of a microdrop, can be approximated by a linear one.The resulting linear flow can have shear and rotational components.Systems that may involve such simultaneous components can be found also in multiphase polymer processing such as injection moulding or cases where two impinging jets collide.Further examples of industrial devices in which domains that have a combination of such rotation and extension components are rotary pumps Dmytriv et al. (2021) and cyclone separators Li, Huang & Li (2023).
There has been an enormous body of publications concerning the deformation of a drop embedded in a viscous domain and subject to a linear flow field.A comprehensive review of the early studies can be found in Stone (1994).In particular, when a viscous drop, initially spherical at rest, is deforming in an extensional or shear flow, the works of Taylor (1932), Rallison & Acrivos (1978) and Hinch & Acrivos (1980) suggest that with increasing flow rate, the drop can reach a critical shape at which it loses its stable form and beyond which its shape becomes unstable in the flow.A similar transition is found by Acrivos & Lo (1978) for a slender body in extensional flow where more than one branch of the solution of the equations of motion, and more than one stationary shape, stable and unstable, are reported for given flow conditions.
The study of rotating liquid drops has been performed by several authors experimentally, theoretically as well as numerically (see e.g.Chandrasekhar 1965;Brown & Scriven 1980;Fontelos, García-Garrido & Kindelán 2011;Lyttleton 2013;Nurse, Coriell & McFadden 2015;Elms et al. 2017;Holgate & Coppins 2018).The computational study of equilibrium shapes of such drops was conducted by Poincaré (1885) and Appell (1932).According to Poincaré (1885), two families of solutions could be obtained, which consist of axisymmetric and asymmetric shapes.A detailed classification of all rotationally symmetric figures of equilibrium corresponding to rotating liquid masses subject to surface tension was given in Elms et al. (2017).Fontelos et al. (2011) studied the stationary shapes, evolution and breakup of a viscous rotating drop, making use of the Stokes model and a boundary integral approach.Multiple axisymmetric solutions were found for a range of governing parameters.Ellipsoidal, flat dimpled and toroidal shapes may exist for the same flow conditions.
The multiplicity of stationary shapes of drops deforming in a linear compressional flow was also reported in Zabarankin et al. (2013), Zabarankin, Lavrenteva & Nir (2015) and Ee et al. (2018).There, the reported branches of solutions were either stable stationary flat drops or unstable yet stationary toroidal drops.These predictions comprise dual branches of stationary solutions, resulting at a given capillary number, but they appeared to be unconnected, in contrast to the drops deforming in a rotational flow field, where the stationary branches of flat and toroidal drops are connected by a branch of flat discs having ever-growing dimples (Fontelos et al. 2011).Lavrenteva et al. (2021) proposed to approximate stationary shapes of the drops in compressional flow by using generalized Cassini ovals.The analysis reproduced the branches with shapes of stationary stable flat drops and stationary unstable toroidal drops, available from numerical calculation, and predicted the transition branch between oval and toroidal shapes, consisting of flat drops with ever-growing dimples up to collapse.Malik, Lavrenteva & Nir (2020) considered deformation of an immiscible drop in a linear flow that is composed of two simultaneous components, rotational flow and axisymmetric flow, either compressional or extensional.The dynamic deformation, stationarity and stability with respect to axisymmetric perturbation were characterized using dimensionless parameters such as the Bond number and capillary number.A domain in the phase plane of the Bond number and the capillary number within which the deformed drop will be stable is obtained.Malik, Lavrenteva & Nir (2021) reported on the branch of toroidal stationary shapes existing under the same conditions as simply connected ones.Malik et al. (2022) used a feedback-control approach to extend the branch of toroidal drops up to near collapse.Yet the branches of singly connected and toroidal shapes remained unconnected in cases with and without rotation.
The main goal of this study is to construct the connecting branch in the case when the outer linear flow consists of rotational and axisymmetric (extensional or compressional) components.The method that was employed is based on the use of classical control theory similar to that employed in Malik et al. (2022).Methods of control theory are employed extensively and successfully to stabilize various types of flow, e.g. for controlling the onset of turbulence, and turbulence suppression in channel and pipe flow.A comprehensive review of the modern state of this subject can be found in Jovanović (2021).Fikl & Bodony (2021) developed optimal control methods for stationary and time-dependent deformation of a drop in viscous flow.
In the present study, a simple proportional-integral controller is used.The intensity of the extensional component of the outer flow is chosen for a control signal.The controller is designed using a simplified model of the system and employed for the full model, demonstrating sufficient robustness.In § 2, we give the mathematical model used for the analysis.The design of the controller based on a simplified model of the process is described in § 3. Section 4 presents the implementation of the controller to the full (boundary integral) dynamic problem.A discussion of the results concerning stationary solutions is given in § 5, followed by concluding remarks in § 6.

Drop dynamics in linear viscous flow
We consider a viscous drop of fixed volume V = 4πl 3 /3 embedded in an immiscible unbounded viscous fluid, where l is the radius of the undeformed spherical drop.Let the domain occupied by the drop be denoted by Ω 1 , and let the unbounded external domain be defined by Ω 2 .The viscosity and density of the drop and the ambient fluid are represented by μ 1 , ρ 1 and μ 2 , ρ 2 , respectively.The velocity and pressure fields in Ω i , i = 1, 2, are denoted as u i , p i , respectively.The drop and external fluid are supposed to rotate in an axisymmetric manner around a common axis with angular velocity ω.
In this study, the external flow is considered to be extensional or compressional in combination with rotation.In the reference frame rotating with the angular velocity ω, the velocity of external flow, in the absence of a drop, is where Here, G denotes a shear rate that characterizes the intensity of the flow, which may be positive or negative, representing the compression and extension modes of the flow, respectively.
In what follows, throughout the paper, non-dimensional variables are used, where l, U = γ /μ 2 , Π = γ /l and T = l/U are the length, velocity, pressure and time scales, respectively.The governing parameters of the problem are the viscosity ratio λ = μ 1 /μ 2 , the capillary number Ca = μ 2 Gl/γ , and the rotational Bond number Bo = (ρ 1 − ρ 2 )ω 2 l 3 /γ .In the present study, the drop and the ambient fluid are considered to have the same viscosity, i.e. λ = 1.
In this investigation, we explored the realm of creeping incompressible flow conditions, where the viscosities of the drop and the ambient fluid are considered to be same.
When inertia and sedimentation effects are ignored, the Stokes system in a rotating frame of reference is given as Here, P i is the non-dimensional form of modified pressure as given in Fontelos et al. (2011).The stress balance at the drop interface is maintained as where the σ i denote the interfacial stresses, n is the unit normal vector pointing outwards from the drop interface, and κ is the surface mean curvature.The non-dimensional parameters Ca and Bo are the capillary and Bond numbers, respectively.At the interface of the drop (δΩ), the velocity (u s ) is assumed to be equal to the velocity of the ambient fluid, i.e. u 1 = u 2 = u s . (2.5) The dynamics of this surface is determined by the kinematic condition where u n denotes the velocity of the surface in the normal direction.Far from the drop, It should be noted that the parameters Ca and Bo can assume positive as well as negative values.Here, the positive and negative values of Ca represent the compressional and extensional flow, respectively.When Bo is positive, the drop density exceeds the density of the ambient fluid, which we call a heavy drop; and when Bo is negative, the density of the ambient fluid exceeds the drop density, which we call a light drop.

Boundary integral formulation
For a given shape of the drop, the problem (2.2)-(2.5),(2.7) can be formulated as a system of integral equations for the velocity at the interface ∂Ω as given by Rallison & Acrivos (1978) and Pozrikidis (1992).In the case of equal viscosity of the phases, the velocity at the interface is given by an explicit expression, where x p denotes the position vector of any point p on the surface, and Recall that u ∞ (y) = ( y 1 , y 2 , −2y 3 ), and the positive and negative signs of Ca refer to compressional and extensional flow, respectively.Here, δ ij is the Kronecker delta.In this study, y 3 = z, and the z-axis is considered to be the axis of rotation.The pressure and velocity are independent of azimuth angle φ in the cylindrical coordinates (r, z, φ).Equation (2.8) can be integrated over the angular coordinate and reduced to expressions for the components of the velocity at the interface in cylindrical coordinates containing curvilinear integrals over the cross-section of the drop interface L: (2.10) The resulting expressions for the kernels M ij can be found in Pozrikidis (1992).The simulation process of dynamic deformation starts with the solution for some given shape of a singly connected drop of volume 4π/3.At each time step, the velocity at the interface is computed from (2.10), and the drop's surface is updated accordingly in the normal direction in a quasi-stationary manner with the help of the kinematic condition (2.6).The detailed literature concerning this approach can be found in Zabarankin et al. (2013), where it was used for simulating the dynamic deformations of a drop in a compressional flow in the absence of rotation, Ca > 0, Bo = 0.The same approach was employed recently by Malik et al. (2020) to study deformation of a drop under the combined action of rotation and compressional or extensional flow.It was demonstrated that for every pair (Bo, Ca), the drop either stabilizes to a stationary shape or eventually breaks up.

Stationary drop shapes
The drop is of stationary shape if the normal velocity at its interface vanishes, that is, (2.11) In the dynamic computations (Zabarankin et al. 2013;Malik et al. 2020), the solution was assumed to be stationary if (2.11) was established to a desired degree of accuracy.We note that, following Malik et al. (2020Malik et al. ( , 2022)), with a mesh of 100 surface elements, the criterion for stationarity is set at max(|u n |) ≤ O(10 −3 ).The obtained singly connected stationary shapes do not change for an indefinitely long time, and are addressed as stable with respect to axisymmetric perturbations.The region of stable drop shapes in terms of Bo and Ca has been presented in Malik et al. (2020).
In the special case when Ca = 0, the velocity field in the laboratory reference frame is just that of a solid body rotation, and the shape of the interface can be found via solving a system of ordinary differential equations (see e.g.Fontelos et al. 2011).Analysis of these equations reveals that the solution is not unique.For a certain region of Bo, oval shapes coexist with toroidal and flattened, dimpled ones.
An alternative method to find stationary solutions is to use an iterative procedure to minimize a norm of the normal velocity at the interface.For a stationary shape, this minimum should be zero.This method is not necessarily related to the actual dynamics of the drop and can be advantageous in finding unstable shapes.In particular, the effectiveness of such an approach is anticipated when the shape of the drop can be approximated by a simple geometric shape defined by a small number of parameters.Numerous such approximations are available in the literature.Zabarankin et al. (2013) suggested several such approximations in the case Bo = 0 and λ = 1, which were shown to provide excellent fitting to all the stable shapes obtained by the dynamic simulations.Lavrenteva et al. (2021) used generalized Cassini ovals to approximate drop cross-sections and find stationary shapes by minimizing a convex functional of normal velocity at the interface.The analysis reproduced the branches with shapes of stationary stable flat drops and stationary toroidal ones available from numerical calculation in the case Bo = 0 and λ = 1.Furthermore, it predicts the point of loss of stability of the flat drop to exhibit the transition branch that leads into the formation of the toroidal shapes, and shows that this branch shows stationary, yet unstable, flat drops with ever-growing dimples up to collapse.However, as it was shown in Malik et al. (2020) that in the presence of rotation, the cross-sections of stationary drops did not always resemble Cassini ovals.Thus another method is required to find the branch of unstable dimpled shapes.
In this paper, we employ a method based on the use of a classical control theory, where the problem is reduced to a dynamical system with a low number of variables (states).We find a stationary solution of this system, linearize the system in the vicinity of this stationary solution, and on the basis of the linearized model, design a linear feedback controller.This controller is then used on the original nonlinear dynamic model of the system.For a control signal, we used the capillary number Ca that is proportional to the intensity of the outer flow.A similar approach was used in Malik et al. (2022) to stabilize toroidal drop shapes.Detailed descriptions of the algorithm and obtained results are presented in the forthcoming sections.

Simplified model for axisymmetric drop dynamics
Consider a two-parameter family of singly-connected bodies of rotation having the same volume as a unit sphere.Let the cross-section L of the surface of a body of this family be described by the parametric equations where θ ∈ [0, π/2], and the coefficients a k satisfy to ensure that the volume of the drop is equal to that of the unit sphere.
In what follows, we use the notation To construct the dynamic equations for these two parameters, we use the kinematic conditions at θ = 0 and π/2, i.e. at the points on the interface with coordinates (R max , 0) and (0, Z 0 ), and take R max (t) and Z 0 (t) for the states of the dynamic system.Due to the symmetry of the interface, at (R max , 0) we have u z = 0 and n z = 0, while at (0, Z 0 ) we have u r = 0 and n r = 0. Thus the dynamic system takes the form Here, u r and u z are computed from (2.10), with L defined by (3.1), where coefficients a 1 , a 2 and a 3 are given by Note that the right-hand sides of (3.4) are nonlinear functions of the model variables R max and Z 0 , and depend also on the flow characteristics defined by Bo and Ca.The latter will be utilized when the drop shape control is addressed in the sequel.

Feedback control to determine stationary solutions
A linear feedback controller is designed using an approximate model of the drop deformation presented in the previous subsection.Consequently, this controller is used for the nonlinear dynamic model formulated in § 2. Since the controller is robust only to limited variations in the model parameters, i.e. in the flow conditions, once the limit of flow variations for which a stationary solution can be found is reached, a new controller is designed for the new stationary solution that is close to that limit.The process is repeated and covers the entire range of flow conditions of interest.
Note that the dynamic model, even the simplified one, is nonlinear.However, since our main goal is to find stationary solutions, and not to perform advanced nonlinear regulation and tracking control tasks, we will employ a linearization-based classical linear control design, using the simple model and a single-input-single-output controller.
In the present study, for any fixed Bo, the control signal is chosen to be Ca.In a previous work (Malik et al. 2021) it was found that for a fixed Bo, R max dictates a unique drop shape.Therefore, the dictated control input and the feedback signal that is used by the controller are chosen as R max .
The control is designed on the basis of the approximate two-state model described in § 3.1.The states of this model are x(t) = (R max (t), Z 0 (t)), and their nonlinear dynamics is governed by (3.4).For the controller design, this model is linearized around a stationary solution denoted by x 0 , corresponding to Ca 0 .The linearized model is given by ) where ẋ denotes dx/dt, and The partial derivatives in the expression for A are calculated numerically using a second-order central difference scheme.
The open-loop transfer function between the control input δCa and the output δR max , which will be used in the controller design, is given by where s is the complex Laplace variable, and I is a 2 × 2 identity matrix.
The uncontrolled process corresponds to Ca = Ca 0 , i.e. δCa = 0.If at least one eigenvalue of A in (3.6), or equivalently, one of the poles of G ol (s), has a positive real part, then the stationary solution of (3.4) is unstable.In the current study, we chose R max to be the controlled output variable.
In order to design a controller that stabilizes the system of (3.6) and also can generate a neighbouring solution defined by δR D max = R D − R 0 max , we introduce the error term 983 A5-7 between the desired and calculated outputs: The control input δCa is chosen in a such a manner that the transfer function in (3.8) does not have a pole at the origin, and to ensure type one characteristics of the closed loop system, a proportional-integral controller is used to stabilize the system (Ogata 2010), given by δCa = −K p e − K I e I , (3.10a) where ėI = e. (3.10b) Here, K p and K I are proportional and integral gains designed to stabilize the closed loop system.The integral component of the controller is introduced for robustness and to better account for the modelling errors caused by the simplified model.Also, this component will allow us to move from one stationary solution to the neighbouring one.
In the present study, the parameters K p and K I are designed using the classical root-locus tuning techniques.The transfer function of the controller is given by The corresponding closed loop transfer function of the system is . (3.12) As discussed earlier, the above controller is used in the vicinity of the stationary solution x 0 to find neighbouring stationary shapes.This is carried out by applying this control logic to the full nonlinear model of the system (2.10).Since the linear controller computes only deviations in the capillary number, the value used in the simulations is given by Ca = Ca 0 + δCa, as is discussed next.

Control implementation to the full boundary integral formulation
For a chosen pair (Ca 0 , Bo), the controller gains K p and K I are determined as described in the previous section.Once the controller for the simplified model is established, it is applied to the full nonlinear dynamics model of (2.8) or (2.10) to determine the stationary solution with R max = R D .Namely, in (2.10), the capillary number is changed dynamically as For the initial conditions, we chose the established shape S 0 corresponding to particular (Ca 0 , Bo), and set some R D that is relatively close to R 0 max = max (r,z)∈S 0 (r).If the deviation in R max is sufficiently small and the controller is robust enough, then simulation will converge to a neighbouring solution of the nonlinear model.An example of the controlled dynamics is presented in figure 1 and Ca 0 = 0.1969, which is close to the critical point reported in Malik et al. (2020).The corresponding stationary shape has R max 1.70648.The gains were chosen as K p = 0.5473 and K I = 0.0274.Figures 1(a) and 1(b) show the dynamics of R max − R D and Ca, respectively.
In figure 1(a), the dynamics of R max − R D is presented.In this presentation, the desired stationary shape of the drop is assumed to be attained when |R max − R D | approaches zero while the capillary number tends to some constant value.The dynamics for R D = 1.706472 (solid blue line) begins with R max − R D = 0, and the modification in Ca(t) (see figure 1b) corrects the R max (t) value in order to approach the stationary solution corresponding to the desired shape.Evidently, the difference R max − R D increases initially, and after a certain time it starts decreasing towards R max − R D = 0.In the cases described by the dashed red, dash-dotted yellow and dotted purple lines in figure 1(a), the desired values of R max are obtained successfully in a similar manner.A previously obtained stationary shape is used as an initial shape for simulating the next desired shape.Note that the overshoot and the settling time of R max − R D are increasing with increasing R D values.The reason behind this may be the decreasing impact of a particular control on higher values of R D .In figure 1(b), in all cases, Ca(t) settles down quickly, and eventually stabilizes to the Ca values (i.e.Ca D ) that maintain R max = R D as depicted in figure 1(a).Note that all cases of the dynamics illustrated in figure 1 are obtained with the same control.
For non-zero values of the Bond number, an algorithm of receiving stationary solutions is similar: a control designed on the basis of the two-state model for a given pair (Bo, Ca 0 ), corresponding to a stationary shape with R 0 max , when applied to a full system with this very stationary shape as initial condition, converges to a stationary solution with R D ∈ (R 0 max ) i.e. in some vicinity of R 0 max .When the same initial stationary state fails to generate the new desired shape, we use the last obtained stationary shape of the stabilized drop as the initial shape for simulation, with the same controller gain values.As in the case Bo = 0, in cases where Bo / = 0, we begin with (Bo, Ca 0 ), close to the critical point report Malik et al. (2020).Note that the instability in drop shape develops faster with the advance along the unstable branch of the solution in any particular case of Bo.The stabilization of such 'more unstable', i.e. faster diverging, shapes may not be attained with the previously designed controller, K P and K I (this is a typical case for Bo / = 0), and is likely to require higher controller gains.At this point, we design a new controller for the updated stationary values corresponding to the new Ca and the corresponding shape.The simulation and stabilization of drops for neighbouring flow conditions are performed considering this new stationary shape.With the help of the procedure described above, we succeeded in obtaining stationary shapes close to those of collapse at the centre of the drop with Z 0 as small as 0.01.
Examples of such controlled dynamics for Bo = −5 are presented in figure 2. The dynamics for R D = 1.72 (dashed red lines) begins with the shape given in Malik et al. (2020) with R 0 max just a little lower than R D .Then the modification in Ca(t) (see figure 2b) corrects the R max (t) value in order to get the stationary solution corresponding to the desired shape.Evidently, after a short time, R max approaches R D .In the cases described by the dotted yellow lines, dash-dotted purple lines and solid magenta lines (figure 2a), the differences between the initial and desired values of R max are relatively high.For example, in the case described by the dotted yellow line (R D = 2.4), the control designed for R 0 max = 1.72 was employed.In this case, it took much longer time to achieve the stationary state, and one can see oscillations of R max (t) and Ca(t) at the initial period of time reflecting complex eigenvalues of the controlled system.Consequently, in the cases corresponding to the dash-dotted purple and solid magenta lines, new controllers were designed based on the updated stationary states and capillary numbers.This results in faster convergence and less noticeable oscillations.
The calculations described below were run up to t = 100, then the normal velocity at the interface was checked.The observed max(|u n |) at t = 100 was typically O(10 −5 ) or less, while in more complex geometry, such as deep dimples or in the vicinity of the region of transition from simply connected to toroidal branches, it was O(10 −3 ) or less.The particular levels of accuracy are indicated in the discussions following the results displayed below in § 5. Malik et al. (2020) reported on dynamic deformation of an immiscible singly connected drop in a linear flow that is composed of two simultaneous components, rotational flow and axisymmetric flow.It was demonstrated that if the Bond number and the capillary number lie within a certain domain in the phase plane, then the drops stabilize to stationary shapes.It has been claimed that outside these bounds, instabilities dominate the drop shapes; hence no stable stationary shapes were obtained.In the present work, we focus on finding new stationary unstable drop shapes.These are attained and stabilized with the help of a control algorithm as described in § § 3 and 4.

Stationary axisymmetrically rotating drops
Recall that the control approach starts with an approximation of drop shape close to the desired one that is deforming in the dynamically modified flow field to obtain the shape with a desired characteristic (R max = R D in our case).In other words, the control method aims to attain a desired drop shape, defined in terms of length of radial axis (R D ), by modifying the ambient flow field, represented by capillary number Ca.
The obtained stationary drop shapes of the cases depicted in figure 1(a) are shown in figure 3. Drop shapes in these cases of the outer flow corresponding to capillary numbers Ca = 0.196843, 0.195015, 0.189963, 0.183298 are shown as solid lines.Approximate solutions of similar cases obtained by Lavrenteva et al. (2021), using generalized Cassini ovals, are shown by dashed lines with Ca = 0.196, 0.195, 0.19, and 0.1831 for figures 3(a,b,c,d), respectively.The shapes are almost indistinguishable.However, the normal velocities at the interface approximated by the generalized Cassini ovals are of the order 10 −3 (see Lavrenteva et al. 2021), being several orders of magnitude higher than those associated with the shapes of the present paper in figure 3, i.e. max |u n | = 8.5 × 10 −6 , 1.5 × 10 −5 , 2.0 × 10 −5 , 3.3 × 10 −5 , respectively.
The computations of Malik et al. (2020) revealed that the heavy drops (Bo ≥ 0) lose stability with the appearance of a slight dimple at the central part of the drop.In figures 4(a,b,c,d), we present the controlled stationary dimpled shapes of heavy drops  for Bo = 4, and Ca = 0.021656, 0.012215, 0.001945 and 0.000003, respectively, with |u n | max ≤ O(10 −4 ).Note that the case presented by figure 4(d) resembles the shape of the drop for Bo = 4.0004 with Ca = 0 as presented in figure 3(d) of Malik et al. (2020).This regenerated result is represented by a dashed line overlapping the present result.Figures 3 and 4 suggest that shapes of rotating heavy drops are less extended in the r direction than those deformed by compressional flow without rotation.In both cases, with the decrease of Ca, an ever-growing dimple in the middle is evident, and at some critical Ca, when the thickness at r = 0 is totally collapsed, the drop may transform to a toroidal shape.
Unlike heavy drops, light drops (Bo < 0) lose stability at shapes without dimples (see Malik et al. 2020).Figure 5 presents the case Bo = −4, where this effect is demonstrated.Shapes depicted by dashed, dotted, dash-dotted and solid shapes correspond to (Ca, R max ) = (0.351324, 1.8), (0.326788, 2.2), (0.302938, 2.8) and (0.301842, 2.85), respectively.The stationary shape at R max = 1.8 is thick in the mid-section and is relatively thin at the radial endpoints.At lower Ca and higher radial extension, i.e.R max = 2.2, the sides of the drop rather than the mid-section are compressed and extended.This is due to the opposite effects of rotation and compressional flow.The thickness of the mid-portion decreases at R max = 2.8.In this case, most of the internal fluid is transferred to the radial endpoints of the drop.The drop achieved an almost collapsed shape when Ca = 0.301842 with R max = 2.85.This transformation described in the previous paragraph is typical for light drops that emerged in the fluids with a higher rate of rotation.A further example, at Bo = −5, is depicted in figure 6, where dashed, dotted, dash-dotted and solid shapes correspond to (Ca, R max ) = (0.3831, 1.72), (0.32815, 2.4), (0.32417, 2.8) and (0.31423, 3.19), respectively.The stationary shape at R max = 1.72 is thick in the mid-section and is relatively thin at the radial endpoints.At lower Ca and higher radial extension, the sides of the drop are more compressed than the mid-section, and a dimple in the centre appears only for highly compressed shape R max = 3.19, close to collapse.In the cases described in figures 5 and 6, |u n | max ≤ O(10 −3 ), slightly higher than those in figures 3 and 4.

Branching into unstable stationary flat drops
The obtained results on the stationary drop shapes are collected in figure 7, where the Taylor deformation factor D = (R max − Z max )/(R max + Z max ) is plotted versus capillary number for various values of the Bond number.There, the newly obtained controlled stationary shapes, beyond those obtained numerically by Malik et al. (2020), are represented by different coloured circles.
An interesting phenomenon, which is evident in figure 7 for all cases of Bo, is that there is a maximal value of the capillary number, Ca max (Bo), beyond which no stationary shapes were observed.These are marked by full diamonds.A further increase in R max is associated with a decrease in the capillary number toward the point of centre collapse at Ca = Ca collapse , and is expressed in the monotonous increase in D. The points next to collapse are marked by solid blue dots.Note that in the region Ca collapse < Ca < Ca max , there exist two different stationary solutions associated with different shapes, R max and D values, for a given Ca and Bo.The shapes with lower D are stable, while those with higher D are unstable and are stabilized by control.Note also that the extreme transition in shape of lighter drops, shown in figures 5 and 6, results in non-monotonous change of the curvature of deformation curves in figure 7.
A summary of this phenomenon is presented in figure 8(a), where Ca max (the critical region) and Ca collapse (the transition phase) are shown for varying Bo in the region −5 < Bo < 4. The upper curve corresponds to critical Ca where drops lose their axisymmetric stability.The lower curve corresponds to the total collapse of the centre region in the drop shape where transition to toroidal shapes is expected.The region between these two curves provides a subdomain of the (Bo, Ca) plane, where the singly connected stationary unstable shapes coexist with the stable stationary drops.Hence no solution exists above the upper curve, while only stable solutions prevail below the lower curve, as reported Also, in Malik et al. (2022) it was observed that there exists a thin domain on the phase plane (Bo, Ca) where two unstable toroidal shapes are possible.Now one can see that this domain lies inside the domain of double singly connected shapes.Thus under certain outer flow, up to four stationary shapes of the submerged drop coexist.This effect was observed  Malik et al. (2021Malik et al. ( , 2022)).Parts of the curves, above full circles, correspond to toroidal shapes that are stable with respect to axisymmetric disturbances.
previously for a rotation fluid mass (Fontelos et al. 2011).Also, approximate solutions of Lavrenteva et al. (2021) suggested that this is the case for the drop in compressional flow without rotation.
The shapes near collapse of the torus and near breakup of the singly connected drop at Ca collapse for the cases Bo = 4, 0 and −4 are depicted in figures 9(a), 9(b) and 9(c), respectively.The overlap of most of the corresponding surfaces is evident.Thus in figure 10, we combine the deformation plots of all four branches of stationarity: a stable branch of singly connected flat drops in the range Ca < Ca max , an unstable branch (stabilized by control) of highly deformed flat drops in the range Ca collapse < Ca < Ca max , and two unstable (stabilized by control) branches of toroidal drops in the range Ca < Ca collapse .We consider this figure as a major integral result of our efforts in Malik et al. (2020Malik et al. ( , 2021Malik et al. ( , 2022) ) and the current paper.

Concluding remarks
In this paper, we have addressed stable and unstable stationary branches of axisymmetric deformed flat drops.The unstable branches were rendered stable with the aid of a control scheme used by us in Malik et al. (2022) to stabilize unstable stationary toroidal drops.The linked description of branches given in figure 10 reveals a complex combination of possible deformation patterns of drops subject to rotation and compression or extension fields.It is noted that for any discussed Bond number, there appear to exist up to four stationary deformed shapes, with only one with the lowest deformation being naturally stable.The four branches are separated by Ca max and Ca collapse , being the bifurcation points, having values associated with the particular Bond number.We note here that this

Figure 1 .
Figure 1.Controlled dynamics of (a) R max − R D and (b) Ca, when Bo = 0.

Figure 8 .Figure 9 .
Figure 8. Lines marked by diamonds and circles correspond to the critical (loss of stability) and transition (collapse to toroidal shapes) points.

Figure 10 .
Figure 10.Deformation curves for various Bond numbers.Solid lines indicate singly connected shapes by Malik et al. (2020) and the present work; dashed lines indicate toroidal shapes byMalik et al. (2021Malik et al. ( , 2022)).Parts of the curves, above full circles, correspond to toroidal shapes that are stable with respect to axisymmetric disturbances.