Influence of Oil Viscosity on Hysteresis Effect in Electrowetting Displays Based on Simulation

1. Introduction

The optical switching in electrowetting display (EWD) technology is achieved through electrical control of colored oil dynamics in pixel cells [1], where contact angle modification is induced by voltage application at solid-liquid interfaces. This principle has been applied to various fields, including microfluidic manipulation in lab-on-chip systems [2], focal length adjustment in optical devices [3], and energy conversion through liquid motion control [4]. Among these implementations, EWD-based displays are currently considered to possess the most promising commercial viability [5]. EWD technology has the characteristics of simple structure, fast response, high reflectivity, high contrast, and rich colors [6]. However, due to the rough surface of the pixel structure and the viscous resistance of the oil–water interface, the hysteresis effect will appear during the pixel driving process [7,8]. As a result, the grayscales of EWDs cannot be accurately controlled. Numerous scholars have devoted substantial efforts to mitigate the hysteresis effect.

The hysteresis effect of EWDs refers to the difference in optical response during the driving process of pixel switching [9]. It was first discovered in 2006 by Van Dijk R of Liquavista when he improved an EWD to display video content with 4 bits of grayscales and a resolution of 170 ppi [10]. In 2012, a hysteresis-free pixel switching structure was designed by adding hydrophilic patches or staircases [11]. The hysteresis effect could be better eliminated by this method, but it would reduce the aperture ratio of pixels and bring unpredictable difficulties to the manufacturing process. In the same year, the electrowetting process of electrolyte drops on smooth and rough surfaces was studied by using the lattice Boltzmann model [12]. It was observed that the hysteresis and saturation of the contact angle were more obvious on the rough surface. In 2015, a theoretical model was proposed to describe the evolution of electrowetting on substrates with contact angle hysteresis, and the relationship between the apparent contact angle, voltage, and other parameters was quantified [13]. The results show that the theory and equation based on energy balance could successfully describe the electrowetting response with obvious contact angle hysteresis. In 2017, a method and principle for improving contact angle hysteresis was given based on the improved Young’s equation [14]. The relationship between lens focal length and applied voltage was studied by using three kinds of oil with different viscosity ratios as insulating liquid. The results show that the focal length hysteresis could be reduced by reducing the oil viscosity. In 2020, a pixel structure was proposed for achieving accurate control of oil rupture position and move direction by adding additional pinning structures or spacing arrays [15,16]. This scheme could effectively reduce the hysteresis effect of EWD.

To mitigate hysteresis effects through oil viscosity modulation, using the COMSOL Multiphysics simulation platform, the oil movement process in a pixel cell was successfully simulated by coupling the laminar two-phase flow phase field and electrostatic field, and the influence of oil viscosity on the hysteresis effect in EWDs was studied based on the proposed simulation model.

2. Principles of Hysteresis Effect in EWDs

2.1. Driving Principle of EWDs

An EWD panel is mainly composed of a glass substrate, indium tin oxides (ITOs), an electrode, a hydrophobic insulation layer (HIL), a pixel wall, colored oil, polar liquid, and a top plate [5,15], as shown in Figure 1. The gate control mechanism of the thin-film transistor (TFT) is realized by dual voltage regulation:

• Color display mode: When gate signal V_G1 is biased at a low level, TFT₁ enters the cutoff state, resulting in zero potential storage (C_S1 = 0 V). This passive configuration enables spontaneous oil spreading to achieve full-pixel coverage through interfacial tension dominance.
• Substrate exposure mode: Activation of V_G2 at a high level triggers TFT₂ conduction, enabling capacitive charging (C_S2 = V_S). The resultant electrostatic actuation induces oil film retraction to pixel corners, with contraction magnitude being voltage-dependently regulated for precise grayscale modulation.

Optical states correlate with oil contraction levels quantified by the aperture ratio ($W_A$), defined as the exposed substrate area proportion relative to total pixel area ($S_{pixel}$) [17], with computational expression provided in Equation (1).

$$W_A(V)=\left(1-{\displaystyle \frac{S_{oil}(V)}{S_{pixel}}}\right)\times 100\%$$

(1)

where $S_{oil}$ represents the oil-covered surface area in a pixel, $V$ represents the driving voltage applied to EWDs, and the area of pixel wall can be excluded from area calculations. Pixel walls are transparent grid structures that partition EWDs into discrete pixels.

