Objective function choice for control of a thermocapillary flow using an adjoint-based control strategy

International Journal of Heat and Fluid Flow 56 (2015) 28–42

Contents lists available at ScienceDirect

International Journal of Heat and Fluid Flow journal homepage: www.elsevier.com/locate/ijhff

Objective function choice for control of a thermocapillary flow using an adjoint-based control strategy Frank H. Muldoon ⇑, Hendrik C. Kuhlmann Institute of Fluid Mechanics and Heat Transfer, Vienna University of Technology, Resselgasse 3, A-1040 Vienna, Austria

a r t i c l e

i n f o

Article history: Received 7 April 2014 Received in revised form 17 February 2015 Accepted 2 June 2015

Keywords: Thermocapillary flow Marangoni effect Control Optimization Impurities Striations in crystal growth

a b s t r a c t The problem of suppressing flow oscillations in a thermocapillary flow is addressed using a gradient-based control strategy. The physical problem addressed is the ‘‘open boat’’ process of crystal growth, the flow in which is driven by thermocapillary and buoyancy effects. The problem is modeled by the two-dimensional unsteady incompressible Navier–Stokes and energy equations under the Boussinesq approximation. The goal of the control is to suppress flow oscillations which arise when the driving forces are such that the flow becomes unsteady. The control is a spatially and temporally varying temperature gradient boundary condition at the free surface. The control which minimizes the flow oscillations is found using a conjugate gradient method, where the gradient of the objective function with respect to the control variables is obtained from solving a set of adjoint equations. The issue of choosing an objective function that can be both optimized in a computationally efficient manner and optimization of which provides control that damps the flow oscillations is investigated. Almost complete suppression of the flow oscillations is obtained for certain choices of the objective function. Ó 2015 Elsevier Inc. All rights reserved.

1. Introduction One goal behind mathematical or experimental modeling of fluid and heat transfer problems is the design or control of engineering systems. Such design or control inherently carries with it the concept of optimization, as the goal is to select control inputs or design parameters such that the engineering system is optimal according to a prescribed metric. Optimization of engineering systems does not require a mathematical model of the system, however, the existence of such a model can greatly aid optimization in two ways. One way is a much shorter time with which a solution satisfying the mathematical model can often be obtained compared to that required to carry out experiments. The other is that a mathematical model potentially can contain information such as derivatives and sensitivities with respect to the quantities being varied that are very difficult, if not practically impossible, to obtain from experiments. The numerical methods used in the area of Computational Fluid Dynamics to obtain approximate solutions to mathematical models such as the Navier–Stokes and energy equations and computing resources have advanced to a stage where many fluid and heat-transfer problems of engineering interest may be reliably solved quickly enough that they can be used in all stages of the design cycle. This has enabled Computational Fluid ⇑ Corresponding author. E-mail address: [email protected] (F.H. Muldoon). http://dx.doi.org/10.1016/j.ijheatfluidflow.2015.06.001 0142-727X/Ó 2015 Elsevier Inc. All rights reserved.

Dynamics to replace expensive physical prototype models and has greatly reduced design cost and time. This corresponds to the first way that mathematical modeling can aid optimization, in which the mathematical model is used as a ‘‘drop in’’ replacement for experimental models. A more efficient means of using a mathematical model is its combination with an optimization method which uses higher order information about the physical system such as derivatives that can be obtained from the mathematical model. A field with a high demand of optimization is crystal growth. A major problem in crystal-growth processes with a free surface is the appearance of microscopic striations which are the results of impurity segregation. The striations arise due to unsteadiness in the growth process at the crystallization front which can arise from temperature fluctuations in the molten material (Müller and Ostrogorsky, 1994). Flow and temperature fluctuations, despite time-independent boundary conditions, are typically due to instabilities of the flow which is driven by thermocapillary and buoyancy forces which cannot be avoided. Thermocapillary surface forces arise along the interface (free surface) between the molten material and the surrounding atmosphere due to the thermocapillary effect which is caused by the dependence of the surface tension on temperature (Kuhlmann, 1999). The goal of the present work is to investigate issues involved with combining a mathematical model of the unsteady Navier– Stokes and energy equations with an optimization method with the aim of controlling a thermocapillary flow including the

F.H. Muldoon, H.C. Kuhlmann / International Journal of Heat and Fluid Flow 56 (2015) 28–42

29

Nomenclature

a

distance along line in search direction aspect ratio, horizontal over vertical distance / control variables /max ð/min Þ maximum (minimum) value of the control variables X the spatial computational domain @X the boundary of X @ Xfs the part of @ X containing the free surface, @ Xfs # @ X H objective function used by the minimizer Hp penalty part of objective function h vertical distance p pressure

C

temperature field that is a model of a manufacturing process for crystals. The control in question is a spatially and temporally varying normal temperature gradient at the free surface of the thermocapillary flow. The objective is to find control that suppresses an unsteady flow due hydrodynamic instabilities which arise for Reynolds numbers greater than a certain critical value. The time scales of slightly supercritical thermocapillary flows are much larger than those of high-Reynolds-number turbulent flows. This makes them an ideal target for real-time control, since the time available to carry out the computations needed to find optimal control over some period of time will be substantially smaller than the real (wall clock) time available to carry out the computations. Note, however, that real crystal growth is much more complicated than the model flow studied here (se e.g. Hurle, 1994). Therefore, our study of a very much simplified model flow can only be a first step towards a real industrial application. A literature review concerning the combination of optimization/control with Computational Fluid Dynamics in general exists in Muldoon (2013). To avoid repetition, we briefly mention here only the work on thermocapillary flows. Thermocapillary systems that have been the subject of investigation of flow control are thermocapillary liquid layers (as a model for the open-boat technique) and thermocapillary half-zones (as a model for the floating-zone technique). Benz et al. (1998) measured the free-surface temperature variation caused by flow oscillations in thermocapillary liquid layers. They have experimentally proven the feasibility of suppressing these oscillations by heating the free surface along lines parallel to the isolines of constant phase of the oscillations. In a series of papers Shiomi and co-workers have experimentally demonstrated the feasibility of feedback control of thermocapillary flows by suppression of oscillations in an annular pool resembling the Czochralski process (Shiomi et al., 2001; Shiomi and Amberg, 2002) and in the half zone-model model of the floating-zone technique (Shiomi et al., 2003). In all experiments temperature signals were recorded from two sampling points on the free surface and the control was applied by point heating the free surface at one or two other points. For the annular pool, Shiomi and Amberg (2005) confirmed their experimental results by simulating some representative cases. In the present work we consider the control of the two-dimensional thermocapillary flow in an open boat. In contrast to previous investigations we make use of the full information provided by the free surface temperature and control is applied at the whole free surface. In a work closely related to the present work Muldoon (2013) demonstrated control of flow oscillations in this geometry. The mathematical model, the two-dimensional unsteady incompressible Navier–Stokes equations and energy equations, was the same as in the present work. It was found that the length of the ‘‘event horizon’’ (i.e., the distance ahead in time over which optimal control was found), played a significant role

T t t EH tA u u v x y

temperature time event horizon length temporal advancement length, the extent of the temporal dimension of /; t A 6 t EH velocity vector first component of velocity vector (i.e., u1 ) second component of velocity vector (i.e., u2 ) spatial coordinate in horizontal direction spatial coordinate in vertical direction

