Adjoint-based optimization of the actuator velocity profile in an inkjet print head

We consider the thermo-viscous acoustic flow inside the narrow channel and nozzle of an inkjet print head. We define a cost function to be the sum of the acoustic energy in the channel and the surface energy of the spherical cap of ink at the end of the nozzle. We derive the adjoint equations and obtain the sensitivity of this cost function to boundary forcing from the piezo-electric actuator opposite the nozzle. We use this forcing to eliminate residual oscillations after a droplet is ejected. We use a gradient-based optimization algorithm to find the time-varying boundary forcing that minimizes the cost function at various final times and for geometries with increasing complexity. For all geometries, the actuator first extracts fluid so that the ink/air interface becomes flat. This unavoidably sends an acoustic wave upstream, which reflects off the inlet manifold. The actuator subsequently moves to absorb this returning wave without reflection. The optimal boundary forcing and the final energy depend on the channel length, the actuator length, the forcing’s temporal resolution, and the available optimization time. The minimum time required to dampen residual oscillations is the time taken for waves to travel from the actuator to the inlet and back. For times greater than this, the total energy inside the microchannel can be reduced by a factor of 1000 compared to the uncontrolled case. This method is general and can be applied to other cost functions and initial conditions. Successful application of this method could lead more repeatable droplets at higher ejection frequencies.


Introduction
An inkjet print head contains several hundred microchannels, each typically 1 mm in length and 70 µm in width and height (Fig. 1).Each channel feeds a short nozzle, typically 50 µm in length and 20 µm in diameter.There is a free surface at the end of each nozzle, where the ink meets the air.A droplet is expelled through a nozzle when the piezo-electric actuator in the corresponding channel pushes in one side of the channel by a few nanometres.After this, acoustic oscillations reverberate through the channel until damped by visco-thermal mechanisms, which strongly influence wave propagation inside narrow channels [1,2].Ref. [3] contains a detailed review of the process, the challenges, and the open questions in inkjet printing.
Drop-on-demand inkjet printing is widely used for printing paper, packaging, and textiles.
It is increasingly used for production of 3D electronic components, MEMS devices, and other applications [4,5].Manufacturers would like to increase the droplet ejection frequency, while retaining or improving the reproducibility of the droplets.There is, however, a trade-off between the droplet ejection frequency and the droplet reproducibility.This is because, as the time between ejections decreases, each droplet becomes increasingly affected by the residual oscillations from the previous ejection.Manufacturers therefore try to design print head shapes that passively damp residual oscillations.They also alter the electrical waveforms sent to the print heads microchannels in order to damp residual oscillations with open loop control.These waveforms are currently adjusted by trial and error in extensive experimental campaigns.Our first paper on this subject [6] used adjoint methods in the frequency domain.That paper showed how to optimize the shape of the print head microchannel in order to passively damp residual oscillations.This paper uses adjoint methods in the time domain in order to damp residual oscillations with open loop control by optimizing the actuating waveform: the velocity at the actuator channel wall.The adjoint methods in both papers reveal the physical mechanisms being exploited by the optimization algorithms.
Early studies [7,8] performed experimental and theoretical analysis of acoustic phenomena inside a drop-on-demand piezoelectric print head.Acoustic motion inside tubes or channels with simple geometries can be successfully approximated using reduced analytical models, such as the narrow channel model [1,2,9] and the lumped elements method [10].For channels with smooth boundaries, the acoustic boundary layer model [11] can be applied.While these models are computationally cheap and useful for initial designs, they can fail to predict performance of print heads with complex geometries.In this paper we impose no restrictions on the channel geometry and solve the full system of thermo-viscous acoustic equations in two dimensions.The method can easily be extended to three dimensions.
The flow in the nozzle, which is short, is typically modelled as incompressible and axisymmetric [12,13], driven by the pressure [14] or velocity inflow [15] from the channel.Separating the nozzle from the channel allows the two domains to be considered separately and, if analysed numerically, to be discretized independently.In turn this allows the analysis in the nozzle to be focused on droplet formation.Various numerical methods have been proposed to model the free surface development, droplet formation, and pinch-off in inkjet print heads.One-dimensional models of droplet formation [16,17,18] and jet break-up [19] have shown good agreement with experiments.The volume of fluid [13] and level-set methods [20,15] have been successfully applied to compute the position of free surface boundary and droplet formation.In the above methods, the nozzle is assumed always to be filled with fluid.This assumption can be relaxed and the dynamics of the contact line between liquid, solid, and gas can be simulated using a variational approach [14,21].In this paper, however, we assume that the nozzle is always filled with fluid and that the fluid/gas interface is a spherical cap.This is a realistic assumption, and it has the advantage that it allows gradients of useful quantities with respect to the inflow conditions to be derived analytically [22,23].It could readily be relaxed in subsequent studies, at the expense of extra computational time.
Typical objectives for optimization in inkjet printing are droplet velocity and volume [24], damping of residual oscillations, damping of cross-talk between channels [25,26], and high frequency jetting [15].Reducing the nozzle diameter allows jetting at megahertz frequency and smaller droplet size [15].We do not calculate the waveform of the applied voltage because this would require us to specify and then simulate the piezo-electric actuator, which varies for different devices.We choose instead to focus on the optimization and physical understanding of the fluid motion alone, in order to produce a general result for the motion of the fluid boundary.
With this knowledge, the required voltage can then be found for a given actuator.
The most basic waveform is the unipolar (trapezoidal push in) waveform.This waveform has a number of drawbacks: the volume of the ejected droplets is large, the amplitude of residual reverberations is high, and satellite droplets are formed in addition to the main droplet [27].
Bipolar (trapezoidal push in then trapezoidal pull out) waveforms are used to eliminate satellite droplets [24].Compared with unipolar and bipolar waveforms [12], the W-shaped waveform type can significantly reduce the volume of the ejected droplet, and eliminate the residual acoustic waves from the last ejection cycle [24].Droplets formed from complex waveforms are, however, more sensitive to changes in the waveform shape.This means that optimal waveforms become harder to find as the waveform type becomes more complex.Given that waveforms are usually found by trial and error during extensive experiments, this leads to considerable experimental cost as the waveforms become more complex and motivates the more systematic approach in this paper.
Systematic waveform optimization can be approached in several ways.A feed-forward control method [25] can be used to eliminate residual reverberations by flattening the response of the meniscus velocity to the pulse frequency.If numerical models are not accurate enough, or are too computationally expensive to predict the droplet characteristics, then the waveform parameter space can be explored with model-free methods by combining an automated experimental rig with an optimization algorithm.The waveform shape is the experimental input and the droplet characteristics are the experimental output.This method has been used with a genetic algorithm [28] and a swarm-intelligence based technique [10].Alternatively, a highly efficient instanteneous adjoint-based approach to control the free surface inside the nozzle has been developed by [14], applied to numerical simulations.Our approach is similar, in that it considers a systematic approach to waveform optimization by using adjoint-based optimization.The differences are that this study includes the acoustics in the channel in addition to the flow in the nozzle, but makes an assumption on the spherical shape of the fluid/gas interface.
Adjoint-based optimization can be applied only to numerical simulations.The adjoint methods provide, in a single calculation, the gradient of an objective function with respect to all of the control parameters [29].This gradient is then used within a gradient-based optimization algorithm in order to converge to a local optimum.This gradient information greatly speeds up optimization.Adjoint-based optimization is faster than non-gradient-based methods when the number of control parameters exceeds the number of objective functions, which is the case here.It has been used in aerodynamics optimization [30], triggering in thermoacoustics [31], hydrodynamic stability [29] and shape optimization for hydrodynamic stability [32].
In this paper, we start from a moment just after pinch-off of the expelled ligament.We use adjoint-based optimization to optimize a velocity waveform at the actuator in order to reduce residual oscillations to zero within a given time.In order to be effective, this time needs to increase as the length of the printhead microchannel increases so that waves have time to reflect off the opposite end and return [7].In practice, this time would be the desired period between droplets.Similarly to previous studies we derive different governing equations for the channel and the nozzle.We develop a general approach to couple the channel to the nozzle through the boundary conditions on the surface between the two.We then derive the adjoint of this coupled channel-nozzle system and optimize the velocity waveform at the actuator for various channel shapes and optimization times.In doing so we reveal the physical mechanisms that are exploited in order to reduce the residual oscillations and define the minimum time between droplets.In this paper, we do not consider jetting, droplet coalescence, or entrainment of air bubbles into the nozzle.We also do not consider the actuation cost required to reduce these residual oscillations because this cost is much smaller than the actuation cost of jetting.This is because the jetting pulse creates a large amount of surface area and gives the jet a large kinetic energy, while the residual control movement has to absorb just a small amount of surface energy.