2.2. Hysteresis Effect of EWDs

As shown in Figure 2, it is the relationship between driving voltage and aperture ratio when the conventional square wave is used to drive an EWD panel. The blue curve is the driving voltage rising stage, and the red curve is the driving voltage falling stage. It can be seen that there is a significant difference in the aperture ratio between the driving voltage rising stage and the driving voltage falling stage, that is, hysteresis effect. Generally, the absolute value of the aperture ratio difference under the same voltage value in the rising stage and the falling stage is called the hysteresis difference, and the aperture ratio curves formed in the rising stage and the falling stage are collectively called the hysteresis curves [18].

Combined with the morphological characteristics of the hysteresis effect, through theoretical analysis and literature review [19,20,21,22], it can be seen that the main causes of the hysteresis effect in EWDs are the existence of oil rupture voltage, oil contact angle hysteresis, and oil contact angle saturation. The oil rupture voltage is caused by the fact that the electric field force is not enough to overcome the static friction and capillary force when the driving voltage is too low. Oil contact angle hysteresis is due to the inconsistency between the oil forward contact angle and the oil backward contact angle caused by dynamic friction and viscous resistance. Oil contact angle saturation is caused by the weakening of the electric field force due to charge trapping.

3. Numerical Methodology

The numerical method was realized by establishing and solving numerical equations. The purpose was to track the morphological changes and motion characteristics of the oil–water interface in the limited pixel cell under the action of the external electric field force. In this paper, COMSOL Multiphysics was employed to simulate electrohydrodynamically modulated laminar biphasic flow through multiphysics coupling of hydrodynamic, interfacial, and electrostatic phenomena. Governing equations, encompassing the Cahn–Hilliard phase-field formulation, Laplace’s electrostatic equation, and Navier–Stokes momentum transport, were solved numerically using the finite element framework within a rigorously validated numerical scheme. First, the Maxwell stress tensor equation was solved by adding an electrostatic field module, and the electric field force output from the electrostatic field module was fed back to the hydrodynamic module. Then, the Navier–Stokes equation and phase field equation were solved by the hydrodynamics module according to the specified fluid boundary conditions and the solution results of the electrostatic field.

3.1. Governing Equations

There are two kinds of fluids in EWD pixel cells: water and oil. During calculations, the fluids are assumed to be an incompressible and incompatible Newtonian fluid, and there are no chemical reactions in the flow process. The Navier–Stokes equation with momentum conservation and the continuum equation with mass conservation are used to describe the laminar flow field [23]. The governing equations are as follows:

$$\rho \left({\displaystyle \frac{\partial u}{\partial t}}+u·\nabla u\right)=-\nabla p+\nabla ·\left(\mu \left(\nabla u+{\left(\nabla u\right)}^T\right)-{\displaystyle \frac{2}{3}}u\left(\nabla ·u\right)I\right)+F$$

(2)

$$F=F_{st}+\rho g+F_{vf}$$

(3)

where $u$ is the fluid velocity, $p$ is the fluid pressure, $\rho$ is the fluid density, $\mu$ is the hydrodynamic viscosity. And each item corresponds to inertial force, pressure, viscous force and external force acting on the fluid. The external force is mainly composed of surface tension, gravity and volume force. $F_{st}$ represents surface tension (N/m³), $\mathrm{g}$ represents gravitational acceleration (m/s²), $F_{vf}$ is the volume force (N/m³).

The fluid flow in the EWD system is laminar flow ($Re=5\times {10}^{-2}\ll 2000$), and there is no diffusion effect between the two fluids. Moreover, the gravity can be ignored ($Bo=1.4\times {10}^{-6}\ll 1$). Therefore, electrostatic field force is the main factor causing fluid flow [24]. The parameters requiring coupling are the volume force $F_{vf}$ and the surface tension $F_{st}$ caused by the applied electric field. The electrostatic volume force (N/m³) can be expressed by the divergence of the Maxwell stress tensor [24].

$$F_{vf}=\nabla T_{ik}$$

(4)

The Maxwell stress tensor can be expressed as Equation (5).

$$T=E{D}^T-{\displaystyle \frac{1}{2}}\left(D·EI\right)$$

