Energy stability of Bénard-Marangoni conduction at infinite Prandtl number¶
Bénard-Marangoni convection describes the flow of a fluid layer driven by surface tension forces of magnitude described by the Marangoni number \(M\). In this example, we consider the idealized case of a unit-depth layer which is infinite in the horizontal directions, and compute the maximum \(M\) for which the conductive state (when the fluid is at rest) is stable with respect to sinusoidal temperature perturbations with frequency magnitude \(k\). In particular, the energy method shows that conduction is stable if
for all functions \(T(z)\) that satisfy the boundary conditions \(T(0)=0\) and \(T'(1)=0\), where
For more details, see e.g. Hagstrom & Doering, Phys. Rev. E 81, 047301 (2010).
The aim of this example is to demonstrate how to solve problems in which unknown boundary values of the dependent variables appear explicitly in the integral inequality constraints. Moreover, we show how QUINOPT can be used to approximate problems with data that is non-polynomial using a truncated Legendre transform.
Download the MATLAB file for this example
1. Set up the variables¶
First, we clear the workspace, as well as QUINOPT’s and YALMIP’s internal variables:
>> clear
>> yalmip clear
>> quinopt clear
Then, we define the integration variable \(z\in[0,1]\), the dependent variable \(T(z)\), and the Marangoni number \(M\), which is the optimization variable:
>> z = indvar(0,1);
>> T = depvar(z);
>> parameters M
2. Construct a polynomial approximation to \(F_k(z)\)¶
QUINOPT can only solve integral inequalities with polynomial data, but the function \(F_k(z)\) is clearly not a polynomial. Yet, QUINOPT can be used if we replace \(F_k(z)\) with a polynomial approximation of degree \(d\). To do so, we first construct the exact function \(F_k(z)\) as a function handle, and we compute its first \(d+1\) Legendre expansion coefficients using the fast Legendre transform command flt()
. Finally, we use these coefficients to build a polynomial approximation of degree \(d\) using the command legpoly()
. Below, we set \(k=\pi\) and \(d=10\) (this gives an accurate approximation at least up to \(k=5\)).
>> k = pi;
>> d = 10;
>> Fk = @(z)k*sinh(k)/(sinh(2*k)-2*k).*(k*z.*cosh(k*z)-sinh(k*z)+(1-k*coth(k))*z.*sinh(k*z));
>> leg_coef = flt(Fk,d+1,[0,1]); % Compute the Legendre expansion coefficients
>> FkPoly = legpoly(z,d,leg_coef); % Build a degree-d legendre polynomial with coefficients specified by leg_coef
We can easily plot both \(F_k(z)\) and its polynomial approximation using the function plot()
, which is overloaded on polynomials built using legpoly()
:
>> clf; % clear current figure
>> plot(0:0.01:1,Fk(0:0.01:1),'Linewidth',1.5); hold on; % plot Fk(z)
>> plot(0:0.01:1,FkPoly,'.','MarkerSize',12); hold off; % plot the polynomial approximation
>> xlabel('$z$','interpreter','latex','fontsize',12);
>> legend('F_k(z)','Polynomial approximation','Location','southwest');
>> axis([0 1 -0.2 0]);
Note
The syntax LCOEF = flt(FUN,NCOEF,DOMAIN)
computes the first NCOEF
coefficients of the Legendre series expansion of the function specified by the function handle FUN
over the interval specified by the input DOMAIN
. These can be used to construct a polynomial approximation to the function specified by FUN
of degree NCOEF-1
(recall that a polynomial of degree \(d\) has \(d+1\) coefficients). Note that DOMAIN
should be a bounded interval in the form [a,b]
, and the function handle FUN
should only take one input argument with values in the range specified by DOMAIN
.
3. Maximize \(M\)¶
Once a polynomial approximation of \(F_k(z)\) has been constructed, the maximum Marangoni number \(M\) satisfying the integral inequality at the top of the page can be computed with QUINOPT. First, we define the integrand of the integral inequality, and the boundary conditions on the dependent variable:
>> EXPR = T(z,1)^2 + k^2*T(z)^2 + M*FkPoly*T(z)*T(1);
>> BC = [T(0); T(1,1)]; % The boundary conditions T(0)=0, T'(1)=0
Then we maximize \(M\) by calling
>> quinopt(EXPR,BC,-M)
>> value(M)
(note the negative sign in the objective function, which is needed because QUINOPT minimizes the specified objective). The optimal solution is found to be \(M\approx 78.55\).
4. Summary¶
In summary, the maximum Marangoni number \(M\) for which a sinusoidal perturbation to the Benard-Marangoni conduction state is stable can be computed with the following lines of code:
>> % Clear the workspace, plus YALMIP's and QUINOPT's internal variables
>> clear
>> yalmip clear
>> quinopt clear
>> % Define the problem variables
>> z = indvar(0,1);
>> T = depvar(z);
>> parameters M
>> % Set k and build a polynomial approximation to Fk(z) of degree d=10
>> k = pi;
>> d = 10;
>> Fk = @(z)k*sinh(k)/(sinh(2*k)-2*k).*(k*z.*cosh(k*z)-sinh(k*z)+(1-k*coth(k))*z.*sinh(k*z));
>> leg_coef = flt(Fk,d+1,[0,1]); % Compute the Legendre expansion coefficients
>> FkPoly = legpoly(z,d,leg_coef); % Build a degree-d legendre polynomial with coefficients specified by leg_coef
>> % Set up and solve the optimization problem
>> EXPR = T(z,1)^2 + k^2*T(z)^2 + M*FkPoly*T(z)*T(1); % the integrand of the inequality
>> BC = [T(0); T(1,1)]; % The boundary conditions T(0)=0, T'(1)=0
>> quinopt(EXPR,BC,-M)
>> value(M)