Formulation of the Direct Problem
We label the spatial domain Ω ⊂ R d with boundary ∂Ω = Γ i ⊂ R d−1 , where Γ i are the labels of sections of the boundary.We label the temporal domain T = {t : 0 < t < t f }, where t f is the final time.We label the mixed spatio-temporal domain of the problem Σ = (Ω × T ).We define the following scalar products in the spatial, temporal, and mixed domains: Here * denotes complex conjugation.An integral of a single variable over a domain is equivalent to the corresponding scalar product with a ≡ 1, ⟨b⟩ Ω ≡ ⟨1, b⟩ Ω .
We split the printhead microchannel into the non-overlapping channel Ω c and the nozzle Ω n domains (figure 1).The channel domain consists of the top part of the nozzle and a horizontal channel connected to vertical outlets.The nozzle domain is the bottom part of the nozzle, which includes the liquid/gas interface.The two domains intersect at a (virtual) flat boundary Ω c Ω n .
We couple the flows inside the channel and nozzle domains through velocity and stress boundary conditions on the shared boundary.On the nozzle domain shared boundary Γ N−C ⊂ Ω n , we prescribe the velocity, which is determined by the flow in the channel domain [13,33].On the channel flow shared boundary Γ C−N ⊂ Ω c , we prescribe the stress, which is determined by the force applied from the flow in the nozzle.
2.1.Governing equations for the acoustic oscillations in the channel 120 Using a two-parameter low Mach number asymptotic expansion [6,34] we separate the compressible Navier-Stokes equations into a stationary ambient state, an incompressible steady flow with no acoustic oscillation, and an acoustic oscillation with no steady flow.Then we consider only the acoustic oscillation with no steady flow.We expand the dimensional flow pressure, velocity and temperature (P * , u * , T * ) in terms of the acoustic Mach number ϵ ≡ ∥u * ∥ /c s ≪ 1: where P b P (0) is the ambient pressure, T b T (0) is the ambient temperature, and (u, P, T ) are the nondimensional acoustic velocity, pressure and temperature.We use Q c ≡ (u, P, T ) to denote the acoustic state.The speed of sound is set to c s = 1000ms −1 .
The acoustic state Q c inside the inkjet channel, Σ c ≡ (Ω c × T ), is governed by the linear nondimensional conservation equations for mass, momentum, and energy, with viscous and thermal losses (3).We include the acoustic density ρ = γ th P − T and entropy s = T / (γ th − 1) − P for convenience: The Reynolds and Peclet numbers are based on the speed of sound: Re = ρ b Lc s /µ, P e = ρ b Lc s c p /K, where ρ b is the fluid density, µ is the dynamic viscosity, K is the thermal conductivity, 125 and c p is the specific heat at constant pressure.The ratio of specific heats is γ th = 1.017, which is characteristic of water at 25 • [35].The viscous stress tensor is given by We rewrite the state equations in matrix form, where either the velocity, U i , or force, f i , is prescribed at the boundary.We apply homogeneous no slip and stress-free boundary conditions by setting U i = 0 on no slip boundaries Γ w , and f i = 0 on stress-free boundaries Γ open : where σ ij ≡ −P δ ij + Re −1 τ ij is the stress tensor.We apply isothermal boundary condition T = 0 on Γ w and Γ vel , and adiabatic boundary conditions ∂T ∂n = 0 on Γ open . 130 The total acoustic energy E ac [36], volume dissipation, R, and boundary energy flux, F of the thermoviscous acoustic problem (4) are: The volume averaged energy balance is: The boundary condition at the interface between the channel and the nozzle domains, Γ C−N is derived in section 2.2.The nozzle flow state is defined by the free surface curvature κ(t) which changes due to the mass flow through the shared boundary as where |Ω n (κ)| is the volume of the nozzle domain, h CL is the distance between the boundary Γ C−N and the static contact line in the nozzle, and γ is the surface tension coefficient.The channel flow boundary condition on Γ C−N is defined by the free surface curvature κ:

Reduced order model for the flow in the nozzle
In this paper, we start the optimization process from a moment just after an ejected ligament has pinched off and the fluid-gas interface has returned to the vicinity of the print head.In practice, the time between the end of the jetting pulse and the moment described above would need to be measured experimentally beforehand, and may be different for different jetting pulses.We model the boundary condition on Γ N−C as a prescribed velocity, determined by the acoustic oscillations inside the channel.Using the same strategy as for (2), we expand the dimensional nozzle flow variables in terms of the oscillating flow Mach number, ρ * ≃ ρ b 1 + ϵρ (1)

Conservation laws in the nozzle subdomain
We assume that the contact line between the solid, liquid and gas phases is stationary: u = 0 on ∂Γ free = Γ free Γ w , meaning that the energy of only the liquid-air interface Γ free contributes to the free surface energy.A vector field w defines the domain deformation due to the movement of the free boundary Γ free .The relationship between w and u depends on the type of boundary.
On Γ w both the domain and the fluid velocities have zero normal component (9a).On Γ free the flow velocity and the domain velocity are equal in the normal direction (9b).On Γ N−C there is no deformation (although there is mass flow through this boundary) (9c).w • n = ϵu (1) • n = 0 on Γ w . (9a) The nondimensional nozzle radius and the distance to the contact line are By choosing ϕ ≡ 1, the functional is the domain volume, J = ⟨1⟩ Ωn ≡ |Ω n |.After applying the boundary condition ( 9), the time derivative of the nozzle volume equals the normal flow at the free surface.
We apply the Reynolds transport theorem (10) to the continuity equation: We integrate the divergence term and apply the boundary conditions (9).The nondimensional mass balance equation is The force applied from the free surface to the nozzle flow is proportional to the curvature of the free surface, σ ij n j = −γκn i , where γ ≡ 2γ * / ϵρ b c 2 s R n is the nondimensional surface tension coefficient, and κ ≡ κ * R n /2 is the nondimensional curvature of the free surface.The 155 pressure drop along the nozzle due to the inertia of the fluid can be estimated as ρ b H CL ω ∥u∥, where ω ≃ 10 6 s −1 is the characteristic oscillation frequency of the flow.For ∥u∥ = 0.1 ms −1 , the nozzle pressure drop is 10 3 Pa which is equivalent to the pressure of a free surface with curvature radius R = 10 −4 m = 10R n .This shows that the pressure drop along the nozzle is small, but may not be negligible.
We need to account for the kinetic and potential energy of the fluid when the free surface curvature is small.The total energy of the nozzle system E n is a sum of the free surface energy E free , and kinetic K n and potential P n energy of the nozzle flow: Ωn .We assume that the temperature is uniform in the nozzle and neglect the thermal energy, so P n = 1 2 P (1) • P (1)  Ωn (and ρ (1) = P (1) ).The nozzle energy changes due to the viscous dissipation The time derivative of the free surface energy equals the energy flux through the free surface where κ ≡ κu (1) • n Γ free / u (1) • n Γ free is the effective curvature of the free surface.
Volume integrals in (12,13) scale as ⟨•⟩ Ωn ∼ h n r 2 n , and surface integrals scale as {•} ∼ r 2 n .In this study we consider cases in which the nozzle height is much shorter than the reference length, and expand the mass and energy conservation equations in terms of the small parameter h n ≪ 1.If h n L ≃ L, however, then the terms proportional to h n cannot be ignored and the flow should be modelled directly in the nozzle as well as in the channel.We approximate the flow velocity and the pressure by their values on the boundaries: which gives ρ (1) Ωn = P (1)  Ωn ≃ |Ω n | P (1)   Γ free Substituting that into (12) gives a mass balance approximation to second order in h n : The second term on the right hand side is proportional to h n and accounts for the nozzle flow compressibility and potential energy.
Substituting ( 14) into (13) and using (11,15) allows us to derive an energy-consistent approximation of the stress on the shared boundary Γ N−C :

Free surface parametrization
The above analysis is valid for fluid-gas interface of any shape.In this paper, inspired by experimental observations [39,40,41], we approximate the free surface as a spherical cap, Γ free ≃ Γfree , neglecting the presence of large wavenumber capillary waves k ≫ π/R n on Γ free .This simplifies the analysis by providing a single parameter (curvature) to quantify the relationship betwen pressure and nozzle mass, but would need to be revised for non-spherical interfaces, which are observed at the end of the jetting process ([3], section 3).The approximate surface curvature is a uniform function and coincides with the effective surface curvature: γ| Γ free ≃ γ| Γfree = κ(t).
The approximated free surface is a hemisphere when κ = 1.The surface area Γfree (κ) equals where cos θ(κ) ≡ √ 1 − κ2 .The nondimensional energy of the free boundary with uniform curvature equals Êfree = ϵ −1 γ rn 2 Γfree .The nozzle volume |Ω n | is a sum of the volume between the shared boundary Γ N−C and the plane of the free surface contact line πr 2 n h CL , and the volume enclosed between the plane of the free surface contact line and the free surface Ωn : The derivatives of the surface area and enclosed volume are functions of κ: The contact line is static d dt h CL = 0, and therefore The nozzle volume conservation equation (15) then becomes an ODE for the uniform curvature κ: subject to an initial condition κ(t = 0) = κ0 .Similarly to the stress boundary condition (17), the compressibility effects are accounted for only inside the static part of the nozzle domain.
The change of the nozzle system energy E n is consistent with the acoustic energy flux through Γ C−N : By combining the channel and the nozzle flows into one system, the total energy of the system changes due to the energy fluxes through the channel boundaries with prescribed velocity Γ vel or force Γ force , and the viscous and thermal dissipation R:

Direct problem discretization
We discretize the thermoviscous acoustic problem (4) in space using the finite elements package FEniCS [42].We choose continuous Lagrange elements of order 2 for the velocity and temperature components, and of order 1 for the pressure component of Q c [43,44], to construct the corresponding weak forms.
We divide the time domain T into N intervals of length ∆t.Following [45], the mid-point rule θ = 1 2 is chosen to discretize the unsteady thermoviscous acoustic equation by a finite difference scheme.This choice of θ is equivalent to the non-dissipative, dispersive, second order accurate Crank-Nicolson scheme [46].We choose a non-dissipative scheme because we want to minimize the error in the acoustic energy due to the numerical effects.
In order to reduce the cost of the unsteady computations, we choose such discretisation of the bilinear weak form that it is independent of time, and the matrix factorization could be stored and re-used each time step.This greatly reduces the computational time required to run a direct solver, and, as shown later, the adjoint solver backwards in time.Moreover, if one performs a direct-adjoint looping multiple times for optimization, the matrix still has to be factorized only once.The resulting linear system is solved with a direct solver MUMPS [47], for each time step.

Governing equations for the adjoint problems
Here we derive the adjoint counterpart of the coupled acoustic and nozzle flow system.The direct acoustic state Q c is defined in Σ c and governed by the thermoviscous acoustic equations (4) and boundary conditions (5).The nozzle flow state Q n is characterized by the free surface curvature ODE (21).The acoustic and nozzle flows are coupled via the boundary conditions (17) and the flow through the shared boundary.The acoustic energy and the flow energy inside the nozzle are defined by ( 6) and (22), respectively.We choose the total energy E = E ac + E n (23) at the final time t f as the objective function J : We multiply the acoustic state equations ( 4) by the adjoint acoustic variables Λ † c ≡ u † , P † , T † .We multiply the nozzle state equations ( 21) by ϵ −1 γκ † , where κ † is the adjoint curvature variable.
The augmented objective function is: We set the variation of the Lagrangian δL with respect to the direct state to zero.After successive integration by parts and applying the initial and boundary conditions from (5), we obtain the adjoint acoustic equations, and the adjoint nozzle flow equations, The adjoint acoustic operator B † c is defined as The adjoint acoustic stress tensor is defined as The adjoint no slip boundary conditions, are equal to the homogeneous boundary conditions of the direct problem [6].
In section 2.3, the direct channel and nozzle states were discretized in time using the Crank-Nicolson scheme.The adjoint state is approximated using the same method, resulting in a consistent time-discrete dual problem [48].The symmetry between the direct and adjoint problems is discussed in Appendix A.

The augmented gradient
The top boundary of the channel contains a piezo-electric actuator, which we model as a prescribed velocity boundary condition U(t).We minimize the objective function J (24) by optimizing the velocity profile on the actuator boundary Γ act .This control is described through a velocity Dirichlet boundary condition: The only non zero term in the Lagrangian variation (25) is the adjoint stress boundary integral on Γ act .This equals the objective variation with respect to the control, and therefore the objective gradient is: In other words, the distribution of the adjoint stress along the control boundary is the sensitivity distribution.

One dimensional test case
In order to illustrate the method and to discuss the physical mechanisms it exploits, we first apply the optimization technique in section 3 to a one-dimensional test case.A viscous acoustic flow inside a unit length domain Ω c = {x : 0 ≤ x ≤ 1} is initially at rest Q c (x, t = 0) = 0.The boundary at x = 1 is set as a control boundary.We prescribe a velocity profile u(x = 1) = U(t) on this boundary.The boundary at x = 0 is a free surface (21).This free surface is initially deformed κ(t = 0) = 0.05, and therefore possesses non-zero initial energy We define the nondimensional acoustic timescale t ac to be the time taken for a wave to travel from x = 0 to x = 1.In these units the time taken for a wave to travel from one side to the other and back is t L = 2.This is a key quantity that will be referred to later.We discretize the time domain T into equal intervals with time step ∆t = 10 −4 .The non-dimensional parameters of the experiment are provided in table 1.The optimization search space consists of control velocity values at each discrete time point U n for n = 1 . . .N − 1.The values of U 0 , U N are fixed to zero.We use the scipy.minimize(method='TNC')[49] implementation of the truncated Newton method as the gradient-based algorithm to minimize the objective function (24).In all experiments, we use the Taylor remainder convergence test to ensure that the sensitivities calculated with the adjoint method are correct.
We start with the case in which the final time is set to t f = 2.5 and use the gradient-based method to find the optimal waveform (figure 3a). Figure 3b shows the time history of the free surface energy E n (green line) and the acoustic energy E ac (blue line) in the optimally controlled case, normalized by the initial total energy value E(t = 0).The final energy at E(t = t f ) is 10 5 time lower than the initial value.In this simple case, the physical mechanism that it exploits can be clearly identified.The optimal waveform consist of three stages.The first stage is a pulse lasting τ p = t f − t L = 0.5 that withdraws half the volume stored in the nozzle domain.In the second stage the actuator remains inactive while the front and back of the pulse reach the free surface at t = 1 and t = 1 + τ p respectively.The pulse velocity amplitude doubles as it reflects from the free surface.The amount of fluid transferred through the free surface is therefore equal to the volume initially stored in the nozzle domain.At t = 2 the front of the reflected pulse reaches the actuator.Between these two times the free surface relaxes to zero curvature and, in doing so, reflects an acoustic pulse back towards the actuator.The third stage is a pulse lasting τ p = t f − t L = 0.5 that withdraws more fluid such that the reflected pulse leaves the channel without further reflection and returns the fluid in the channel to its initial state at exactly t = t f .
The two pulses, when combined, withdraw exactly the mass of fluid between the initial and final positions of the surface.
For comparison we then examine the case in which the final time is set to t f = 3.0 (figure 4). Figure 4b shows the time history of the free surface energy E n (green line) and the acoustic energy E ac (blue line) in the optimally controlled case, normalized by the initial total energy value E(t = 0).The optimal waveform is qualitatively identical to that found when t f = 2.5 but with pulses lasting τ p = t f − t L = 1.0 rather than τ p = 0.5.As before, the actuator withdraws exactly the half mass of fluid between the initial and final positions of the surface, but this time with a longer pulse and, consequently, with a smaller actuator velocity.Consideration of the physical mechanism shows that the optimization time has to be greater than t f = 2 and that, beyond that, its lower limit will be determined by the maximum speed of the actuator.In both cases, the spurious oscillations in U appear due to step-like control at initial time (Gibbs phenomena).

