Numerical Demonstration of In-Tube Liquid-Column Migration Driven by Photoisomerization

Droplet manipulation by light-induced isomerization was numerically demonstrated and investigated regarding the driving mechanism. Such a non-invasive manipulation of a droplet in a microchannel can be realized, for example, by the use of watery solution of photoresponsive surfactant that exhibits the isomerization. Due to variable fluid properties between the cis and trans isomers, one-side light irradiation on a liquid column in a tube would lead to some kind of imbalance between the two ends of the liquid column and then drive droplet migration. The present numerical simulations of air–liquid two-phase flow and its scalar transport of the isomer, considering the variable static contact angle, agreed quantitatively with the experimental results in terms of the migration speed. This fact supports the contention that the droplet migration is more likely to be driven by an imbalance in the wettability, or the contact angle. The migration speed was found to be less dependent on the liquid-column length, but proportional to the tube diameter.


Introduction
Droplet manipulation techniques for microfluidic devices and lab-on-a-chip have attracted much attention in various fields such as medicine, chemistry, and biology [1][2][3][4][5]. This is because such an analytic operation on a micro scale increases the surface-to-volume ratio, thus obtaining many advantages, for instance, minituarization of samples, high-speed reaction, and downsizing of devices. Hence, several researches have devised a method of fluid driving in a micro-scale channel by changing the surface tension and wettability, e.g., a plasma-etched polymer nanostructure that enhances the droplet mobility [6], dielectrophoresis that employs an electric dipole moment by immersing an electrocode in a channel [7], and EWOD (Electrowetting on Dielectric) that generates the wettability gradient due to static electricity [8]. However, the challenges related to those methods include not only the fabrication of the channel but also the contamination by contact and the difficultly of flexible and selective manipulation. An alternative method using light as an external stimulus has several advantages such as not requiring the fabrication of a channel, simple adjustment of the stimulus, flexible and selective manipulation, and less contamination due to the non-contact [9]. In addition, it can be utilized on a living body due to its non-invasiveness.
There exist various ways of driving liquid by light irradiation: optical tweezers [10,11], photothermal Marangoni flow [12][13][14][15], and photoisomerization [16][17][18][19][20][21][22][23]. Optical tweezers is a method that uses the radiation pressure of light, but its corresponding force is very small, in the order of pico-newtons. The use of photothermal Marangoni flows should be accompanied by heat, together with the light irradiation. This method is flexible and selectively operable, and the resulting force is sufficient for manipulation; however, it cannot be applied to a substance that is denatured by heat. On the other hand, the fluid manipulation by utilizing cis-trans photoisomerization can avoid the use of heat, although it is necessary to change the fluid properties in response to light.
As shown in Figure 1, the cis-trans photoisomerization is a property in which the cis and trans isomers are reversibly changed by light of a specific wavelength such as ultraviolet (UV) light [24]. The isomers are represented by the same molecular formula but different molecular structures. As the molecular structure changes, fluid properties such as the contact angle and surface tension are varied. Here, the reaction rate constant k represents a rate at which the reactants increase or decrease in a chemical reaction. In the photoisomerization, k represents the degree of photo-induced change between the cis isomer and the trans isomer. Azobenzene is a representative example [25][26][27][28]. In previous studies, droplets on a substrate were manipulated by performing photoisomerization on the substrate and changing wettability [16][17][18][19][20], or on photoresponsive surfactant [20][21][22][23]. Recently, Muto et al. [23] demonstrated the manipulation of a liquid column in a millimeter-scale glass tube. Its droplet manipulation was done by UV-light irradiation on a one-side surface of the liquid column generating differences in the contact angle and surface tension between both sides of the finite liquid column. The mobility of the migrated liquid column was reported to depend on the liquid column length. However, such an experimental demonstration might often suffer from a pinning effect that cannot be avoided, and its flow dynamics and developments of each isomer distribution are not fully understood, thus increasing the difficulty of experimental measurement. To the authors' knowledge, no numerical simulation of such a liquid driving has been performed.  [23].
The cis isomer irradiated with visible light should change into the trans isomer, while the trans isomer irradiated with ultraviolet (UV) light would change into the cis isomer. These two isomers may exhibit different fluid properties. The rate at which the number of isomers increases or decreases by photo-induced isomerization is expressed as a reaction rate constant k.
In order to reveal the light-induced droplet-migration phenomenon due to cis-trans isomerization in a watery solution of photoresponsive surfactant, we carried out a numerical study of in-tube liquid columns, considering the variable contact angle and surface tension. We used a framework of OpenFOAM R (version 2.3.1) [29], which is an open source software and has been verified by benchmark tests for multiphase flows [30][31][32][33][34]. The Volume-of-Fluid (VoF) method [35] is used as an interface-capturing method, and the Continuum-Surface-Force (CSF) model [36] is applied to calculate the surface tension on liquid-air interfaces. In addition, we have implemented the Continuous-Species-Transfer (CST) method [37] to express the cis/trans-isomer transports more accurately. We discuss a validation with the experiment [23], and numerically investigate the isomer distribution and the effects of the liquid column length and radius on the liquid-column migration.