in the success of finding control that suppressed the flow oscillations. The present work represents a continuation of Muldoon (2013) focusing on investigation of different objective functions and issues involved with solving the optimization problems required to determine the optimal control. The present work uses a conjugate-gradient algorithm to find optimal control. Such a method requires the gradient of the objective function to be minimized with respect to the control input. Therefore, a key issue is the means by which this gradient is obtained. In the present work this gradient is obtained by solving an additional set of adjoint equations derived from the governing mathematical model. Aside from gradient-based optimization methods, there exist methods such as simulated annealing and genetic algorithms that do not require a gradient (Davis, 1987). Such methods are very useful if the gradient of the objective function cannot be obtained or does not exist, but since they lack the information about the function being minimized contained in the gradient, they typically require many more evaluations of the objective function than gradient-based methods. This represents a serious disadvantage for the solution of the mathematical problem considered in the present work, which, as it involves solving the unsteady Navier–Stokes and energy equations over time periods of significant length, is extremely expensive to compute. In the present work, the computational expense in terms of arithmetical operations required for computing the gradient by means of solving the adjoint equations is similar to that required to compute the objective function. For this reason, a gradient-based optimization method is a more efficient approach for optimization for the present problem. 2. Formulation of the problem A schematic with dimensions and boundary conditions of the physical system to be investigated is given in Fig. 1. This system is known as an ‘‘open boat’’ in the literature. The problem contains a free surface, at which, due to gradients in the surface tension

Fig. 1. Schematic of physical domain for the minimization problem showing boundary conditions and dimensions.

30

F.H. Muldoon, H.C. Kuhlmann / International Journal of Heat and Fluid Flow 56 (2015) 28–42

arising from its dependence on temperature, there exists a stress tangential to the free surface, i.e., the thermocapillary effect. In addition, there exists a buoyancy force due to the effect of gravity on changes in the density arising from density gradients due to thermal expansion. The driving force for the flow is the temperature difference between the vertical walls, with the cold left and hot right walls maintained at T C and T H , respectively. The problem is simplified by assuming the surrounding gas is sufficiently dilute such that its effect on the liquid can be neglected and that the liquid is incompressible with constant fluid properties except for the surface tension and density. The surface tension

r ¼ r0  cðT  T 0 Þ;

ð1Þ

where T 0 ¼ ðT C þ T H Þ=2, is assumed to depend linearly on the temperature with surface tension coefficient c. Eq. (1) is a good approximation for typical temperature differences in experiments. The assumption is made that the mean surface tension r0 is large enough that flow-induced deformations of the free surface are vanishingly small and can be neglected. This implies that the capillary number cðT H  T C Þ=r0 is vanishingly small. These assumptions remove the mean surface tension from the mathematical model. A more detailed explanation can be found in Kuhlmann (1999). The density of the liquid is also assumed to depend linearly on the temperature and for its variations to be important only when multiplied by the gravitational acceleration, i.e., we model the flow using the Boussinesq approximation of the Navier–Stokes and energy equations. A gross feature of the flow is a surface-tension-gradient driven flow along the free surface from the hot towards the cold wall and a return flow near the lower wall. Upon increasing the driving force from zero, a time-invariant multicellular structure arises with a row of cells stretching from the hot to the cold walls. For Pr ¼ 14:8 (Shevtsova and Legros, 2003) the strength of the co-rotating vortices increases towards the hot wall (for low Prandtl numbers, see, Ben Hadid and Roux, 1990; Laure et al., 1990). Upon increasing the thermocapillary driving force beyond a critical value the flow becomes time-periodic (Peltier and Biringen, 1993; Xu and Zebib, 1998) with the cells and the associated temperature variations moving back and forth parallel to the free surface. Movie 1 shows an example for the uncontrolled oscillatory flow. These flow oscillations are the target of the control which we intend to suppress by imposing a suitable temperature gradient normal to the free surface. 2.1. Optimization problem In the present work, square brackets ½  denote functional dependence, while parentheses ðÞ denote multiplication; i.e., a½b þ c means a is a function of b þ c , while aðb þ cÞ means a times b þ c. The velocity vector is denoted by u, the individual components of which are defined by u ¼ ðu; v ÞT . We consider a rectangular domain X of height h and width-to-height ratio C (Fig. 1) occupied by the liquid and with a free surface @ Xfs . The control is to be applied along the entire free surface. The goal is to minimize the strength of oscillations in X arising from buoyancy forces and the thermocapillary stress by   varying the applied normal temperature gradient / @ Xfs ; t at the free surface. Mathematically this problem can be posed as Minimize an objective function

    H q / @ Xfs ; t on X  ðt1 ; t 2 Þ

ð2Þ

which quantifies the strength of oscillations in a yet to be specified manner (see Section 4.4) and where q ¼ u; p; T, with respect to a   control function / @ Xfs ; t defined on a boundary @ X/ and time interval ðt 1 ; t 2 Þ such that the constraints (the Boussinesq approximation of the Navier–Stokes and energy equations)

@u þ u  ru þ rp  r2 u þ GrTe ¼ 0; @t r  u ¼ 0; in X @T 1 þ u  rT  r2 T ¼ 0; in X @t Pr u ¼ 0; on ð@ X  @ Xfs Þ

in X

ð3Þ ð4Þ ð5Þ ð6Þ

u  n ¼ 0;

on @ Xfs  ðt1 ; t 2 Þ    T ðI  nn Þ  ru þ ðruÞT  n þ RerT ¼ 0;

ð7Þ on @ Xfs

W½T  ¼ 0; on ð@ X  @ Xfs Þ   rT  n  / @ Xfs ; t ¼ 0; on @ Xfs   /min 6 / @ Xfs ; t 6 /max

ð8Þ ð9Þ ð10Þ ð11Þ

are satisfied on the same temporal interval of ðt 1 ; t2 Þ. The duration t2  t1 of the time period over which the optimization problem is defined is assumed to be large compared to the period of the oscillation. In Eq. (3) e is a unit vector in the direction in which gravity acts (in the present work e ¼ ð0; 1Þ). The unit vector normal to the free surface is n. The function

W½T; C=2; y ¼ T W½T; þC=2; y ¼ T  1 W½T; x; 1=2 ¼

ð12Þ

@T @y

in Eq. (9) defines the boundary conditions for the energy equation Eq. (5) on the part of the boundary which is made by solid walls, with Eq. (10) defining them at the free surface where control is applied. In the absence of control the free surface is considered adiabatic. In the presence of control heat is supplied along the entire free surface creating a normal temperature gradient in the liquid. Such a heating could be realized, e.g., by laser irradiation (Kamotani and Ostrach, 1994; Benz et al., 1998). Eq. (11) is an inequality constraint on the control variables which quantifies the realizable limits of control input. The above equations have been made dimensionless using the 2

2

scales (see also Fig. 1) h; h =m; m=h; qm2 =h , DT ¼ T H  T C , for length, time, velocity, pressure and temperature, respectively. The resulting dimensionless groups are the Prandtl, Grashof and thermocapillary Reynolds numbers

m

Pr ¼ ; k

3

Gr ¼

gbDTh

m2

;

and Re ¼

cDTh : qm2