Two dimensional straight channel
Having shown that the optimization algorithm works for a simple 1D case, and having identified the physical mechanism in that case, we now examine a 2D straight channel with a nozzle placed at the centre of one wall, an actuator along the opposite wall, an outlet boundary at the left side, and a symmetry boundary at the right side (figure 5).The dimensional parameters of the nozzle domain and acoustic constants are given in table 2. The spatial domain is discretized into 120 • 10 3 triangular elements [50].In the remaining sections, we present the time and space scales in dimensional format because this enables easier comparison with the literature.
We assume that the nozzle state Q n shortly after the droplet has been expelled is a free surface with uniform curvature.For ease of demostration, we assume that the acoustic energy in the channel is zero.This assumption could easily be relaxed in practice.The simulation therefore starts from the zero acoustic state and non zero curvature: The objective is to minimize the total energy at the final time Typically the time between droplets is between 2 µs and 20 µs.In this paper, we calculate the open loop forcing that would need to be applied for periods of 1 µs, 2 µs and 3 µs in order to minimize the total energy at the end of the forcing period.We discretize the time domain    T into equal time intervals with time step ∆t = 10 −3 µs.This time step is chosen in order to have sufficient time resolution of the acoustic motion inside the narrowest part of the channel Ω c , near the nozzle boundary.
The actuating waveforms in inkjet printing are constrained by the limitations of the driving electronics and response of the piezo-electric actuator.In this section, we assume that the piezoelectric actuator moves as a solid plate in the direction normal to the channel's wall.We model this as a boundary velocity U that is spatially uniform along the control boundary Γ act .There is no mass flow through the actual physical actuator boundary, but we prescribe a velocity boundary condition on it.So there is effective mass flow through the fixed boundary used in computations.
In 4.4, we compare the spatially uniform boundary control with a parabolic actuator velocity profile.
In practical devices the electric signal that forces the piezo-electric actuator is piecewise linear with a temporal resolution, w, between 0.01 µs and 0.1 µs.We use a continuous piecewise linear function in time to describe the boundary velocity.
In our model, shown in figure 5, the length of the control element (actuator boundary) varies from L act = 20µm to L act = 200µm.The waveform time resolution is fixed at w = 0.1µs.
The left boundary is stress-free.The right boundary of the computational domain (shaded) is a symmetry plane.Figure 6 shows the optimized total energy at the final time E n (t f ) normalized by the uncontrolled total energy at the final time E * (t f ).For 1µs optimization time there is Normalized energy, almost no energy reduction because, as will be shown later, there is insufficient time to control the wave reflected by the left boundary.For 2µs optimization time the final energy reduces by one order of magnitude.For 3µs optimization time the final energy reduces by a further order of magnitude.The actuator size has little influence as long as it exceeds 100µm.The physical reasons for this are explored next.
For illustration, we will examine the results with final time t f = 2µs and actuator length L act = 200µm.Figure 7a shows the nozzle surface energy (green) and the total energy (blue) for the uncontrolled (dashed) and controlled (solid) cases.Figure 7b shows, for the uncontrolled case, the integrated acoustic energy dissipation (red dashed) and, for the controlled case, the total energy E(t) (blue), the integrated energy flux through the actuator boundary F act (τ )dτ (purple), the integrated acoustic energy dissipation R(τ )dτ (red solid).Due to the energy balance (7), these are related by Figure 7c shows the mass flux through the actuator boundary (black) and the nozzle boundary (red).Figure 8 shows snapshots of the pressure field of the controlled case at times corresponding to the empty circles in figure 7a.
For the uncontrolled case, the nozzle surface energy E * n (green dashed line in figure 7a) reduces smoothly as the free surface relaxes.As for the one-dimensional test case, this sends an acoustic wave down the channel, increasing the acoustic energy.The total energy E * (blue dashed line, figure 7a), which comprises the nozzle surface energy and the acoustic energy, reduces gently as the wave dissipates due to thermo-viscous mechanisms.
At time t = t f = 2.0 µs, the controlled case has almost 20 times lower energy than the uncontrolled case (compare the solid blue and dashed blue lines in figure 7a).The optimal waveform (black line in figure 7c) consists of three phases.During the first phase A + : 0 ≤ t ≤ 0.38 µs the actuator pulls fluid upwards and creates a negative pressure wave (figure 8a).This wave moves down towards the nozzle and left along the channel.The wave reaches the nozzle boundary at t = 0.12µs, at which point mass starts to be pulled out of the nozzle (red solid line in 7b), and the nozzle surface energy starts to reduce (green solid line in 7a).The wave reflects back off the nozzle, reaching the actuator at t = 0.24µs.Reverberations at this timescale, which is that of the height of the channel, continue during the controlled period.This behaviour is   similar to that of the 1D test case but is more complicated because the flow is 2D.This wave relaxes the free surface but, unavoidably, produces a large amplitude wave moving left along the channel (figure 8b).Indeed the energy flux through the control boundary during the first phase, A + , is large and the total energy rapidly increases to nearly three times that of the uncontrolled case.
During the second phase 0.38 ≤ t ≤ 1.62 µs the integrated mass flux through the actuator, M act ≡ {u • n} Γact , is almost identical to the integrated mass flux away from the nozzle, M n (black and red solid lines in figure 7c).Compared with the first pulse, this motion is relatively slow, shown by the fact that the energy flux through the actuator boundary is small ∂ t E ≃ 0.
The nozzle surface energy reduces to nearly zero during this phase (green solid line in 7a).
Meanwhile, in the channel, the pressure wave generated during the first phase reflects off the stress-free boundary and a positive pressure wave travels back towards the nozzle and the actuator (figure 8c).The total energy E steadily decreases due to viscous and thermal dissipation of the acoustic wave (figure 7b, red line).
The third phase A − : 1.62 ≤ t ≤ 2.0 µs is the counterpart of the first phase A + .When the positive pressure wave reaches the symmetry plane of the channel the actuator quickly moves out of the domain again (black line in figure 7b) and optimally absorbs the acoustic energy (blue  line in figure 7b) by moving to make the wave do work on the actuator boundary (purple line in figure 7b).The acoustic pressure quickly reduces and remains small thereafter (fig.8f).This acts on the timescale of the channel: 2L/c s = 2×900µm/ (1000m/s) = 1.8µs.If the optimization time is 1µs, the nozzle free surface energy could be reduced by the actuator but there would then be insufficient time to absorb the acoustic wave that is reflected off the stress-free boundary.The optimal solution is to do almost nothing.
We now investigate the effect of the actuator size L act on the final energy and the optimal waveform.Figure 6 (orange line) shows the final total energy of the controlled cases, normalized by the final energy of the uncontrolled case.Figure 9 shows the optimal mass flux through the actuator boundary for different actuator sizes, L act .All waveforms have the three-phase shape described above and exploit the same mechanism.The mass flux increases with the actuator length.This is because the actuators produce a wave that both pulls mass out of the nozzle and propagates to the stress free boundary, where it pulls mass through that boundary.The nozzle mass flux is determined by the mass in the nozzle and, with that mass flux fixed, the boundary mass flux then increases with the actuator length.Shorter actuators are therefore more efficient.