Problem Setting: A Photoisomerizable Liquid Column in a Tube
We focus on an experimental demonstration performed by Muto et al. [23] and employ a similar problem setting for our numerical analysis. Figure 2 shows the present analysis object, which is a liquid column given in an infinite straight tube with a constant radius R. The liquid column was initially placed at the center of the computational domain. A UV light is assumed to be irradiated on the right half of the domain: see Figure 2. The irradiation is started from the state of all trans isomers in the liquid of interest; that is, the initial C cis was 0 in the entire domain. Since the tube is on the millimeter-scale, as tested by Muto et al. [23], the UV light is assumed not to decay throughout the liquid [38]. In the simulation, we used a wedge mesh which consisted of a single grid cell in the circumferential direction by assuming an axisymmetric flow with respect to the z axis. This allows us to reduce the computational cost by keeping fine meshes in the other directions. While Muto et al. [23] used a sufficiently-long open tube for their measurement, it is practically difficult to simulate such a system rigorously. Two different boundary conditions in the z direction were tested in this study: the periodic boundary condition (PBC) and the inlet/outlet boundary condition. On the tube surface, the no-slip boundary condition was applied. The initial length of the gas phase on both sides of the liquid column was set at L g = 15 mm, regardless of the liquid column length L c . The liquid column is initially located at the center of the tube domain and in the form of a simple column shape, being surrounded by air. The length of the liquid column is denoted as L c , while that of the gas phase is L g . Table 1 shows the fluid properties of our present targets. Since we simulated both the air and liquid phases simultaneously, the air properties at room temperature were given. Essentially important properties in the present problem are the contact angle ϑ and the surface tension σ, as explained in Section 1. Those values vary depending on the isomer, and the liquid of our interest reveals rather hydrophilic features with the cis-isomer: the contact angle of purely cis-isomer liquid, ϑ cis , is slightly lower than that for the trans-isomer, ϑ trans . We referred to the values of ϑ cis and ϑ trans measured experimentally by the extension/contraction method, and σ cis and σ trans measured by the pendant drop method [39]. The reaction rate constant k for the photoisomerization was identified by 1 H-NMR measurement. The concentration diffusivity D for the isomer diffusion in each fluid was chosen as a typical value: cf., Ref. [28].

