Casting Computational Fluid Mechanics into a Convex Quadratic Optimization Framework

Although variational methods have demonstrated success in many fields of mechanics (Lanczos, 2012), their implementation in fluid mechanics remains relatively limited (Penfield Jr, 1966). This stems from the difficulty of incorporating dissipative force such as viscosity using the standard Hamiltonian formulation (Bretherton, 1970; Salmon, 1988; Morrison, 1998; Morrison et al., 2020). As a result, traditional computational methods in fluid mechanics remain predominantly confined to Newtonian mechanics where the primary focus is on the solution of the Navier-Stokes equation (NSE) or one of its derivatives (Gresho & Sani, 1987).

In a departure from traditional variational principles that are based on least action, Taha et al. (Gonzalez & Taha, 2022; Taha et al., 2023) used Gauss’ principle of least constraint (Gauß, 1829; Papastavridis, 2014; Udwadia & Kalaba, 1992, 2002, 1996; Udwadia, 2023) to develop a new variational principle in fluid mechanics known as the principle of minimum pressure gradient (PMPG). This principle asserts that an incompressible flow evolves from one instant in time to another such that the total magnitude of the pressure gradient over the domain is minimized. Taha et al. (Taha et al., 2023) proved that the necessary condition to minimize the cost of the pressure gradient is guaranteed to satisfy the NSE, which transforms a fluid mechanics problem into an optimization problem. Taha et al. (Taha et al., 2023) employed PMPG to find analytical solutions to some classical problems in fluid mechanics, including unsteady laminar flow in a channel and Stokes’ second problem. The attained solutions provided additional insights that are otherwise difficult to infer using traditional approaches (Gonzalez & Taha, 2022; Shorbagy & Taha, 2024).

The development of computational frameworks based on PMPG remains limited and is still in its early stages. Recent work has utilized Physics-Informed Neural Networks (PINNs) to solve the optimization problem resulting from PMPG (Sababha et al., 2024; Alhussein & Daqaq, 2024; Atallah et al., 2024). The proposed method, validated against benchmark examples, demonstrated two key advantages. First, it eliminates the pressure-velocity coupling, thereby alleviating the computationally expensive step of solving the Poisson equation at every time step, which often dominates the runtime in traditional algorithms. Second, it removes the need to impose nonphysical outflow boundary conditions, allowing for a reduction in the computational domain required for most fluid mechanics problems. Despite their advantages, the use of PINN for optimization presents several challenges. First, accurately capturing the complexity of fluid flows often necessitates large neural networks, significantly increasing the number of trainable parameters, such as weights and biases. Second, constraints are typically enforced in a soft manner, making it challenging to guarantee their satisfaction throughout the entire domain. Finally, the optimization process remains inherently non-convex, which can result in multiple, potentially inconsistent solutions, as highlighted in previous studies (He et al., 2023; Wang et al., 2023).

In this work, we transform the PMPG computational framework into a quadratic optimization framework. By focusing on the local acceleration as the primary variable, we reformulate the minimization problem into a convex structure subject to linear constraints. This permits the implementation of well-established quadratic programming tools for optimization or a direct solution through the well-known Karush-Kuhn-Tucker (KKT) condition. We demonstrate the effectiveness of this approach and highlight its simplicity using three benchmark examples; namely the flow in a led-driven cavity, the Poiseuille flow, and the flow past a backward step.

The remainder of this paper is organized as follows. Section 2 provides a review of the theoretical foundation of the PMPG. Section 3 presents the formulation of the problem. Section 4 presents the results. Section 5 presents concluding remarks that discuss the advantages and limitations of the current approach, as well as potential extensions.

In the context of the constrained motion of particles, Gauss’ principle of least constraint states that the constrained system accelerates so that it minimizes, in a weighted least-squares sense, the deviation between its current acceleration and that of the free motion (Gauß, 1829). Thus, a particle adjusts its path only to the extent necessary to meet the constraints, ensuring the least possible deviation from the unconstrained trajectory; i.e., its free/natural motion.

In mathematical terms, for $N$-constrained particles, each of mass, $m_{i}$, whose motion is described by q generalized coordinates, Gauss’ principle asserts that the quantity

also known as the Gaussian cost, is a minimum with respect to the generalized accelerations, $\ddot{\mathbf{q}}$, at every instant of time. Here, $\mathbf{a}_{i}$ and $\mathbf{F}_{i}$ are, respectively, the acceleration of the $i^{th}$ particle, and the net (non-constraint) force acting on it. Gauss’ principle, therefore, turns the mechanics problem governed by Newton’s Equation into the instantaneous minimization problem.

Using the NSE and the definition of the Gaussian cost, Taha et al. (Taha et al., 2023) extended Gauss’ principle of least constraint to a continuum of fluid particles forming an incompressible fluid. They showed that, in the case of some fluid domain, $\Omega$, the Gaussian cost can be written in Eulerian coordinates as:

where, $\rho$, is the density of the fluid, $\mathbf{u}$ is the fluid velocity vector, and $\nu$ is the kinematic viscosity.