Waveform optimization for t
For illustration, we will examine the results with final time t f = 3µs and actuator length Normalized energy Normalized energy  The first phase A + : 0 ≤ t ≤ 0.55 µs is similar to the A + phase of the t f = 2µs case.The actuator pulls fluid upwards and generates a negative pressure wave (figure 11a).As before, the wave reaches the nozzle boundary at t = 0.12µs and mass starts to be pulled out of the nozzle.This wave reflects back off the nozzle.The reverberations in the fluid are as strong as in the previous case, but the actuator is wider.The wider actuator is unable to reduce the reverberations' amplitude without simultaneously creating high amplitude waves elsewhere.

Timescale of length of channel
Therefore the actuator moves less in order to avoid creating these high amplitude waves elsewhere.As before, this motion produces a large amplitude wave moving left along the channel.
During the second phase B + : 0.55 ≤ t ≤ 1.2µs the mass flux through the nozzle reaches a maximum.The actuator is still moving upwards, partially compensating for the flow from the nozzle and slowly absorbing the acoustic energy.A positive pressure wave moves left along the channel (fig.11b).(The optimization algorithm does not create this wave for the 2µs case because there is insufficient time to cancel it.)During this phase, almost all of the fluid is quickly transferred from the nozzle to the channel.By the end of the B + phase, the free surface has nearly reached its final low energy state.
During the third phase 1.2 ≤ t ≤ 1.75µs the mass flux from the actuator broadly cancels that from the nozzle.No new pressure waves form during this phase (fig.11c).The negative pressure left-running wave reflects from the open end (11d) and becomes a positive pressure right-running wave.This middle phase is similar to the middle phase of the t f = 2µs case.
The fourth phase A − : 1.8 ≤ t ≤ 2.3µs is the same as the A − phase in the t f = 2µs case.The positive right-running wave is absorbed by the actuator (fig.11e, 11f).This results in a rapid decrease in total energy.Meanwhile the positive left-running wave has reflected from the open end and has become a negative right-running wave.During the fifth phase B − : 2.3 ≤ t ≤ 3.0µs, the negative left-running pressure wave reaches the actuator, and is optimally absorbed by doing work on the actuator boundary.Figure 12 shows the mass flow at the actuator M act as a function of time for different L act .The shapes of the waveforms for t f = 3µs with different actuator lengths are similar.The increase of mass flux with actuator length is less pronounced because the combination of A + , B + waves pulls less fluid through the stress-free boundary than the A + wave alone.
In summary, the most effective control to minimize the energy at the final time is achieved when the actuator moves out of the domain and provides a negative pressure pulse to accelerate the flow at the nozzle, and then adapts to the large mass flux during the B + phase.This combination efficiently transfers all fluid from the nozzle to the channel and leaves sufficient time for the actuator to absorb the reflected waves afterwards.At lower optimization times, reasonably effective control can be achieved with a short negative pulse, followed by a long middle period in which the channel slowly absorbs the fluid from the nozzle, leaving sufficient time to cancel its reverberation from the ends of the channel.If there is insufficient time to cancel the reverberation from the ends of the channel then any control is ineffective.This shows that the minimum optimization time (i.e. the minimum time between droplet ejections) is 2L/c s .