Governing Equations for Fluid Motions
Although the actual system of our interest consists of incompressible liquid and compressible air contained in a tube, we considered the air phase as incompressible for simulating the fluid behavior. This assumption allows us to use the governing equations of incompressible and immiscible gas-liquid two-phase flows: the equation of continuity and the Navier-Stokes equation which includes a surface tension force f σ that works on the liquid-air interfaces based on the Continuum-Surface-Force (CSF) model [36]. In the present study, all simulations were under the zero-gravity condition. We used the Volume-of-Fluid (VoF) method [35], which is a well-known way to capture an interface between two different fluids. In this method, the advection equation of the VoF function α can be written as The VoF function α describes the liquid fraction in each computational grid cell to determine the two-fluid allocation. In the present air-liquid two-phase flow, we defined it as : Liquid phase 0 : Air phase (0, 1) : Interface (4) Instead of using simply Equation (3), one may compute the following equation with an artificial term with the aim of sharpening the liquid-air interface.
The third term on the left hand side of Equation (5) represents the artificial interface compression term, in which the compression velocity u r is calculated, as follows: Here, n f is the normal vector of the interface in a grid-cell control volume, φ is the mass flux, S f is the area of the interface in the control volume, and C α is an adjustable coefficient. While C α can be changed arbitrary, we decided to set C α = 0 in this study to make unavoidable spurious currents less pronounced, as discussed in Appendix A.1.
The local density ρ and viscosity µ at each grid cell depend on α: The surface tension term in Equation (2) is expressed as (9) in the CSF model [36]. The unit normal vector n and the curvature κ of a local interface can be expressed by the VoF function α, respectively, as Note here that Equation (9) includes only the normal force on the liquid-air interface, neglecting the tangential force. The normal force induces the Laplace pressure, while the Marangoni convection may be triggered by tangential force on the interface. Compared to the normal force, the expected tangential force would be much smaller in this study, because the difference between σ cis and σ trans is almost negligible relative to its absolute value, as given in Table 1.

Representation of Cis-/Trans-Isomer Liquid
Since the driving force exerted on a liquid column by the cis-trans photoisomerization is induced by the imbalance of the surface tension and contact angle (or Laplace pressure), it is necessary to express numerically the two different states of either the cis or trans isomer and to model the exchange between them due to the photoisomerization. Then, let us introduce C cis , which is defined in each computational grid cell as the volume-fraction ratio of the cis isomer: : All cis isomer 0 : All trans isomer (0, 1) : Mix of cis & trans isomers (11) Both the contact angle and surface tension were determined as a function of C cis : A validation test was performed on a simple droplet on the flat wall with a variable contact angle that depends on C cis , and we confirmed a reasonable C cis -dependence of the contact angle as well as a response of the droplet wetting to UV-light irradiation, as presented in Appendix A.2.

Transport of the Cis/Trans Isomer with the CST Method
The transport equation of the isomer ratio C cis can be written as follows: where R is the source term due to the photochemical reaction and D is the arithmetic mean diffusivity: In the actual phenomenon, no transport of the isomer across the liquid-air interface should occur. However, unavoidable spurious currents, i.e., numerically-artificial flows, in simulation may lead to scalar transport across the interface. To suppress such a leakage of scalar (or C cis in this study) across the liquid-air interface, we implemented the Continuous-Species-Transfer (CST) method [37]. Then, Equation (14) is reformulated into where Φ is the discontinuity term and it works at the interface, as follows: The degree of suppression can be calibrated by the Henry coefficient H. We set H = 100: see also Appendix A.1.
It should be noted that, in the present simulation, only the change from the trans to cis photoisomerization by UV-light irradiation was expressed by but the opposite change from the cis to trans photoisomerization (by visible light) was not considered. Because the latter photoisomerization would not provide any contribution to the liquid-droplet migration under the present condition. Here, k in Equation (18) is the reaction rate constant of change from trans to cis, and δ UV is Equation (18) expresses that the UV-photoisomerization occurs only in a liquid portion that is irradiated with UV light, by multiplying α and δ UV . Basically, C cis = 0 should be kept in the air and/or non-irradiated part. Figure 3 shows our obtained migration distance z m , normalized by the liquid column length L c , from the initial position as a function of the irradiation time T. The instance of T = 0 corresponds to the beginning of the simulation as well as that of the UV-light irradiation. The migration distance z m from the initial position z 0 was calculated from the axial shift of the center of gravity of the liquid column.