Taha et al. (Taha et al., 2023) further demonstrated that the solution, u, which satisfies NSE and the continuity, $\nabla.\mathbf{u}=0$, subject to any boundary conditions, is the same as the one that minimizes the functional, $\mathcal{A}$, in Equation (2) with respect to the local acceleration, $\frac{\partial\mathbf{u}}{\partial t}\equiv\mathbf{u}_{t}$ subject to $\nabla.\mathbf{u}_{t}=0$ and the same boundary conditions. Thus, for an incompressible fluid domain, $\Omega$, bounded by the surface, $\Gamma$, the problem of finding $\mathbf{u}$, which satisfies continuity and NSE can be cast in the following form:

Find $\mathbf{u}_{t}$ such that

A fundamental advantage of this formulation is that the expression in Equation (2) is quadratic in $\mathbf{u}_{t}$, and the constraints are linear. Thus, the solution scheme lends itself naturally to the quadratic optimization family whose solution can be obtained using standard optimization techniques that are straightforward to implement and computationally efficient.

To solve the minimization problem, we consider a structured domain with equispaced collocation points arranged in an $n_{x}\times n_{y}$ grid. At each node, the optimization variables are $\mathbf{u_{t}}=(u_{t},v_{t})$, representing the local accelerations in the $x$ and $y$ directions, respectively. Thus, the total number of optimization variables is $2n_{x}n_{y}$. The optimization vector, here denoted as $\mathbf{w}$, is constructed by first flattening the $u_{t}$ and $v_{t}$ components from the grid, and then concatenating them; transforming the representation from a structured grid into a single optimization vector as:

Using the mean rule and the definition of the vector $\mathbf{w}$, the integral of the objective function $\mathcal{A}$ can be re-written in a quadratic form as

where $H$ is an $2n_{x}n_{y}\times 2n_{x}n_{y}$ positive definite matrix often referred to as the Hessian matrix, and $f$ is a $2n_{x}n_{y}\times 1$ vector representing external forces arising from the convective and viscous terms. It is worth noting that the objective function in this form excludes constant terms, as they do not influence the optimization process.

Next, we rewrite the continuity constraint (Equation (2.4)) in terms of the vector $\mathbf{w}$ using the finite difference method. In particular, we approximate the continuity equation at an interior point $(x_{i},y_{i})$ by using the forward difference formula:

Using Equation (8), the continuity equation can be expressed as $D\textbf{w}=0$, where the elements of the matrix $D$ at each interior point $i,j$ have non-zero entries at $D_{i,i-1}=-\frac{1}{h}$, $D_{i,j-1}=\frac{1}{h}$ and $D_{i,j}=\frac{1}{h}$. The rest of the elements are zero making the system extremely sparse.

Finally, the boundary conditions are represented by the condition $B\textbf{w}=0$, where the matrix $B$ is chosen to ensure that $\mathbf{w}$ vanishes along the four boundaries. By combining the continuity and boundary constraint matrices, the augmented system of constraints becomes:

Equations (7) and (9) represent the standard form of a quadratic optimization problem, which can be solved in a finite amount of computational time (Nocedal & Wright, 2006b). The computational effort required to find a solution is influenced by the characteristics of the objective function. In this case, the matrix $H$ is positive definite, making the problem strictly convex. Strictly convex problems are computationally similar to linear programming problems, which are well-known for their efficient solution methods (Nocedal & Wright, 2006a).The solution methodology is outlined in the flowchart 1.

To demonstrate the effectiveness of the proposed approach, we consider three benchmark examples in fluid mechanics: $i)$ the unsteady flow inside a lid-driven cavity (Koseff & Street, 1984), $ii)$ the Poiseuille flow (the unsteady laminar flow in a channel of uniform cross-section), and $iii)$ the unsteady flow in a backward facing step (Biswas et al., 2004). As illustrated in Fig. 2, the characteristic diameter, $D$, is normalized to 1, the inlet velocity, $U_{0}$, is set to 1, and the kinematic viscosity, $\nu$, is chosen such that the Reynolds number is 100 for all examples.

The optimization can be carried utilizing well-established quadratic programming. Here we use MATLAB’s quadratic programming toolbox, employing the interior-point-convex algorithm. This algorithm is tailored for large-scale optimization problems, utilizing an advanced sparse linear solver to efficiently address the substantial computational demands. The reliance on sparse solvers is particularly beneficial for problems of this nature, where the number of variables is high, but the matrices $H$ and $A$ exhibit significant sparsity (Nocedal & Wright, 2006a).

Alternatively, the vector $\mathbf{w}^{*}$ which minimizes Equation (7) subject to Equation (9) is given by:

where $H$ is positive definite¹¹1For an equispaced grid similar to the one considered in this analysis, $H$ turns out to be the identity matrix and is therefore guaranteed to be positive definite., and $\boldsymbol{\lambda}^{*}$ is the vector of Lagrange multipliers. Equation (10) is a consequence of the general result of the first-order optimality condition commonly referred to as the Karush-Kuhn-Tucker (KKT) condition (Nocedal & Wright, 2006b).