Two dimensional U-shaped print head channel
Having shown that the optimization algorithm works for a straight channel, and having highlighted the importance of wave reflections at the ends of the channel, we now examine a realistic case, in which the ends of the channel bend upwards to make a U-shape and the length along the centreline is L = 1235µm (figure 13).We set the actuator length to L act = 400µm.
From the snapshots (figure 15) we see that the waves disperse slightly as they travel round the corner.By comparing figure 14 with figure 10, however, we see that the optimal profiles for the U-shaped channel are qualitatively identical to those for the long straight channel in section 4.2.2.From this, and the snapshots, we deduce that the optimization method is exploiting the same physical mechanism.Normalized energy Finally we investigate how the objective value changes when we increase the waveform time resolution from w = 1.0 to w = 0.5, and 0.25µs.These correspond to the time resolution of stateof-the-art piezoelectric controllers.We project the optimal waveform with w = 1.0µs (figure 16b, solid line) to a higher dimensional space with w = 0.5µs, and use the projected solution as the initial guess for a new optimization problem.The optimal waveform for w = 0.5µs (figure 16b, dash-dotted line) results in 25% lower objective value (figure 16a).The optimal waveform for w = 0.25µs (figure 16b, dotted line) results in further 4% reduction in the objective value.These improvements are quite small, showing that the rather basic sinusoid-type waveform for the w = 1.0µs case provides a good trade-off between efficiency and complexity.This waveform is an outward moving pulse lasting 2.7µs, and an inwards moving pulse lasting 2 µs.The outward pulse causes the free surface to relax but generates acoustic waves that travel down the channel.These reflect and first arrive back at the actuator at t = 2.4µs.The trailing edge of the outward pulse and then the inward pulse absorb these reflected waves optimally.Given that this waveform would be imposed just after ejection of the droplet, it is reassuringly similar to the W-shaped waveform that, by trial and error, has been shown to remove residual acoustic waves arising from the previous ejection cycle [24].The total energy reduces to nearly zero, so if there is another local minimum, it cannot be significantly better than this one.The advance in this paper is to show this rigorously with adjoint-based optimization in the time domain, to identify the physical mechanisms that this waveform exploits, and to employ a general method that can be applied to other geometries.

Optimization with a parabolic actuator velocity profile
In previous sections, the flow inside a straight and a U-shaped channels driven by an actuator with a uniform actuator velocity profile has been optimized.This section examines a system with a parabolic actuator velocity profile, and two cases are considered: a straight channel with L act = 200µm, t f = 2µs (similar to 4.2), and a U-shaped channel with L act = 400µm, t f = 5µs (similar to 4.

3). The velocity boundary condition on Γ
x is the coordinate along the actuator boundary, and x act is the position of the leftmost point on the actuator boundary.
The same optimization strategy as descried in section 4.1 is applied to find an optimal waveform.Figure 17 shows an optimal waveform, and a mass and energy fluxes for a straight channel (figure 5) and t f = 2µs (dashed lines).Solid lines denote the optimized results for a flat actuator velocity profile (from section 4.2).The optimal waveforms (figure 17c) for parabolic and flat velocity profiles are very similar, and the energy transfer between the nozzle, channel, and the actuator follow the same pattern (figure 17b).
Figure 18 shows optimization results for a U-shaped printhead channel (the same as in section 4.3) with a parabolic actuator velocity profile.Again, the parabolic profile results (dashed lines) match very closely the optimization results for a uniform profile (solid lines).This means that the energy damping mechanism is independent of the exact shape of the actuator velocity profile, and is governed by the mass flow through the actuator boundary M act .

Conclusions
In this paper we develop a gradient-based approach that uses adjoint methods in the time domain to optimize the deformation of the piezo-electric actuator in an inkjet microchannel.
We show that an optimally controlled actuator can reduce the total energy inside the printhead microchannel geometry by 1000 times, compared with the uncontrolled case.The actuator's initial movement withdraws fluid from the nozzle.The liquid/gas surface relaxes towards the zero curvature state.This initial movement sends an acoustic wave down the microchannel, which reflects off the open end and returns some time later.The actuator then moves in order to perfectly absorb the reflected acoustic wave.The minimum time over which optimization can be successful, which is therefore the minimum time between droplets, is the time taken for an acoustic wave to travel from the actuator to the open end and back.The duration of the waveform itself must be added to this time.Short waveforms (e.g.A + in figure 7c) are in one direction only and reduce the energy by just over one order of magnitude.Longer optimization times allow sufficient time for the waveform to have two components in opposite directions (e.g. A + and B + in figure 10c).This waveform reduces the free surface energy more quickly and by over a further order of magnitude.This paper demonstrates the success of a more systematic approach than that currently used in industry: adjoint-based optimization of the waveform using numerical simulations, which can be interpreted physically.At a minimum, this physical understanding is useful in narrowing down the range of waveforms tested experimentally by trial and error.At best, this technique shows how to use numerical simulations in order to systematically and efficiently find the optimal waveform and the minimum time between drops in inkjet print heads.
This method can easily be applied to microchannels with complex geometry in three dimensions.The same method can be used to perform waveform optimization for the droplet ejection, in order to control the droplet volume and momentum.This would require a new nozzle model and a new objective function, but the channel model and optimization method would remain the same.Validation of the formulation and the numerical solution is the subject of future work.

Figure 1 :
Figure 1: Channel (yellow) and nozzle (blue) domains of a single microchannel within an inkjet print head.In this paper the channel is modelled as 2D planar and the nozzle is modelled as 2D axisymmetric.The piezo-electric actuator sits on the channel wall opposite the nozzle.

Figure 2 :
Figure 2: Fluid inside the droplet domain is bounded by the solid walls Γw and the free surface Γ free , and connects to the channel through Γ N−C .

Figure 2
Figure 2 shows the nozzle domain Ω n , which is an axisymmetric tube of radius R n .The distance between the shared boundary and the multiphase interface at rest (the contact line) is H CL = 10µm.Both R n and H CL are small in comparison to the characteristic channel length, L = 100µm.The droplet domain is bounded by three types of boundary: the boundary Γ N−C between the nozzle and the channel domains, the no slip wall boundary Γ w , and the moving free surface Γ free .The boundary Γ free is a multiphase interface between the liquid inside the nozzle, the nozzle walls and the outside air.For a typical nozzle radius R n = 10µm, fluid density ρ b = 10 3 kg/m 3 , and surface tension coefficient γ * = 50 × 10 −3 Nm −2 , the capillary time scale is t γ ∼ ρR 3n /γ * = 10 −5 s[37,38].The acoustic time scale, t ac ∼ L/c s = 10 −7 s, where c s ≈ 1000ms −1 is the speed of sound.The acoustics therefore occur two orders of magnitude more quickly than the capillary motion, indicating that the capillary motion is quasi-static compared with the acoustic motion.In this paper we model both the acoustic and nozzle flows on the acoustic timescale.
and h CL ≡ H CL /L.The nondimensional nozzle volume |Ω n (t)| is the sum of a static component, πr 2 n h CL , and a time-varying component Ωn (t), where |•| means surface or volume measure according to the context.A nondimensional effective height of the channel is defined as h n ≡ |Ω n | /(πr 2 n ).The algebraic relation between the domain volume and the domain deformation, although obvious, can be derived rigorously with the Reynolds transport theorem by considering the time derivative of the volume functional J = ⟨ϕ⟩ Ω :