The Prandtl, Grashof, thermocapillary Reynolds numbers and aspect ratio of C (defined in Fig. 1) investigated in the present work are 10; 300; 2500 and 8, respectively. They are chosen based on the linear stability results of Kuhlmann and Albensoeder (2008) such that supercritical flow oscillations exist. An estimate of how supercritical the state is, is given by  ¼ ðRe  Rec Þ=Rec ¼ :28 using data from Fig. 3 of Kuhlmann and Albensoeder (2008) for C 6 7. 3. Numerical solution of the governing equations Eqs. (3)–(5) are discretized using a spatially second-order finite difference method. The temporal variation of convection is treated using the explicit third-order Adams–Bashforth scheme, while diffusion is treated implicitly using the Euler method to avoid the viscous stability constraint on the time step. The Marangoni boundary condition Eq. (8) and the normal temperature gradient boundary condition Eq. (10) at the free surface are treated using a temporally fully implicit method, meaning they are evaluated at the current time level only. The resulting linear systems are solved at each time step using the MA41 direct sparse matrix solver routine of the Harwell Subroutine Library (Numerical Analysis Group, 2007). Since the velocity and pressure are solved together as a

F.H. Muldoon, H.C. Kuhlmann / International Journal of Heat and Fluid Flow 56 (2015) 28–42

coupled system, the factorization error existing in projection methods does not exist. Solutions are obtained on a cartesian mesh, with grid stretching to increase resolution at the walls and free surface. Unless otherwise mentioned, all results were obtained using a 266  50 grid in the x and y-directions respectively. In the interests of brevity and to avoid duplication, we refer the reader to Muldoon (2013) for the details of the method used to discretize the governing equations and solve the resulting algebraic system. 4. Solution procedure for the optimization problem 4.1. Adjoint equations and gradient of objective function In the present work, a conjugate gradient method is used to find the minimum of the objective function and, therefore, the gradient of the objective function with respect to the control variables must be obtained. Note that the objective function Eq. (2) does not   explicitly depend on the control variable function / @ Xfs ; t , but does so only implicitly through the constraints defined by Eqs. (3)–(11). As a result, obtaining the gradient using a computationally efficient method is not a simple task. While a perturbation approach (Squire and Trapp, 1998) is simple to implement, such an approach is computationally efficient only if there exists a handful of control variables, since a perturbed objective function must be evaluated for each control variable. Since in the present work, there are up to 1:4  107 control variables (see Section 5.3.3), the inefficiency of this approach makes it unfeasible. A computationally efficient means of obtaining the gradient involves solving a set of adjoint equations which, in the present work, are derived from the discretized versions of Eqs. (3)–(10). The gradient used in the present work is obtained by this method which provides a means to obtain the gradient of the objective function with respect to the control variables, the computational cost, as measured by the number of arithmetical operations, of which is independent of the number of control variables. The means by which the adjoint equations are derived and solved in the present work are described in Muldoon (2013). Derivations of the adjoint equations can be found in Muldoon (2008), Giles et al. (2003), Nielsen (1998). Details of issues involved in deriving and solving the adjoint equations can be found in Nadarajah and Jameson (2000, 2001), Nielsen et al. (2009, 2004), and McNamara et al. (2004), and references therein. Two issues are worth mentioning, however. The first is that the adjoint equations have a form very similar to the mathematical model that they are derived from and, therefore, solution strategies can be very similar to those used in solving the mathematical model. The second is that the required data storage scales linearly with the temporal dimension of the control problem and can be a severe computational restriction. The mathematical model described, the adjoint equations derived from them and a conjugate gradient optimizer have been implemented in a Fortran 2003 computer program called Tetra, developed by the first author. All results in the present work have been obtained from this computer program. Additional details concerning Tetra can be found in Muldoon (2004). Validation of the results obtained from Tetra for the problem investigated in the present work by means of comparisons against another code, a linear stability analysis and an extensive study of grid and time step dependence may be found in Muldoon (2013).

in Eq. (16) and therefore would be expected to have a smaller and smaller effect on the objective function. For this reason, the present work investigates the effect of finding control over a portion (termed the temporal advancement length and denoted by tA where tA 6 t EH ) of the event horizon. This control, however, is determined such that it minimizes the objective function over the entire event horizon for objective functions using average temporal weighting Eq. (16) or at its end for the terminal objective functions Eq. (17). Such shrinking of the control variable space through removal of control variables to which the objective function is less sensitive should also aid solution of the resulting minimization problem by making it better conditioned. The total control from t1 to t 2 is then obtained by piecing together without any breaks the control over the temporal advancement length found within each event horizon. An illustration of this temporal piecing together of the control is given in Fig. 2 where control is initiated at t1 . An objective function defined from ½t1 ; t1 þ tEH Þ is minimized in a control variable space defined from ½t1 ; t1 þ tA Þ. Using these control variables, the flow is then advanced by solving Eqs. (3)–(10) over the interval ½t1 ; t1 þ t A Þ. A new objective function defined from ½t1 þ tA ; t1 þ tA þ tEH Þ is then minimized in a control variable space defined from ½t 1 þ t A ; t 1 þ 2tA Þ. Using these control variables, the flow is then advanced over the interval ½t1 þ tA ; t1 þ 2tA Þ. And so forth. In the continuous problem of Section 2.1 the control is the continuous function /½x; t. When the continuous problem is discretized according to Section 3 the continuous function /½x; t becomes the discrete temperature gradient normal to the free surface at each grid point at the free surface (excluding the two points on the left and right walls) at each time step. As a result the control variables are these discrete temperature gradients and possess spatial and temporal dimensions. The grid spacing varies such that the grid spacing normal to the boundaries is substantially smaller, this means that the influence of control variables in regions with small grid spacing generally would be smaller than those in regions with large grid spacing since the amount of heat added by the control is proportional to the grid spacing parallel to the free surface. 4.3. Optimization algorithm The optimization strategy in the present work uses the Hestenes–Stiefel conjugate gradient method (Nocedal and Wright, 1999) along with a one-dimensional line search to find the minimum along the search direction. The search direction is reset to the gradient direction every 10 iterations. The results and efficiency of a conjugate gradient method applied to a non-linear problem depend heavily on the termination criteria and the line search method and are described in detail in

4.2. Control variables The control variables are functions of space and time. As a result of the one-way nature of time, as the location of the control in the temporal dimension moves towards the end of the event horizon, the control affects a smaller and smaller portion of the integrand

31

Fig. 2. Schematic of the control scheme.

32

F.H. Muldoon, H.C. Kuhlmann / International Journal of Heat and Fluid Flow 56 (2015) 28–42

Muldoon (2013). At the beginning of each new event horizon the initial guess for the control variables is ð/max þ /min Þ=2. Assuming that it is desired to find optimal control from t 1 to t2 and N is such that Nt A 6 t 2  t 1 < ðN þ 1ÞtA , a flow chart of the optimization algorithm is given in Table 1.

the temporally varying spatially averaged fluctuating squared temperature !2 Z Z t1 þðm1ÞtA þtEH 1 1 0 0 T T ½t; m ¼ T ½x; y; t dt  T ½x; y; t  dA A½X n @ X Xn@ X t EH t1 þðm1ÞtA

ð14Þ