(5)

where $I$ is the identity matrix, $E$ is the electric field, and $D$ is the potential shift field.

$$E=-\nabla V$$

(6)

$$D={\epsilon }_0{\epsilon }_{r}E$$

(7)

In the two-dimensional structural model, the Maxwell stress tensor can be expressed as Equation (8).

$$T=\left[\begin{array}{cc}{T}_{xx}& T_{xy}\\ T_{yx}& T_{yy}\end{array}\right]$$

(8)

The stress tensor in each direction can be obtained by substituting the parameters.

$$T=\left[\begin{array}{cc}{\epsilon }_0{\epsilon }_r{E}_x^2-{\displaystyle \frac{1}{2}}{\epsilon }_0{\epsilon }_r\left(E_x^2+E_y^2\right)& {\epsilon }_0{\epsilon }_r{E}_x{E}_y\\ {\epsilon }_0{\epsilon }_r{E}_x{E}_y& {\epsilon }_0{\epsilon }_r{E}_y^2-{\displaystyle \frac{1}{2}}{\epsilon }_0{\epsilon }_r\left(E_x^2+E_y^2\right)\end{array}\right]$$

(9)

The change in volume force caused by electrostatic field only occurs at the interface of two-phase fluid, and its derivative can be obtained.

$$F=\left[\begin{array}{cc}{\displaystyle \frac{\partial \left(T_{xx}\right)}{\partial x}}& {\displaystyle \frac{\partial \left(T_{xy}\right)}{\partial y}}\\ {\displaystyle \frac{\partial \left(T_{yx}\right)}{\partial x}}& {\displaystyle \frac{\partial \left(T_{yy}\right)}{\partial y}}\end{array}\right]$$

(10)

After an electric field is applied, the contact angle of the three-phase contact line will be changed by the Maxwell stress. The contact angle refers to the external angle during the movement of water on the hydrophobic insulation layer, which is expressed by the Lippmann–Young equation [9]. As shown in Equation (11), where $\theta$ is the contact angle when the driving voltage is $V$, and ${\theta }_0$ is the contact angle when the driving voltage is $0$, ${\epsilon }_0$ and ${\epsilon }_r$, respectively, represent the vacuum dielectric constant and the relative dielectric constant of the dielectric layer. And $d$ is the thickness of the dielectric layer, and ${\gamma }_{LG}$ is the liquid–gas contact line.

$$\cos \theta =\cos {\theta }_0+{\displaystyle \frac{{\epsilon }_0{\epsilon }_r}{2d{\gamma }_{LG}}}{V}^2$$

(11)

For two-phase fluid movement, the volume function method (VOF), level set method, or phase field method are usually used to track the oil–water interface [25,26]. Based on the diffusion interface model, the phase field method characterizes the interface evolution of the system through the changes in order parameters and concentration field. In the phase field theory framework of COMSOL, the comprehensive action of physical interface motion and chemical potential change is expressed by a differential equation so that the instantaneous change in interface is obtained. The phase field method does not directly track the change in the interface, but through the phase field variable $\varnothing$, it describes the fluid distance phase interface $\varnothing =0$. In the region of two fluids, the phase field variable is taken as $\pm 1$, while in the transition region of the interface, the phase field variable $\varnothing$ changes smoothly between $-1$ and $1$. The oil–water interface is defined as a diffusion interface, the change in the oil–water interface is described by Cahn–Hilliard convection equation, and the phase field variables are obtained by solving the phase field equation [27,28,29]. The phase field governing equations are as follows:

$${\displaystyle \frac{\partial \varnothing }{\partial t}}+u·\nabla \varnothing =\nabla ·{\displaystyle \frac{\gamma \lambda }{{\epsilon }^2}}\nabla \psi$$

(12)

$$\psi =-\nabla ·{\epsilon }^2\nabla \varnothing +\left({\varnothing }^2-1\right)\varnothing +({\displaystyle \frac{{\epsilon }^2}{\lambda }}){\displaystyle \frac{{\partial f}_{ext}}{\partial \varnothing }}$$

(13)

$$\sigma ={\displaystyle \frac{2\sqrt{2}\lambda }{3\epsilon }}$$

(14)