Grid Resolution: Comparative Validation with Experimental Results
According to the experimental results shown also in Figure 3, the liquid column would be transported at a constant speed toward the UV light side (i.e., to the positive z direction). As z m /L c approaches 0.5, the migration speed appears to decelerate and finally the liquid column stops moving. Note here that z m /L c = 0.5 corresponds to the situation that the entire liquid column has just entered the UV area. To demonstrate such phenomena in our simulation, we first examined the required grid resolution, or the number of grids in the axial-(N z ) and radial-(N r ) directions. Our parametric study for N z = 500-2000 and N r = 25-300 revealed that a combination of (N z , N r ) = (512, 256) was required for the case of L c = 20 mm, as shown in Figure 3a. The linear motion in the initial stage could be reproduced even with lower resolutions, but the migration speed was overestimated remarkably. Furthermore, regarding the z-direction boundary condition, the periodic boundary condition (PBC) was found to provide a better result, compared to the inlet/outlet condition. As for a longer liquid column of L c = 30 mm, we achieved good quantitative agreement with the experimental results, as shown in Figure 3b. Similar to the shorter case mentioned above, the liquid column exhibited an acceleration period until T < 20 s, and it attained a constant migration speed at T ≈ 20 s. Whereas the linear motion obviously terminated before reaching z m /L c = 0.5 in the experiment, the present simulation shows a continuous linear motion until z m /L c reaches 0.5. Such a deviation in this late stage of the droplet migration must be attributed to the pinning effect that cannot be simulated with the present code. At least, in the initial stage including the linear motion, the photoisomerization-induced droplet migration has been reasonably demonstrated by our simulation. Figure 4 shows the developing distribution of C cis in the z-r section of the domain including both the air and liquid phases. The top panel in the figure represents the initial condition. At T = 0.01 s, the liquid-air interface is already concave due to the wettability with ϑ < 90 • and surface tension. This implies that the timescale of surface deformation is much faster than the reaction rate of the photoisomerization of interest. From T = 20 s, one can recognize that the isomers in the UV-light irradiation side changed gradually from trans to cis with time, but it can also be recognized that an axial shift of the liquid column already occurred. Figure 5 shows the temporal variation of the distribution of C cis on the z-axis inside the liquid column, which is marked with a white dashed line in Figure 4. The horizontal axis represents the z subtracted by the moving distance z m of the liquid column: z − z m = 0 corresponds to the center of the liquid column and the air-liquid interfaces locate at about z − z m = ±10 mm because L c = 20 mm. As seen in the figure, C cis increases firstly on the UV-light irradiation side because of the photoisomerization. As time progresses, the plateau of the C cis profile becomes high and narrow (z − z m = 2-8 → 6-7 as T = 20 → 80), while the values gradually start to increase also in the column that is half on the non-irradiated side. This is the result of competition between the photoisomerization reaction and the incoming trans-isomer liquid from the non-irradiated side. At T = 80 s, the non-zero C cis region reaches the interface on the non-irradiated side. This is to be expected due to the fact that the liquid column starts to decelerate its migration at about T = 80 s. From this, we may draw the conclusion that the liquid column would stop when the cis isomer concentrations at both sides of the liquid column are comparable.    Figure 4 for L c = 20 mm. The horizontal axis represents the relative position with respect to the center of the liquid column. The interface on the UV-light irradiation side is located at z − z m ≈ 10 mm, and the opposite one is at z − z m ≈ −10 mm.