4.4. Objective functions Since the goal of the present work is to suppress the flow oscillations using an optimization algorithm it is necessary to quantify them in the form of a scalar objective function which will then be minimized. There is no unique way to do this, with any measure that captures the unsteadiness of the flow in some sense quantifying them. The present work investigates measures based on the fluctuating kinetic energy, the fluctuating temperature and the fluctuating pressure to do this. Ideally, the objective function used by the minimizer should be defined by Eq. (2). However, due to finite computer resources, this is impractical as the time period ½t 1 ; t 2  in Eq. (2) is quite large and therefore evaluation of this equation, in particular the storage requirements of the adjoint equations, are too computationally expensive. Therefore, the objective function can be defined only over a much smaller time period. This finite period of time, commonly called the event horizon and abbreviated in the present work as tEH , is the distance forward in time over which the flow is evolved and in which new control variables are generated to minimize the objective function. It is therefore the length of time over which the gradient of the objective function with respect to the control variables must be computed. Clearly, the longer the event horizon, the more control exists over the objective function. However, this must be balanced against the increased computational expense of longer event horizons. An important computational limitation is that the velocity and temperature must be stored at each time step in the event horizon for use in solving the adjoint equations, the solution of which is needed to determine the gradient. 4.4.1. Spatial measures The spatial measures used in defining the objective functions are twice the temporally varying spatially averaged fluctuating kinetic energy, 0  R 2 1 t 1 þðm1Þt A þt EH 1 Z u½x; y; tdt  u½x; y; t  1 B tEH t1 þðm1ÞtA C k½t; m ¼ @  R 2 AdA A½X n @ X Xn@ X t1 þðm1ÞtA þtEH 1 þ tEH t1 þðm1ÞtA v ½x; y; tdt  v ½x; y; t

and the temporally varying spatially averaged fluctuating squared pressure !2 Z Z t1 þðm1ÞtA þtEH 1 1 0 0 p p ½t; m ¼ p½x; y; tdt  p½x; y; t dA: A½X n @ X Xn@ X tEH t1 þðm1ÞtA

ð15Þ where A has dimensions of area and X n @ X indicates the interior of the domain (excluding the boundary). In Eqs. (13)–(15), m ¼ 1; 2 . . . indicates the discrete event horizons and corresponding advancement lengths and therefore k½t; m; T 0 T 0 ½t; m and p0 p0 ½t; m are specific to a unique event horizon. t A is the length of time for which control is sought within each event horizon (see Section 4.2), while t 1 is the starting time for the global control problem (see Section 2.1). The temporal integral of a generic variable q½t

Z

t 1 þðm1Þt A þt EH

q½tdt

t 1 þðm1Þt A

is to be understood as carried out over the individual event horizon defined by the time period from t 1 þ ðm  1Þt A to t1 þ ðmÞt A . Note that the spatial integrals in these equations are only over the interior of the domain, therefore the velocity, temperature and pressure on the free surface are not included. 4.4.2. Temporal weighting In the previous section spatial measures were defined that depend on time and the individual event horizons (numbered by m). In order to obtain a scalar objective function within each unique event horizon, it is necessary to remove the temporal dependence of these spatial measures. Following from this the values of the spatial measures within the event horizons should be weighted in some manner in a temporal integral over the event horizon. In the present work two approaches are used. The first is a constant weighting of unity, which yields the average over the event horizon, the second is to consider only the measure at the final (terminal) time step of the event horizon. Applying the two temporal weightings

ð13Þ

q¼ Table 1 Optimization flow chart.

1 t EH

Z

t1 þðm1Þt A þt EH

q½t; mdt

ð16Þ

t 1 þðm1Þt A

qterminal ¼ q½t1 þ ðm  1Þt A þ t EH ; m

ð17Þ

For m ¼ 1; M 1 2 3 4 5 6 7 8 9 10

Make a guess for the control variables for t1 þ ðm  1ÞtA < t 6 t1 þ mtA Solve Eqs. (3)–(9) and store u and T from t 1 þ ðm  1Þt A to t1 þ ðm  1ÞtA þ tEH Using u; T from 2, solve adjoint eqs. from t 1 þ ðm  1Þt A þ tEH to t1 þ ðm  1ÞtA Using adjoint variables from 3, compute the gradient of the objective function Using the gradient, get a search direction from a conjugate gradient method Do a line search in the search direction and find a minimum along the line Set the control variables to that determined by 6 If the value of the objective function at the minimum along the line found in 6 satisfies the stopping criteria, go to 9, if not return to 2 Using control found in 6, advance the solution from t1 þ ðm  1ÞtA to max½t1 þ mtA ; t 2  Next m

where q stands for one of the three spatial measures k; T 0 T 0 and p0 p0 yields six objective functions studied in the present work. Note that the terminal objective functions mean that the optimizer is free to increase the spatial measures Eqs. (13)–(15) within each event horizon as long as so doing results in a lower value of the measure at the final time step of the event horizon. 4.5. Constraints For the optimization problem to be physically realistic, constraints must be imposed on the control variables. The control in the present work is the normal temperature gradient at the free surface which leads to the constraint that it be greater than or equal to zero, i.e., /min ¼ 0, precluding cooling of the free surface which would be difficult to implement in a physical system. In

F.H. Muldoon, H.C. Kuhlmann / International Journal of Heat and Fluid Flow 56 (2015) 28–42

addition, an upper bound of /max ¼ 0:1 on the normal temperature gradient at the free surface is imposed. These inequality constraints are implemented by the inclusion of a penalty term in the objective function used by the optimizer. This is accomplished using the penalty function given by

0 !2105 1 n N I 1 X X /  ð / þ / Þ max min i 2 @ A; Hp ¼ 1 ð/max  /min Þ 2 n¼1 i¼2

ð18Þ

where the inner summation over i ¼ 2 . . . I is from first to the last grid points on the free surface excluding the grid points lying on the vertical walls. This function has the desired property of being very close to zero for all /min < /i < /max and increasing extremely rapidly if any /i > /max or /i < /min . The outer summation on the time index n is present, since the penalty is present at each time step of the event horizon. The objective function used by the conjugate gradient minimizer H is the sum of the desired objective function (Eqs. (16) and (17)) and the penalty term Hp Eq. (18). As the sole purpose of the penalty function is to ensure that the inequality constraints are satisfied, it is desired that the penalty term does not introduce any significant changes to the objective function in the regions in control variable space in which solutions are found, which would thereby ‘‘corrupt’’ the desired objective function (Eqs. (16) and (17)). These is indeed the case in the present work as in all results, the size of the penalty term is negligible compared to the desired objective function when the optimizer has converged. 4.6. Sensitivity of the objective function An important consideration is the behavior of the gradient of the objective function with respect to the control variables. The control variable / varies in both space and time. In the present work, the control variables /i are spatially discretized while the control variables are free to vary in time within the event horizon

Fig. 3. L1 spatial norm of the gradient of the objective function divided by the value of the objective function (vertical axis) versus time (horizontal axis). Temporal extent of the gradient is t ¼ ½45; 65. The region from t ¼ 45 to t ¼ 55 is not shown as it is nearly identical to the region t ¼ ½55  60. Gradient is at / ¼ 0 (i.e., zero normal temperature gradient at the free surface) in control variable space. No penalty term is used in defining the objective functions. 134  26 grid. The green vertical dashed lines indicate the lengths of different temporal extents measured from t ¼ 65 corresponding to the lengths of the event horizons used in the present work. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

33