1 ) i
Ωn and the energy flux through the shared boundary Γ N−C are: ) Here ∆ Γ is the tangent Laplace operator on Γ N−C .The term on the right hand side proportional to h n accounts for the internal and viscous effects in the fluid.The nozzle flow velocity on Γ N−C is determined by the acoustic velocity on Γ C−N .Considering the inertial and viscous effects only in the static part of the nozzle domain h CL , we obtain a nonlinear (through the curvature term κ) acoustic impedance boundary condition:

Figure 3 :
Figure 3: Optimally controlled case of the one-dimensional unit length domain, t f = 2.5.(a) the optimal velocity U of the control boundary, where τp indicates the pulse duration, and t L is the timescale of the length of channel;(b) the nozzle energy (green), and acoustic energy (blue), normalized by the initial total energy value E(t = 0).

Figure 4 :
Figure 4: Optimally controlled case of the one-dimensional unit-length domain, t f = 3.(a) the optimal velocity U of the control boundary, where τp indicates the pulse duration, and t L is the timescale of the length of channel; (b) the nozzle energy (green), and acoustic energy (blue), normalized by the initial total energy value E(t = 0).

Figure 6 :
Figure 6: Optimized objective values for eight actuator lengths, Lact, and three final times, t f .The values are normalized by the total energy at final time in the uncontrolled case.

Figure 7 :
Figure 7: Comparison between the uncontrolled case (dashed lines) and the optimally controlled case (solid lines) with Lact = 200µm and t f = 2µs.(a) Total energy E(t) (blue) and nozzle energy En(t) (green), normalized by E(t = 0).(b) Integrated energy flux through the actuator boundary Fact(τ )dτ (purple line), integrated acoustic energy dissipation R(τ )dτ (red solid line), and the integrated acoustic energy dissipation R * (τ )dτ in the uncontrolled case (red dashed line).(c) Boundary mass flux through the actuator boundary Mact (black) and the nozzle boundary Mn (red).The red dashed line is the mass flux M * n through the nozzle boundary in the uncontrolled case.The coloured patches denote the actuation phases, when the acoustic waves are formed (A + ) and absorbed by the actuator (A − ).

Figure 8 :
Figure 8: Snapshots of the pressure distribution inside an injector channel at different times, with the optimal control applied to the actuator boundary, Lact = 200µm, t f = 2µs.The actuator is the thin horizontal rectangle at the top-right of each frame.

Figure 10 Figure 9 :
Figure 9: Optimal mass flux through the control boundary as a function of time, for different actuator lengths.The control duration is 2 µs.

Figure 11 :Figure 12 :
Figure 11: Snapshots of the pressure distribution inside an injector channel at different times, with the optimal control applied to the actuator boundary, Lact = 100µm, t f = 3µs.channel T 3 L 100.mp4 is a movie of the case (see supplementary materials).

Figure 14 :
Figure 14: Comparison between the printhead uncontrolled case (dashed lines) and the optimally controlled case (solid lines) with Lact = 400µm and t f = 5µs.(a) Total energy E(t) (blue) and nozzle energy En(t) (green), normalized by E(t = 0).(b) Integrated energy flux through the actuator boundary t 0 Fact(τ )dτ (purple line), integrated acoustic energy dissipation R(τ )dτ (red solid line), and the integrated acoustic energy dissipation R * (τ )dτ in the uncontrolled case (red dashed line).(c) Boundary mass flux through the actuator boundary Mact (black) and the nozzle boundary Mn (red).The red dashed line is the mass flux M * n through the nozzle boundary in the uncontrolled case.

Figure 15 :
Figure 15: Snapshots of the pressure distribution inside a U-shaped printhead at different times, with the optimal control applied to the actuator boundary, Lact = 400µm, t f = 5µs.printhead T 5 L 400.mp4 is a movie of the case (see supplementary materials).

Figure 16 :
Figure 16: (a) Optimized objective values (total energy at the final time E(T ) normalized by initial value E(0)) for different waveform resolution w, and (b) corresponding waveforms.

Figure 17 :
Figure 17: Comparison between the optimally controlled straight channel flow with flat (solid lines) and parabolic (dahsed lines) actuator velocity profiles with Lact = 200µm and t f = 2µs.(a) Total energy En(t) (blue) and nozzle energy En(t) (green), normalized by E(t = 0).(b) Integrated energy fluxes through the actuator boundary Fact(τ )dτ (purple line), integrated acoustic energy dissipation R(τ )dτ (red solid line).(c) Boundary mass fluxes through the actuator boundary Mact (black) and the nozzle boundary Mn (red).The coloured patches denote the actuation phases, when the acoustic waves are formed (A + ) and absorbed by the actuator (A − ).

Figure 18 :
Figure 18: Comparison between the optimally controlled printhead flow with flat (solid lines) and parabolic (dahsed lines) actuator velocity profiles with Lact = 400µm and t f = 5µs.(a) Total energy E(t) (blue) and nozzle energy En(t) (green), normalized by E(t = 0).(b) Integrated energy fluxes through the actuator boundary t 0 Fact(τ )dτ (purple line), integrated acoustic energy dissipation R(τ )dτ (red solid line).(c) Boundary mass fluxes through the actuator boundary Mact (black) and the nozzle boundary Mn (red).

Table 1 :
Parameters of the optimization test case 4.1