where $\lambda$ is the energy density (N), and $\epsilon$ is the capillary width varying with the interface thickness (m). These two parameters and the surface tension coefficient $\sigma$ (N/m) are described by Equation (14). $\gamma$ is the mobility parameter, and the relation function between $\gamma$ and $\epsilon$ is $\gamma =\chi {\epsilon }^2$, $\chi$ is the mobility adjustment parameter. The phase field variable $\varnothing$ is expressed as $1$ in oil and $-1$ in water. In order to correctly couple the electrostatic field and laminar flow field, it is necessary to track the dielectric constant, density and viscosity between diffusion interfaces. The calculation formulas are as follows:

$$\rho ={\rho }_1+\left({\rho }_2-{\rho }_1\right)\varnothing$$

(15)

$$\mu ={\mu }_1+\left({\mu }_2-{\mu }_1\right)\varnothing$$

(16)

$${\epsilon }_r={\epsilon }_{r1}+\left({\epsilon }_{r2}-{\epsilon }_{r1}\right)\varnothing$$

(17)

where $\rho$, $\mu$, and ${\epsilon }_r$ represent the density, viscosity, and dielectric constant of the fluid, respectively. The average curvature between two liquid interfaces is calculated as follows [29]:

$$\kappa =2(1+\varnothing )(1-\varnothing ){\displaystyle \frac{G}{\sigma }}$$

(18)

$$G=\lambda \left(-{\nabla }^2\varnothing +{\displaystyle \frac{\varnothing \left({\varnothing }^2-1\right)}{{\epsilon }^2}}\right)+{\displaystyle \frac{\partial f}{\partial \varnothing }}$$

(19)

where $\kappa$ and $G$ represent mean curvature and chemical potential energy, respectively.

3.2. Boundary Conditions

The boundary conditions were the precondition that the governing equations had a definite solution on the boundary of the region [30]. The size of the pixel cell used in the calculation domain is shown in Figure 3a, and the relevant boundary conditions are shown in Figure 3b.

In the EWD structure model, zero-charge boundary conditions were set around the model, that is, $n·D=0$. This also meant that no displacement field could penetrate the boundary at the surrounding inner boundary, and the potential was discontinuous at the boundary. For electrostatic field boundary conditions, the potential Vapp and grounding boundary conditions should also be specified. And the wetted wall, initial interface, inlet, and outlet need to be set for the phase field boundary conditions. The wetted wall boundary conditions were determined by Equations (20) and (21).

$$n··{\epsilon }^2\nabla \varnothing ={\epsilon }^2\cos {\theta }_w\left|\nabla \varnothing \right|$$

(20)

$$n·{\displaystyle \frac{\gamma \lambda }{{\epsilon }^2}}\nabla \psi =0$$

(21)

In the initial interface boundary condition, the two-phase flow contact interface was selected as the position of the initial interface, $\varnothing =0$. Both sides of the fluid of the EWDs were selected as the inlet and outlet boundary conditions, which were set as the outlet end in the actual simulation process, and the pressure constraint term $p=0$ was set at the end point. In addition, the initial value of volume force and velocity should be set, and the initial value of velocity was set to 0 (m/s).

Table 1. Material, geometric, and interfacial parameters of the proposed model [5,17,31].

Parameters	Quantity	Symbol	Value	Unit
Material	density of oil	${\rho }_{oil}$	880	kg/m3
	density of water	${\rho }_{water}$	998	kg/m3
	dynamic viscosity of oil	${\mu }_{oil}$	0.002	Pa·s
	dynamic viscosity of water	${\mu }_{water}$	0.001	Pa·s
	dielectric constant of oil	${\epsilon }_{oil}$	2.2	1
	dielectric constant of water	${\epsilon }_{water}$	80	1
	dielectric constant of hydrophobic insulating layer	${\epsilon }_{hil}$	1.95	1
	dielectric constant of pixel wall	${\epsilon }_{grid}$	3.28	1
Geometric	width of pixel cell	$w_{pixel}$	150	μm
	height of pixel wall	$d_{grid}$	5.6	μm
	width of pixel wall	$w_{grid}$	15	μm
	thickness of hydrophobic insulating layer	$d_{hil}$	1	μm
	thickness of oil	$d_{oil}$	5.6	μm