length. The temporal behavior of the gradient of the objective function can be investigated by applying a norm to the spatial distribution of the control. Fig. 3 shows the spatial L1 norm of the gradient of the objective functions normalized by the value of the objective function as a function of time. Due to the inner integrals over time in Eqs. (13)–(15), the gradient of the objective functions even over the same period of time differs due to the differing lengths of the event horizons. This may be seen in Fig. 4. For all objective functions defined using temporal averaging over the event horizon (i.e., using Eq. (16)), this norm becomes quite small towards the end of the event horizon. This is because of the one-way nature of time, as a result of which, as the distance in time of the control from the end of the event horizon vanishes, the control affects only a vanishingly small portion of the integrand in Eq. (16). The general trend is that the sensitivity increases as the temporal distance from the end of the event horizon (at t ¼ 65) increases until the maximum sensitivity is reached  2 units of time from the end of the event horizon, after which the sensitivity decays slightly and then reaches a periodic state. This behavior of the sensitivity suggests that, as the length of the event horizon increases, the control effectiveness will increase until the peak sensitivity lies within the event horizon (i.e., at 2 units of time), with further increases in the length of the event horizon unlikely to improve the effectiveness of the control in reducing the objective function. The use of terminal temporal weighting Eq. (17) results in a very different temporal distribution of the sensitivity. For these objective functions the sensitivity does not go to zero as the distance in time of the control from the end of the event horizon vanishes. There however exist fundamental differences between the objective function based on the pressure and those based on the fluctuating kinetic energy and temperature. For p0 p0terminal there is a sharp spike in the sensitivity at the end of the event horizon, with the sensitivity there being two orders of magnitude larger than that for t ¼ ½55;  60. In contrast for T 0 T 0terminal the peak sensitivity occurs at t  0:7, while for kterminal the peak sensitivity occurs at t  0:08. For both kterminal ; T 0 T 0terminal the sensitivity reaches a periodic state with the sensitivity not more than an order of magnitude smaller than the peak sensitivity. Fig. 5a shows the Fourier spectrum over t ¼ ½45; 55 of the sensitivity. For this case the fundamental frequency of the flow is 2:4 (Table 3 of Muldoon (2013)). For all objective functions the fundamental frequency of the flow given by the flow oscillations is by far the dominant one, with the other frequencies, with orders of magnitude smaller amplitudes, appearing strictly as multiples of the fundamental frequency, with the exception of the p0 p0terminal objective function where some incommensurate frequencies appear. An explanation for the strong changes in sensitivity near the end of the event horizon can be understood by recognizing that the mathematical model contains two effects possessing different time scales. One is an immediate (elliptic) effect arising from the assumption of an incompressible fluid and hence infinite speed of sound. As a result, control input at a given time step will influence the solution within the entire spatial domain at that and all future time steps. Due to the absence in the mathematical model of a temporal derivative of pressure, pressure at previous time steps has no effect on the solution at future time steps, with pressure playing the role of a Lagrange multiplier which enforces the solenoidal constraint at each time step. However, the pressure at every time in the event horizon enters p0 p0terminal through the inner integral of Eq. (15). Due to this the terminal objective functions are in a sense not truly terminal as would be the case in the absence of the inner integral.

34

F.H. Muldoon, H.C. Kuhlmann / International Journal of Heat and Fluid Flow 56 (2015) 28–42

Fig. 4. L1 spatial norm of the gradient of the objective function divided by the value of the objective function (vertical axis) versus time (horizontal axis) for various event horizon lengths. All event horizons end at t ¼ 65. Gradient is at / ¼ 0 (i.e., zero normal temperature gradient at the free surface) in control variable space. No penalty term is used in defining the objective functions. 134  26 grid.

lengths of 0:5t EH ; 0:75t EH and t EH , giving a total of 72 unique optimization problems. The maximum event horizon length is chosen to be sufficient to capture the peak sensitivities (see Section 3) of the objective function. The instantaneous heat flux through the cold wall

qCOLD ½t ¼

Z

1=2

1=2

@T ½C=2; y; t dy @x

ð19Þ

is used to quantify the effect of the control on the flow oscillations. 5.1. Monitoring the heat flux

Fig. 5. Power density spectrum of the L1 spatial norm of the gradient of the objective function using an event horizon length of 20, at t ¼ 45. Gradient is at / ¼ 0 (i.e., zero normal temperature gradient at the free surface) in control variable space. No penalty term is used in defining the objective function. k is defined as 1=period rather than as 2p=period. The natural frequency for this case is 2:400 (Table 3 of Muldoon (2013)). 134  26 grid.

5. Results The optimization procedure used in the present work has four parameters which define a specific optimization problem, these being the three spatial measures (Section 4.4.1) the two types of temporal weighting (Section 4.4.2) and the event horizon ðt EH Þ and temporal advancement ðtA Þ lengths. In the present work, four event horizon lengths of 0:08; 0:24, 0:72 and 2:16 are used. Each event horizon length is used with three temporal advancement

Figs. 6–9 present the results of the heat flux through the cold wall for the objective functions used in the present work. It is evident that the results strongly depend on the spatial measures, the type of temporal weighting and the event horizon and temporal advancement lengths. In some cases, the flow oscillations are actually strengthened, while in others they are almost completely suppressed. For the event horizon length of 0:08 in no cases are the flow oscillations well suppressed. Results are quite different for the event horizon length of 0:24 where suppression occurs for the five cases ðk; tA ¼ 0:5t EH ; 0:75t EH Þ, ðT 0 T 0 ; tA ¼ 0:5tEH ; 0:75tEH Þ; ðT 0 T 0terminal ; tA ¼ 0:5t EH Þ. For the event horizon length of 0:72 suppression occurs for ðk; tA ¼ 0:5t EH ; 0:75t EH Þ and ðp0 p0 ; t A ¼ 0:5tEH ; 0:75tEH Þ. Noticeable is the complete lack of suppression for all objective functions based on the temperature for this event horizon length. Indeed for these cases, the flow oscillations are sometimes strengthened. For the event horizon length of 2:16 significant suppression (although to a lesser degree than the previously listed   cases) occurs for ðk; t A ¼ 0:5t EH ; 0:75tEH Þ, p0 p0 ; tA ¼ 0:5t EH ; 0:75t EH and ðT 0 T 0 ; t A ¼ 0:75tEH Þ. At this event horizon length, terminal temporal weighting does not lead to suppression for any objective functions. As a general rule, terminal weighting leads to poor suppression of the flow oscillations, with there being only one case  0 0  T T terminal ; t A ¼ 0:5tEH where the oscillations are almost

F.H. Muldoon, H.C. Kuhlmann / International Journal of Heat and Fluid Flow 56 (2015) 28–42

35

Fig. 6. Heat flux through the cold wall Eq. (19). Event horizon length (tEH ) of 0:08.

completely suppressed. Effective suppression of the flow oscillations using objective functions based on pressure occurs only with the two largest event horizons. As poor suppression of the flow oscillations is obtained using the short event horizon length of 0:08, we focus on results using the other event horizon lengths.

be to use a stochastic method to generate searches from different starting guesses. The essential problem is that while increasing the event horizon length provides more ability to control the flow, the increased event horizon length greatly increases the difficulty of the minimization problem solved in each event horizon, making it more likely that the optimal control cannot be found.

