Bounds on energy dissipation for a stress-driven shear flow

As in the previous example, consider an idealized two-dimensional layer of fluid, periodic in the horizontal (\(x\)) direction with period \(\Lambda\), bounded below (\(y=0\)) by a solid wall, and driven at the surface (\(y=1\)) by a shear stress of non-dimensional magnitude \(G\) (known as the Grashoff number).

Suppose that we can find a background field \(\varphi(y)\) that satisfies the boundary conditions

\[\begin{aligned} \varphi(0)&=0, & \varphi'(1)&=G. \end{aligned}\]

and such that the integral constraint

\[\begin{split}\begin{aligned} \int_0^1 &\left\{ \frac{1}{k_n^2}\left[ \vert u''(y) \vert^2 + \vert v''(y) \vert^2\right] +2 \left[ \vert u'(y) \vert^2 + \vert v'(y) \vert^2\right] \right. \\ &\qquad\qquad\qquad \left. + k_n^2 \left[ \vert u(y) \vert^2 + \vert w(y) \vert^2\right] - \frac{2}{k_n}\,\varphi'(y) \left[ u(y)\,v'(y) - v(y)\,u'(y) \right] \right\} {\rm d}y \geq 0 \end{aligned}\end{split}\]

holds for all wavenumbers \(k_n\) and all functions \(u(v)\) and \(v(y)\) that satisfy the boundary conditions

\[\begin{split}\begin{aligned} u(0) &= 0, & u(1) &= 0, & u'(0)&=0, & u''(1) &= 0,\\ v(0) &= 0, & v(1) &= 0, & v'(0)&=0, & v''(1) &= 0.\\ \end{aligned}\end{split}\]

Then, the time-averaged bulk energy dissipation coefficient \(C_\varepsilon\) can be bounded according to

\[C_\varepsilon \leq G \times \left[ 2\varphi(1) - \int_0^1 \vert \varphi'(y) \vert^2 \,{\rm d}y \right]^{-2}.\]

For more details, see e.g. Fantuzzi & Wynn, Phys. Rev. E 93(4), 043308 (2016).

This example shows how QUINOPT can be used to construt a polynomial background field \(\varphi(y)\) to maximize

\[\mathcal{B}(\varphi) = 2\varphi(1) - \int_0^1 \vert \varphi'(y) \vert^2 \,{\rm d}y,\]

and thereby minimize the bound on the time-averaged bulk energy dissipation, when \(G=1000\) and \(\Lambda=2\). For simplicity, we assume that it suffices to check the integral inequality constraint for \(k_1=2\pi/\Lambda=\pi\). The aim of the example is to illustrate how constraints on the optimization variables can be specified in addition to integral inequality constraints.

Download the MATLAB file for this example

1. Define the problem variables

We begin by clearing the workspace and defining the problem parameters.

>> clear
>> yalmip clear
>> quinopt clear
>> lambda = 2;
>> G = 1000;

Then, we define the independent and dependent variables with the commands

>> y = indvar(0,1);
>> [u,v] = depvar(y);

Then, we construct the polynomial background field \(\varphi(y)\). We will take \(\varphi(y)\) to have degree 20. Since QUINOPT represents polynomials in the Legendre basis internally, we define \(\varphi(y)\) in the Legendre basis directly using the command legpoly().

>> phi = legpoly(y,20);

Moreover, the integral inequality constraint depends on the first derivative of \(\varphi(y)\). This is easily found using the function jacobian():

>> D1phi = jacobian(phi,y);

Finally, we need to specify the boundary conditions on the background field, \(\varphi(0)=0\) and \(\varphi'(1)=G\). These can be specified in a vector of constraints like any standard YALMIP constraint, and the boundary values \(\varphi(0)\) and \(\varphi'(1)\) can be accessed using the function legpolyval():

>> CNSTR(1) = legpolyval(phi,0)==0;              % \phi(0)=0
>> CNSTR(2) = legpolyval(D1phi,1)==G;            % \phi'(1)=G

2. Set up the optimization problem

The integral inequality constraint can be set up, as usual, by defining its integrand and the boundary conditions on the dependent variables:

>> k = pi;
>> EXPR = ( u(y,2)^2+v(y,2)^2 )/k^2 + 2*( u(y,1)^2+v(y,1)^2 ) + k^2*( u(y)^2+v(y)^2 ) - 2*D1phi/k*( u(y)*v(y,1) - u(y,1)*v(y) );
>> BC = [u(0); u(1); u(0,1); u(1,2)];        % boundary conditions on u
>> BC = [BC; v(0); v(1); v(0,1); v(1,2)];    % boundary conditions on v

The objective function \(\mathcal{B}(\varphi)\), to be maximised, can be set up using the command legpolyint() to compute the boundary value \(\varphi(1)\), and the command int() to integrate the square of \(\varphi'(y)\) over \([0,1]\):

>> OBJ = 2*legpolyval(phi,1) - int(D1phi^2,y,0,1)/G;

3. Solve and plot the optimal \(\varphi(y)\)

Having defined all variables and constraints, we can maximize the objective function OBJ using the syntax

>> quinopt(EXPR,BC,-OBJ,[],CNSTR)

Note

The first two arguments specify the integral inequality constraint, while the additional constraints are specified in the fifth argument CNSTR. The fourth argument specifies QUINOPT’s options, and here it is left empty to use the default options. Finally, note the minus sign in the objective function, which is needed because QUINOPT minimizes the specified objective function by default.

Once the problem is (successfully) solved, we can compute the upper bound on the dissipation coefficient \(C_\varepsilon\) with

>> UB = G/(value(OBJ))^2;

to find \(C_\varepsilon \leq 7.48\times 10^{-3}\) approximately. Finally, we can plot the optimal background field \(\varphi(y)\) and its first derivative using the command plot(), which is overloaded on polynomials defined using the function legpoly():

>> subplot(2,1,1)
>> plot(0:0.01:1,phi,'-','LineWidth',1.5);
>> subplot(2,1,2)
>> plot(0:0.01:1,D1phi,'-','LineWidth',1.5);

This produces the figure below; note that the boundary conditions \(\varphi(0)=0\) and \(\varphi'(1)=G\,(=1000)\) are indeed satisfied.

../_images/shearflowBF.png