Influence of the Liquid Column Length and Tube Radius
Furthermore, the dependence of the migration motion on the column length L c was investigated. Four column lengths of 5, 10, 20, and 30 mm were examined under the same system with a fixed tube radius R = 1.25 mm, and the initial exposed column length was set to L c /2 for all cases. Figure 6 compares the four cases, where the profiles are plotted in a dimensional form of millimeter versus second. From the gradient of each profile, the migration speed can be determined. The present numerical results obviously reveal the independence of the migration speed on the column length, at least in the initial stage of the migration. This aspect contradicts the experimental observation [23], but this mismatch might be due to a pinning effect on an actual tube in the experiment. With the increasing L c , the migration distance is elongated further. The constant linear motion is found to be at a speed of ≈ 0.15 mm/s, irrespective of L c . This speed could be varied by changing the fluid properties as well as the tube geometry. Further investigation on the fluid-property dependence of the migration speed is potentially interesting for a practical application of the photoisomerization in a droplet manipulation. The dependence on the tube radius was not investigated experimentally [23], so we numerically examined it under the same fluid conditions but with various R. Figure 7 shows results for R = 0.25-2.5 mm. The temporal evolution of the migration distance did not change qualitatively with a different R, but the migration speed was found to vary significantly. Normalized by each tube radius, the profiles are scaled well, as given in Figure 7b. According to the present numerical result, we conjecture that the migration speed is roughly proportional to the tube radius. This aspect is consistent with Muto et al. [23] with respect to their discussion on the net driving force F that was exerted on the liquid column. They suggested a form of where ϑ a-cis and ϑ r-trans denote the advancing contact angle of the cis isomer and the receding one of the trans isomer, respectively. Actually, Equation (21) does not relate to the column length L c , but includes the tube radius R, supporting also the L c -independency observed in Figure 6. Last, but not least, we would argue that the main driving force of the liquid-column migration by photoisomerization is attributed to the difference between the static contact angles ϑ cis and ϑ trans . In the context of Equation (21), the given condition of ϑ cis < ϑ trans reasonably results in F > 0. In terms of the surface tension, σ cis that is slightly larger than σ trans might also contribute to the onset of a positive F, but a comparison such as cos ϑ cis / cos ϑ trans > σ cis /σ trans (cf. Table 1) implies that the contact-angle variation would account for the imbalance between the two ends of the liquid column. On the other hand, Muto et al. [23] declared that the migration was induced by a difference in the surface tension rather than the contact angle, because they observed ϑ a-cis ϑ r-trans ; that is, cos ϑ a-cis < cos ϑ r-trans and this fact should demand σ cis σ trans to derive the positive F from Equation (21). However, they still suffered from difficulties in measuring the dynamics contact angle and surface tension and did not perform quantitative evaluation. Although our conclusion opposes their hypothesis, the quantitative agreement in the migration speed with the experimental result supports the present conclusion that the wettability, or the contact angle, is key to the driving force for the present droplet manipulation.

Conclusions
We performed numerical simulations to investigate the liquid-column migration driven by the cis-trans photoisomerization phenomenon that was experimentally demonstrated by Muto et al. [23].
The liquid-column migration should be induced by an imbalance in the wettability, or the contact angle and surface tension, between the two ends of the liquid column. To track the liquid-air interface, we employed the VoF method in conjunction with the CST model and the CSF model. In order to express the developing distribution of the cis/trans isomer, we also defined a volume-fraction ratio of the cis isomer, on which both the contact angle and surface tension were dependent. Neglecting the gravity, evaporation, thermocapillary, and dynamic contact angle, our simulation successfully demonstrated the photoisomerization-induced droplet migration and achieved good agreement with the experimental results. Our conclusion regarding the mechanism of the present droplet manipulation is that the driving force is caused mainly by the imbalance in the wettability, or the contact angle, between the two ends of liquid column rather than the surface tension. Through a numerical investigation of the cis isomer distribution, which is difficult to measure experimentally, we confirmed that the liquid-column migration terminated when the cis isomer distribution reached the non-irradiated region. We also found that the migration speed was less dependent on the liquid-column length and was proportional to the tube diameter.
Author Contributions: K.N. and T.T. conceived and designed the simulations; K.N. performed the numerical simulations and analyzed the data; K.N. and T.T. wrote the paper. All authors discussed the results and commented on the manuscript.
Funding: This research received no external funding.