Interfacial	surface tension of oil and water	${\gamma }_{ow}$	0.015	N/m
	contact angle of pixel wall	${\theta }_{grid}$	65	deg
	contact angle of hydrophobic insulating layer	${\theta }_{hil}$	150	deg
	contact angle of top plate	${\theta }_{top}$	90	deg

Table 1. Material, geometric, and interfacial parameters of the proposed model [5,17,31].

Parameters	Quantity	Symbol	Value	Unit
Material	density of oil	${\rho }_{oil}$	880	kg/m3
	density of water	${\rho }_{water}$	998	kg/m3
	dynamic viscosity of oil	${\mu }_{oil}$	0.002	Pa·s
	dynamic viscosity of water	${\mu }_{water}$	0.001	Pa·s
	dielectric constant of oil	${\epsilon }_{oil}$	2.2	1
	dielectric constant of water	${\epsilon }_{water}$	80	1
	dielectric constant of hydrophobic insulating layer	${\epsilon }_{hil}$	1.95	1
	dielectric constant of pixel wall	${\epsilon }_{grid}$	3.28	1
Geometric	width of pixel cell	$w_{pixel}$	150	μm
	height of pixel wall	$d_{grid}$	5.6	μm
	width of pixel wall	$w_{grid}$	15	μm
	thickness of hydrophobic insulating layer	$d_{hil}$	1	μm
	thickness of oil	$d_{oil}$	5.6	μm
Interfacial	surface tension of oil and water	${\gamma }_{ow}$	0.015	N/m
	contact angle of pixel wall	${\theta }_{grid}$	65	deg
	contact angle of hydrophobic insulating layer	${\theta }_{hil}$	150	deg
	contact angle of top plate	${\theta }_{top}$	90	deg

4. Experimental Results and Discussion

The oil movement process in an EWD pixel cell was successfully simulated by the established simulation model. Therefore, based on this model, the influence of oil viscosity on the performance of EWDs can be investigated through experiments, including response time, maximum aperture ratio, and hysteresis effect. In experiments, the dynamic viscosity coefficient was used to characterize the oil viscosity, and the influence of oil viscosity on the hysteresis effect of the EWD was emphatically analyzed. All parameters were kept constant, and the temperature was kept at 25 °C. Segment function was used to describe the driving waveform, and the driving voltage ranged from 0 V to 30 V. The dynamic viscosity coefficient was dynamically changed, and the corresponding aperture ratio data were recorded.

According to Equation (1), the formula for calculating the aperture ratio in a two-dimensional model can be obtained. Assume that the pixel cell is square, and the side length is $a$. The contact surface between the oil and the hydrophobic insulating layer is round, with a diameter of $d$. Then, the calculation formula of the aperture ratio in the two-dimensional model can be obtained. The specific description is as follows, where $V_R$ is the rupture voltage:

$$W_A\left(V\right)=\left\{\begin{array}{c}0,V<V_R\\ {\displaystyle \frac{4a^2-\pi d^2}{4a^2}}\times 100\%,V\ge V_R\end{array}\right.$$

(22)

4.1. Simulation of Oil Movement Process

The oil movement process can be divided into an opening stage after applying the external voltage and a closing stage after removing the external voltage. A segment function was used to describe the traditional square wave, as shown in Figure 4a. Based on the established simulation model, the simulation of an entire oil movement process, including an opening stage and a closing stage, can be realized. As shown in Figure 4b,c, it was a simulation of the oil movement process in an EWD pixel cell. When $t=0 \mathrm{m}\mathrm{s}$, an external voltage of 30 V was applied, and the oil started to contract under the action of electric field force. At $t=0.5 \mathrm{m}\mathrm{s}$, the oil began to retract from the pixel wall, forming an exposed substrate region. At $t=2.3 \mathrm{m}\mathrm{s}$, the opening reached the maximum value. At this point, the oil opening process was completed. When $t=4 \mathrm{m}\mathrm{s}$, the applied external voltage was removed, and the oil began to spread under the effect of surface tension. When $t=5.4 \mathrm{m}\mathrm{s}$, the oil was closed and completely contacted the pixel wall. At this point, the oil closing process was completed.

4.2. Influence of Oil Viscosity on Response Time

In the experiment, a 30 V driving voltage was applied first and removed after the oil was completely opened. During this process, the pixel-on response time and the pixel-off response time were recorded when the oil was fully opened and fully closed under different oil dynamic viscosity coefficients. The relationship between the oil dynamic viscosity coefficient and pixel response time could be obtained, as shown in Figure 5. In the figure, the blue curve is the pixel-on response time curve, the red curve is the pixel-off response time curve, and the green curve indicates that the pixel could not be opened within 16 ms. It can be seen that the pixel-on response time and pixel-off response time increased with the increase in the oil dynamic viscosity coefficient. This was because the increase in oil dynamic viscosity would increase the viscosity resistance.

4.3. Influence of Oil Viscosity on Maximum Aperture Ratio

In the experiment, a 30 V driving voltage was applied. During this process, the maximum aperture ratio of a pixel cell under different oil dynamic viscosity coefficients was recorded when the oil was fully opened. The relationship between the oil dynamic viscosity coefficient and the maximum aperture ratio of a pixel cell could be obtained, as shown in Figure 6. The blue curve is the maximum aperture ratio curve of a pixel, and the red curve indicates that the pixel could not be opened within 16 ms. It can be seen that the maximum aperture ratio of a pixel cell was very close, indicating that the change in oil dynamic viscosity could not change the maximum aperture ratio of a pixel cell.

4.4. Influence of Oil Viscosity on Hysteresis Effect

The dynamic viscosity coefficient of oil commonly used in the preparation of EWDs is 2 mPa∙s. In the experiment, the traditional square wave was used to drive a pixel cell, and the hysteresis curves under different dynamic viscosity systems were recorded. The relationship between the oil dynamic viscosity coefficient and the hysteresis curve could be obtained, as shown in Figure 7a. In order to more intuitively reflect the relationship between the oil dynamic viscosity coefficient and the hysteresis effect, the data of aperture ratio in the experimental results were further processed. The absolute value of the difference in the aperture ratio under the same voltage during the rising and falling process was defined as the hysteresis difference. The relationship between the oil dynamic viscosity coefficient and the hysteresis difference curve could be obtained, as shown in Figure 7b.

It can be seen from Figure 7 that during the movement of the oil in a pixel cell, the hysteresis difference first increases and then decreases with the increase in the driving voltage, forming a maximum hysteresis difference in this process, indicating the degree of hysteresis effect. The maximum hysteresis value increases with the increase in oil dynamic viscosity coefficient and decreases with the decrease in the oil dynamic viscosity coefficient, indicating that the oil dynamic viscosity coefficient is positively correlated with the hysteresis effect. This is because the increase in oil viscosity will cause an increase in shear stress during oil movement, that is, an increase in viscosity resistance. Greater viscous resistance will consume more electric field force to do work, resulting in a change in the effective work carried out by the same electric field force during oil movement.

To demonstrate the effectiveness of the proposed model, we compared and analyzed the experimental results with the conclusions of relevant research work, as shown in Table 2. This study primarily investigates the influence of oil viscosity trends on the performance of EWDs, laying a critical foundation for future experimental validation.

While this work focuses on first-principles modeling of viscosity-dependent hysteresis, emerging techniques in human–robot interaction (e.g., adaptive LLM control) and hybrid sensing systems (e.g., multimodal strain analysis) offer promising avenues for closed-loop optimization of EWDs in practical environments [32]. For example, combining simulations with real-time optical feedback to mitigate hysteresis.

5. Conclusions

In summary, based on the COMSOL Multiphysics simulation environment, the driving process of the hysteresis effect in an electrowetting display pixel cell was successfully simulated by coupling the laminar two-phase flow phase field and electrostatic field. Based on the proposed simulation model, the influence of oil viscosity on the hysteresis effect of electrowetting display was studied. Our simulations demonstrated that the maximum hysteresis difference in the hysteresis effect increased with the increase in oil viscosity and decreased with the decrease in oil viscosity. The oil viscosity had little effect on the maximum aperture ratio of EWD. The pixel-on response time and pixel-off response time increased with the increase in oil viscosity. Beyond electrowetting displays, the findings of this study can also be extended to other voltage-driven microfluidic systems. Although the simulation framework offers valuable insights, we acknowledge the following limitations: simplified boundary conditions, temperature dependence, and scalability. These extensions and limitations will be further explored in our future work.