5.2. Effect of the control parameters The fact that the greatest damping of the flow oscillations does not occur using the longest event horizon requires an explanation. As a general rule, the longer the event horizon, the more control exists over the flow. Aside from this, it certainly is true that control found using a given event horizon length should (in theory) also be able to be found within a longer event horizon. This counter-intuitive behavior can be explained by the fact that the objective function depends on the solution to a complicated non-linear system of equations (i.e., the Navier–Stokes and energy equations). Following from this, the objective function can be a complicated function of the control variable space and possess multiple local minima. As the event horizon increases, and therefore linearly increases the dimension of the control variable space, the likelihood of the objective function possessing multiple local minima increases. The result is that the likelihood of a gradient-based method, as used in the present work, becoming trapped in a local minimum increases. A gradient-based method will find a local minimum quickly but the specific one it finds will depend on the starting guess and the details of the (deterministic) search method. A potential, although computationally expensive, means of avoiding becoming trapped in a local minimum, would

5.3. Cases with significant suppression Attention is now focused on the cases for which the flow oscillations were best suppressed. There are six cases for which significant suppression is achieved. These are listed in Table 2. The spatial L1 norm of the control as a function of time for these six cases is shown in Fig. 10. Immediately evident is the only slight temporal variation of this norm for ðk; t EH ¼ 0:24; tA ¼ 0:75t EH Þ compared to all other cases. The cases ðk; t EH ¼ 0:24; t A ¼ 0:5t EH Þ and ðk; t EH ¼ 0:72; tA ¼ 0:5t EH Þ have a much stronger variation, however, the damping of the flow oscillations is nearly the same as for ðk; t EH ¼ 0:24; tA ¼ 0:75t EH Þ. All plots show significant variation between each other, indicating that there are a number of control inputs which will significantly suppress the flow oscillations. For all cases, there is a transient period during which the flow oscillations present in the uncontrolled case are damped. After this occurs, the norm assumes a quasi-periodic behavior, during which suppression of the flow oscillations is maintained. This suppression, however, requires a constant input of heat flux in the form of a temporally and spatially varying normal temperature gradient at the free surface. Fig. 11 gives the Fourier spectrum of Fig. 10. It

36

F.H. Muldoon, H.C. Kuhlmann / International Journal of Heat and Fluid Flow 56 (2015) 28–42

Fig. 7. Heat flux through the cold wall Eq. (19). Event horizon length (t EH ) of 0:24.

can be seen that for all cases, the dominant frequency is that of the temporal advancement length, indicating the dependence of the computed control on this parameter. From this it can be inferred that an improved optimization/control technique could find better control which depends more on the physics of the problem and less on the parameters governing the optimization/control technique. Movies 2 and 3 show the control and the resulting flow field for the six best cases in Table 2 along with the flow with a constant normal temperature gradient at the free surface of / ¼ ð/max þ /min Þ=2 and with the absence of control (/ ¼ 0). From this movie it is evident that suppression of the flow oscillations can be accomplished with control inputs that differ significantly from each other. It is also clear that suppression does not occur due to the ‘‘brute force’’ addition of a constant heat flux coming from / ¼ ð/max þ /min Þ=2 alone, but rather from imposing the normal temperature gradient in certain spatially and temporally varying manners. For all cases, the control fluctuates around ð/max þ /min Þ=2. This fluctuation is very small near the walls, indicating that the addition of heat near the walls has little effect on the flow. This corresponds with the observation in Muldoon (2013) that the details of the flow in the region immediately adjacent to the cold wall are unimportant for the rest of the flow which is plausible, since the proximity of the wall represents a damping by diffusion of momentum or heat of any driving force. For all cases, the fluctuation of the control is most pronounced around x ¼ 2. This is near the region where the oscillatory instability is supposed to be triggered (Peltier and Biringen, 1993).

Aside from suppression of the flow oscillations, the controlled flow differs from the uncontrolled flow in other respects. One is that the free surface velocity for the six best controlled cases is always negative, i.e., from the hot to the cold wall, while the uncontrolled and the case with constant control of ð/max þ /min Þ=2 have small regions with positive velocities for 2 < x < 3. This is evident in Fig. 12. In this figure can be seen regions where the temperature gradient tangent to the free surface is negative and hence the local thermocapillary force is directed opposite the overall hot-to-cold wall thermocapillary force. However, this force, instead of causing an instability as might be expected, instead suppresses one. Another difference is that the controlled flow has five vortices while the uncontrolled flow has four (this can be seen most clearly in Movie 3). However, as the case of constant control of / ¼ ð/max þ /min Þ=2 also exhibits five vortices, the additional vortex appears to arise as an overall effect of heat addition at the free surface. For the case k; tEH ¼ 0:24; tA ¼ 0:75tEH , the temporal variation in the control needed to keep the flow oscillations suppressed is extremely small and a reasonable question to ask is whether the flow has somehow reached a steady state solution that is stable. That this is not the case can be seen in Fig. 13, where the control is turned off after the flow oscillations have been suppressed. The result is that the flow oscillations return in the same form as before the control was turned on. The path by which they return is almost identical to that for the case k; tEH ¼ 0:72; tA ¼ 0:75tEH where the temporal variation in the

F.H. Muldoon, H.C. Kuhlmann / International Journal of Heat and Fluid Flow 56 (2015) 28–42

control needed to keep the flow oscillations suppressed is much larger. 5.3.1. Behavior of objective functions The behavior of the objective function (i.e., how it changes when the control variables change) is of interest, as it can indicate whether the optimization algorithm is appropriate. In particular, the behavior of the objective function along the search direction in control variable space is of great interest. Fig. 14 shows the change in the objective function when moving along the line in control variable space defined by the search direction of the conjugate gradient algorithm for the first iteration within an event horizon near t ¼ 60. There are between eight and forty-five sample points, beginning at the point in control variable space where the search direction is defined and terminating at a distance equal to the starting length used by the line search. Note that the distance a along the search direction is normalized by its maximum value determined from the inequality constraints on the control variables, while the objective function is normalized by its value at the beginning of the line search (i.e., where a ¼ 0). The reason that in some cases data are not plotted up to the maximum value of a ¼ 1 is that the implementation of the optimization algorithm has logic which, depending on the course of the minimization process, speeds up the line search by reducing the distance along the search direction that the search is carried out over. Note that the desired objective function is plotted. Due to the large exponent in Eq. (18), the difference between the desired

37

objective function and the objective function plus the penalty term is minuscule wherever the inequality bounds on the control variables are satisfied. Plotting the desired objective function avoids the irrelevant sharp spike in the penalty term in Eq. (18), which has a value greater than one at the termination of the line, where at least one of the control variables equals the value of the inequality bound on the control variables. The first two rows (Fig. 14a–f) represent the six cases for which best results where obtained (Table 2). It can be seen that for all the cases (except for T 0 T 0terminal , tEH ¼ 0:24; t A ¼ 0:5tEH ), a single local minimum exists between a ¼ 0 and the end of the search direction defined by the inequality constraints. For the cases for which poor results are obtained, the most common behavior is that shown in Fig. 14i–j where the objective function monotonically decays along the search direction until it reaches the maximum allowable value of the search vector determined by the inequality constraints on the control variables. An exception is the longer event horizon of tEH ¼ 2:16 (Fig. 14g–h). where in many cases there is one local minimum far from the maximum allowable value of the search vector. In all cases examined, no more than two local minima were found. As the iterations of the conjugate gradient minimizer proceed within each event horizon, the behavior of the objective function along the search line changes, with the rate of change of the objective function along the search direction greatly decreasing for most cases (Fig. 15g–h). This is accompanied with the minimum almost always being found at the end (determined by the inequality constraints, Eq. (11)) of the search vector. This indicates that further