Acknowledgments:
We are particularly grateful to Masahiro Motosuke, Ken Yamamoto, and Masakazu Muto, at Tokyo University of Science, for their own experimental data and valuable discussions. Some simulations were performed using SX-ACE at Cybermedia Center, Osaka University. We also thank the reviewers for their insightful and constructive comments.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Preliminary Tests
The present simulations were executed in a framework of a customized opensourse solver in OpenFOAM R [29]. By tuning the solver code, we combined several schemes of the CST model and the CSF model to obtain sharp liquid-air interfaces and to suppress scalar mass flux across the interfaces. For those models, we should determine the optimal values of the relevant numerical parameters. In particular, the Henry coefficient H in Equation (17) and C α in Equation (6) are crucial parameters. As a preliminary test to examine H and C α , which provide a reasonable result for the presently-given combination of fluids and scales, numerical simulations were carried out on a square droplet. In addition, the solver code was extended to simulate the photoisomerization in a part irradiated with UV light. To validate our code in terms of the modelling of the photoisomerization, we simulated a droplet on a flat plate, where the relationship between the wetting width and the contact angle is known theoretically. Figure A1 schematically illustrates the configurations of two preliminary simulations, the results of which will be reported in the following subsections.  Figure A1a. The computational domain size and the number of grids were L x = L y = 5 mm and N x = N y = 100, respectively. The surface tension was fixed at σ trans = 40 mN/m, since no UV-light irradiation was considered here. Other physical properties were the same as the liquid column driven by the photoisomerization, as given in Table 1. By systematically changing the Henry coefficient to H = 0.01, 1, and 100, and the interface compression value to C α = 0, 0.5, and 1, we investigated the degree of leakage of C cis through the liquid-air interface and the influence of spurious currents occurring near the interface. Figure A2 shows the cross-sectional views of the droplet, where the interface has been deformed into the circle shape by the surface tension. In the simulation without the CST model, we detected a severe leakage of the cis isomer from the droplet, as shown in (b). Figure A2c also shows a significant decrease in C cis inside the droplet. By applying the CST model with a large Henry coefficient H = 100, we successfully avoid an increase in C cis outside the droplet, as visualized in (d)-(f). Non-zero C α values may provide more suppression of variations in C cis outside of the droplet, but the C cis conservation inside the droplet would obviously suffer from the spurious currents that exist near the liquid-air interface, as shown in (e) and (f).
In order to confirm the suppression of leakage from the droplet by employing the CST model with C α = 0 and H = 100, we define an index representing the conservation of C cis in the liquid.
The following index C cr is a ratio between an arbitrary time T and the initial t = 0, regarding the volume-integrated cis-isomer concentration: Here, A represents the area of the computational domain, since the present simulation was performed in two-dimensional space. This ratio C cr would be time-dependent in simulations, once a leakage of C cis from the liquid to air phase occurs. In the ideal situation, C cr should be 1 at any time. As shown in Figure A3, the case visualized also in Figure A2d indeed exhibits less deviation from the unity value. According to this result, we have chosen the CST model with C α = 0 and H = 100 for the main liquid-column simulation.

Appendix A.2. Changed Static Contact Angle by Photoisomerization
An initially-hemispherical droplet with a radius of R 0 = 0.5 mm was placed on the wall, as shown in Figure A1b. The computational domain size and the number of grids were L x × L y = 3 × 1 mm 2 and N x × N y = 150 × 50, respectively. The surface tension was fixed at σ trans = 40 mN/m without any change due to photoisomerization, but the static contract angle ranged from ϑ trans = 90 • to ϑ cis = 30 • . The reaction rate constant in this droplet simulation was set at a rather large value of k = 1 s −1 . Other physical properties were the same as the main liquid-column simulation, as shown in Table 1. Figure A4 shows the cross-sectional views of the droplet, where the contact angle between the liquid-air interface and the bottom solid wall can be observed to change gradually from the given ϑ trans to ϑ cis . Such a decrease in the contact angle can be interpreted as an increase in the wettability of the liquid. x (mm) Figure A4. Temporal evolution of a droplet due to the photoisomerization-induced contact-angle change. Contour shows C cis . The contact angles for C cis = 0 and 1 are 30 • and 90 • , respectively. Solid cyan line outlines the droplet shape, and the purple-colored regions represent the UV-light irradiation area.
The droplet wet length L can be derived theoretically from the initial droplet size R 0 and the contact angle ϑ. By assuming that the inertia should be negligible compared to the surface tension and viscosity, the expression of L = 2R 0 sin ϑ π 2 (ϑ − sin ϑ cos ϑ) (A2) would be valid. In the present (preliminary) system, ϑ should be changed monotonically from that for the trans isomer to that for the cis isomer as a function of C cis , which is also time dependent in the UV-irradiated area and can be estimated mathematically on the basis of Equation (18). In Figure A5, the predicted evolution of the wet length is plotted with a dashed line. At T = 1 s, C cis in the light-irradiated liquid has attained C cis = 1; that is, a state of all cis isomers in the vicinity of the wall. After that, no expansion of the wet length occurs and the droplet maintains the contact angle of ϑ = ϑ cis . Also shown in Figure A5 is the result obtained by the present simulation. It can be confirmed that the present simulation reasonably demonstrates the temporal variation of the wet length, implying the validity of our simulation code with respect to the photoisomerization and its accompanying variation of the wettability.