Simulations for the lid-driven cavity were performed over the time range $0\leq t\leq 1\,\text{s}$. The minimized cost objective $\mathcal{A}$ was compared with the values obtained from a traditional high-fidelity CFD model (see ref. (Sababha et al., 2024)), revealing excellent agreement between the two approaches. As shown in Fig. 3(a), $\mathcal{A}$ initially decreases sharply, then increases at a diminishing rate, eventually stabilizing to indicate the steady-state condition, where the local acceleration approaches zero. To further validate the proposed method, we compared the predicted flow field $\mathbf{u}$ at $t=1\,\text{s}$ focusing on the middle cross-section, namely (a-a) as highlighted in Fig. 2. A comparison between the values of $\mathbf{u}$ as obtained using the proposed optimization framework (solid line) and those obtained using the high-fidelity model (dashed lines) demonstrates excellent agreement as shown in Fig. 3(b).

Results for the 2D Poiseuille flow are presented in Fig.4. Simulations were conducted over the time range $0\leq t\leq 1\,\text{s}$, showcasing the evolution of the developing flow at various time intervals. The transition from an initially uniform profile to a fully developed state is clearly observed, with boundary layer growth along the walls. Additionally, the predicted flow field at the outlet is shown, demonstrating a parabolic profile that aligns closely with the analytical solution derived from the Navier-Stokes equations (White & Xue, 2003).

Finally, the results for the backward-facing step are presented, with simulations conducted over the time interval $0\leq t\leq 5\,\text{s}$. The temporal evolution of the streamlines is illustrated in Fig. 5, showcasing the transition from an initially uniform velocity profile to a fully developed flow state. The extent of the primary recirculation region is clearly evident, and its measured length is in excellent agreement with the numerical results reported in reference (Biswas et al., 2004).

This study presents a new approach to solving incompressible fluid mechanics problems. In particular, using PMPG, we transform the CFD problem into a convex optimization framework subject to linear constraints and solve it using well-established quadratic programming tools. A direct solution can also be obtained using the well-known Karush-Kuhn-Tucker condition. The approach is implemented using the well-known unsteady lid-driven cavity problem requiring only a few lines of MATLAB code (see Appendix B), in contrast to the lengthy and complex codes often required by conventional finite volume methods, such as those used in Ansys or OpenFOAM. The results are compared to high-fidelity CFD simulations, demonstrating excellent agreement and showcasing the potential of this framework as an efficient and robust alternative to traditional CFD tools. The simplicity of the proposed method not only makes modeling fluid mechanics more accessible but also expands the potential user base, eliminating the need for custom computational frameworks.

In addition, using PMPG to cast CFD problems into a convex optimization framework offers two key advantages. First, by eliminating the pressure from the formulation, the need for pressure-velocity coupling is removed, thus avoiding the computationally expensive step of solving the Poisson equation after every time step. Second, the convexity of the problem guarantees a solution in finite time.

Compared to the previous work of the authors (Sababha et al., 2024), which employed an optimization framework based on physics-informed neural networks (PINNs), this approach offers distinct advantages. First, it preserves the convexity of the problem, avoiding the challenges of non-convex training procedures in PINNs that often lead to sub-optimal performance due to the presence of local minima. Second, it ensures accurate enforcement of constraints, which are critical for obtaining physically valid solutions. PINN-based methods frequently exhibit constraint violations in certain regions, resulting in approximate but not exact solutions.

Despite these advantages, this formulation has certain limitations, as it is currently limited to square domains. Extending the method to handle irregular domains and/or domains with unstructured collocation points is a critical next step in broadening its utility. Moreover, although the advantages of the method are evident, a rigorous comparison with traditional solvers, such as finite-volume methods, is essential to fully assess the performance and applicability of the method. Benchmarking against these established tools will provide valuable information on the accuracy, computational cost, and scalability of the solution, enabling a clearer understanding of the advantages and limitations of this approach. Additionally, exploring the performance of different optimization algorithms within our framework, such as active set methods or gradient projection methods, could uncover further potential that is yet to be fully explored.

Finally, we believe that this work represents a notable step forward in addressing CFD problems. It illustrates that minimizing the pressure gradient transforms them into a convex optimization framework with guaranteed global minima that can be tackled efficiently using a wide array of well-established robust optimization tools. Furthermore, it appears that this new formulation provides a more natural approach to solving incompressible fluid mechanics. For instance, removing pressure from the formulation inherently eliminates several challenges commonly encountered in computational fluid mechanics. These include the need to solve the pressure Poisson equation (Chorin, 1997; Toutant, 2018), the complexity of prescribing appropriate boundary conditions (Papanastasiou et al., 1992), and issues such as spurious oscillations that arise from interpolation in collocated grids (Zhang et al., 2014). Nevertheless, while pressure is removed from the formulation, it can still be recovered from the Lagrange multipliers in Equation (10), since the pressure serves as a Lagrange multiplier enforcing the continuity constraint (Gresho & Sani, 1987).