Fig. 8. Heat flux through the cold wall Eq. (19). Event horizon length (t EH ) of 0:72.

38

F.H. Muldoon, H.C. Kuhlmann / International Journal of Heat and Fluid Flow 56 (2015) 28–42

Fig. 9. Heat flux through the cold wall Eq. (19). Event horizon length (tEH ) of 2:160. Table 2 Parameters for the best controlled cases. k

tEH ¼ 0:24; tA ¼ :50tEH

k

tEH ¼ 0:24; tA ¼ :75tEH

k kterminal

tEH ¼ 0:72; tA ¼ :50tEH

T0T0 T 0 T 0terminal

tEH ¼ 0:24; tA ¼ :50tEH tEH ¼ 0:24; tA ¼ :75tEH tEH ¼ 0:24; tA ¼ :50tEH

Fig. 11. Power density spectrum of data in Fig. 10 from t ¼ 57:5 to t ¼ 70. The order of the cases corresponds to that in Table 2, i.e., (a) – k; t EH ¼ 0:24; t A ¼ :50tEH , (b) – k; t EH ¼ 0:24; t A ¼ :75tEH , etc.

decrease of the objective function is prevented by the inequality constraints.

Fig. 10. Spatial L1 norm of the control variables as a function of time. The order of the cases corresponds to that in Eq. (2), i.e., (a) – k; t EH ¼ 0:24; t A ¼ :50tEH , (b) – k; tEH ¼ 0:24; t A ¼ :75tEH , etc.

5.3.2. Performance of the minimizer Fig. 16, shows the change in the normalized objective function along with the maximum and minimum of the control variables relative to the iterations of the minimizer within a single event horizon. For all cases, the graph of the objective function sharply

F.H. Muldoon, H.C. Kuhlmann / International Journal of Heat and Fluid Flow 56 (2015) 28–42

39

Fig. 12. Graphical description of quantities at the free surface with y-axis indicating time and x-axis the position along the free surface. (I) Control. (II) Black indicates regions where @T=@x < 0, white otherwise. (III) u with black regions (within the white oval) indicating where u½y ¼ 1=2 > 0.

decreases at first, before leveling out and eventually terminating due to the termination criteria described in Muldoon (2013). It can be seen that many fewer iterations are needed for shorter event horizons before the minimizer stops according to the termination criteria. Also evident is the fact that both the minimum and maximum inequality constraints on the control variables are active in approximately equal measure. It can be seen that termination of the minimizer is strongly correlated with the inequality constraints becoming active. While not shown, examination of the

equivalent data shown in Fig. 16 for the objective functions based on kterminal ; T 0 T 0 , T 0 T 0terminal ; p0 p0 and p0 p0terminal has shown that in all cases the inequality constraints become active as the iterations proceed.

5.3.3. Computational requirements In solving non-linear optimization problems by use of gradient-based methods, a major portion of computational time

40

F.H. Muldoon, H.C. Kuhlmann / International Journal of Heat and Fluid Flow 56 (2015) 28–42

Fig. 13. Graphical description of the spatial and temporal development of u at the free surface after control is turned off at t ¼ 70. y-axis indicating time and x-axis the position along the free surface. Black regions indicate where u > 0.

Fig. 14. Objective function along search direction for the first iteration for an event horizon starting at t  60. Thin continuous lines are the objective functions evaluated at 45 locations along the search direction. Squares indicate objective function evaluations during the line search, while the solid triangle indicates the termination point of the line search. The y-axis is the objective function without the penalty term normalized by the value of the objective function at the start of the vector defining the direction of the line search, while the x-axis is the distance in control variable space along the search direction normalized by its maximum value determined from the inequality constraints on the control variables.

is often spent determining the minimum of the objective function along a search direction. The largest (smallest) ratio of the number of objective-function evaluations in the line search to the total number for the optimization problems solved in the present work is 0:85 (0:75). The largest (smallest) number of evaluations of the objective function for the simulations was 1:86  104 (3:07  102 ). The count of these objective function evaluations includes those at the origin of the search direction, which are computed when computing the gradient used to determine the search direction, plus the additional ones along the search direction

Fig. 15. Objective function along search direction for the last iteration for an event horizon starting at t  60. Thin continuous lines are the objective functions evaluated at up to 45 locations along the search direction. Squares indicate objective function evaluations during the line search, while the solid triangle indicates the termination point of the line search. The y-axis is the objective function without the penalty term normalized by the value of the objective function at the start of the vector defining the direction of the line search, while the x-axis is the distance in control variable space along the search direction normalized by its maximum value determined from the inequality constraints on the control variables.

needed to bracket the minimum during the line search. Also included are the objective function evaluations needed in determining the shape of the objective function along the search direction (Section 5.3.1). As each evaluation requires a solution of the discretized versions of Eqs. (3)–(10) over the event horizon, this involves a substantial amount of computational work. The maximum number of iterations used by the golden section line search was 15, the maximum number of iterations of the minimizer within any single event horizon was 23. The largest (smallest) ratio of the time spent solving the linear system arising from the adjoint equations to the total execution time was 0:20 (0:14), by way of comparison, the ratio of the time spent solving the linear system arising from the discretization of the flow and temperature equations to the total execution time was 0:77 (0:67). The non-linear optimization problems solved within each event horizon in the present work are computationally challenging in their own right without considering the necessity of solving them for numerous event horizons to determine the control variables for the entire temporal domain of interest. The dimension of the control variable space in which the minimizer operates in the present work is very large; for the 266  50 grid using the longest temporal advancement length, there are 264 (spatial dimension)  54,000 (temporal dimension) = 1:4  107 control variables in each event horizon. The calculations in the present work have been carried out using a single core of a 2.67 GHz Intel Xeon X5550 processor with a 8192 KB cache. The solution of the governing equations at each time step took 0.01 seconds of wall clock time on the 266  50 grid. The largest (smallest) amount of wall clock time required to solve the optimization problems was 654:8 (49:3) hours; by comparison, running the simulation without any control over the same time period required 2.75 hours. The memory requirements, due to the necessity of storing the velocity for use in solving the adjoint equations, scale linearly with the length of the event horizon. The largest event horizon used in solving the optimization problems required 17.1 GB of memory, the largest used in Fig. 3 required 20.6 GB. All simulations were carried out in double precision.

F.H. Muldoon, H.C. Kuhlmann / International Journal of Heat and Fluid Flow 56 (2015) 28–42

41

Fig. 16. Convergence history of the objective function along with the maximum and minimum of the control variables. Objective function normalized by its value at a ¼ 0. Final event horizon. H ¼ k.

6. Conclusions The present work has investigated issues involved with finding control in the form of normal temperature gradient at the free surface capable of suppressing flow oscillations in a flow driven by buoyancy and thermocapillary forces and modeled by the unsteady Navier–Stokes and energy equations under the Boussinesq approximation, within the context of a gradient-based optimization algorithm. Attention has been focused on the choice of objective function, which includes the event horizon and temporal advancement lengths in its definition. The control has been found by solving an optimization problem in a very high dimensional space with up to 1:4  107 control variables. The ability to suppress the flow oscillations is found to be very sensitive to the event horizon and temporal advancement lengths, with short event horizon and temporal advancement lengths greatly inhibiting their suppression. It has been found that using objective functions defined by a spatial measure temporally averaged over the event horizon generally provides better results than objective functions defined by a spatial measure at the final terminal time level of the event horizon. Basing the objective function on the fluctuating kinetic energy has been found to give better results than basing it on the fluctuating temperature which in turn generally gives better results than one based on the fluctuating pressure. However, objective functions based on fluctuating temperature and pressure are also capable of resulting in control that greatly suppresses the flow oscillations. Examination of the behavior of the objective function when moving in control variable space has shown that a gradient-based optimization method is appropriate for this flow as the number of local minima along the search direction is small. All results have been obtained using a single processor indicating that with the use of parallel computing, it is feasible to apply the method to investigate three-dimensional thermocapillary flows. It should be mentioned that a defect in the present work is that we only consider the liquid phase and therefore can only indirectly predict the affect of the application of a variable normal temperature gradient at the free-surface on the formation of dopant striations. In our opinion the ideal way to address this issue would be to

use as governing equations those that model the crystal growth process to a level that captures dopant transfer and striations and to formulate an objective function that quantifies dopant striations. The main difficulty would be the greatly increased computational cost due to the additional equations and more importantly the introduction of different spatial and temporal scales. This, however, is a direction that we would like to investigate in the future. Acknowledgements The first author is grateful for help concerning animations and visualization provided by Michael Saunders. The second author acknowledges support from the Austrian Space Agency under Project Number 819714. The vast majority of the computations were carried out on the cluster Phoenix and on the CAE system of Zentraler Informatikdienst of Vienna University of Technology. References Ben Hadid, H., Roux, B., 1990. Thermocapillary convection in long horizontal layers of low-Prandtl-number melts subject to a horizontal temperature gradient. J. Fluid Mech. 221, 77–103. Benz, S., Hintz, P., Riley, R.J., Neitzel, G.P., 1998. Instability of thermocapillary– buoyancy convection in shallow layers. Part 2. Suppression of hydrothermal waves. J. Fluid Mech. 359, 165–180. Davis, L., 1987. Genetic Algorithms and Simulated Annealing. Morgan Kaufmann Publishers Inc., San Francisco, CA. Giles, M.B., Duta, M.C., Muller, J.-D., Pierce, N.A., 2003. Algorithm developments for discrete adjoint methods. AIAA J. 41 (2), 198–205. Hurle, D.T.J. (Ed.), 1994. Handbook of Crystal Growth. North Holland. Kamotani, Y., Ostrach, S., 1994. Analysis of velocity data taken in surface tension driven convection experiment in microgravity. Phys. Fluids 6, 3601–3609. Kuhlmann, H.C., 1999. Thermocapillary convection in models of crystal growth. Springer Tracts in Modern Physics, vol. 152. Springer, Berlin, Heidelberg. Kuhlmann, H., Albensoeder, S., 2008. Three-dimensional flow instabilities in a thermocapillary-driven cavity. Phys. Rev. E 77, 036303-1–036303-15. Laure, P., Roux, B., Ben Hadid, H., 1990. Nonlinear study of the flow in a long rectangular cavity subjected to thermocapillary effect. Phys. Fluids A 2, 516– 524. McNamara, A., Treuille, A., Popovic, Z., Stam, J., 2004. Fluid control using the adjoint method. In: ACM Transactions on Graphics (ACM SIGGRAPH 2004). . Muldoon, F., 2004. Numerical Methods for the Unsteady Incompressible Navier– Stokes Equations and Their Application to the Direct Numerical Simulation of Turbulent Flows. Ph.D. thesis, Louisiana State University. .

42

F.H. Muldoon, H.C. Kuhlmann / International Journal of Heat and Fluid Flow 56 (2015) 28–42

Muldoon, F., 2013. Control of hydrothermal waves in a thermocapillary flow using a gradient-based control strategy. Int. J. Numer. Methods Fluids 72 (1), 90–118, http://dx.doi.org/10.1002/fld.3735. Muldoon, F., 2008. Control of a simplified unsteady film-cooling flow using gradient-based optimization. AIAA J. 46 (10), 2443–2458 (October). Müller, G., Ostrogorsky, A., 1994. Convection in melt growth. In: Hurle, D.T.J. (Ed.), Handbook of Crystal Growth, vol. 2b. North Holland, p. 709 (Bulk Crystal Growth). Nadarajah, S., Jameson, A., 2000. A comparison of the continuous and discrete adjoint approach to automatic aerodynamic optimization. In: AIAA 38th Aerospace Sciences Meeting and Exhibit, Reno Nevada, AIAA Paper 2000-0667. Nadarajah, S., Jameson, A., 2001. Studies of the continuous and discrete adjoint approaches to viscous automatic aerodynamic shape optimization. In: 15th Computational Fluid Dynamics Conference, Anaheim, CA, AIAA-2001-2530. Nielsen, E., 1998. Aerodynamic Design Sensitivities on an Unstructured Mesh using the Navier–Stokes Equations and a Discrete Adjoint Formulation. Ph.D. thesis, Virginia Polytechnic Institute and State University. Nielsen, E.J., Lu, J., Park, M.A., Darmofal, D.L., 2004. An implicit, exact dual adjoint solution method for turbulent flows on unstructured grids. Comput. Fluids 33 (9), 1131–1155, . Nielsen, E.J., Diskin, B., Yamaleev, N.K., June 2009. Discrete adjoint-based design optimization of unsteady turbulent flows on dynamic unstructured grids. In: AIAA-2009-3802.

Nocedal, J., Wright, S.J., 1999. Numerical Optimization. Springer-Verlag New York, Inc., 175 Fifth Avenue, New York, NY 10010, USA. Numerical Analysis Group, 2007. Harwell Subroutine Library, A Collection of Fortran Codes for Large Scale Scientific Computation. . Peltier, L.J., Biringen, S., 1993. Time-dependent thermocapillary convection in a rectangular cavity: numerical results for a moderate Prandtl number fluid. J. Fluid Mech. 257, 339–357. Shevtsova, V.M., Legros, J.C., 2003. Instability in thin layer of liquid confined between rigid walls at different temperatures. Acta Astronaut. 52, 541–549. Shiomi, J., Amberg, G., 2002. Active control of a global thermocapillary instability. Phys. Fluids 14, 3039–3045. Shiomi, J., Amberg, G., 2005. Numerical investigation of feedback control of thermocapillary instability. Phys. Fluids 17, 054107-1–054107-12. Shiomi, J., Amberg, G., Alfredsson, H., 2001. Active control of oscillatory thermocapillary convection. Phys. Rev. E 64, 031205-1–031205-7. Shiomi, J., Kudo, M., Ueno, I., Kawamura, H., Amberg, G., 2003. Feedback control of oscillatory thermocapillary convection in a half-zone liquid bridge. J. Fluid Mech. 496, 193–211. Squire, W., Trapp, G., 1998. Using complex variables to estimate derivatives of real functions. SIAM Rev. 40 (1), 110–112. Xu, J., Zebib, A., 1998. Oscillatory two- and three-dimensional thermocapillary convection. J. Fluid Mech. 364, 187–209.