Finite differences on staggered grids preserving the port-Hamiltonian structure with application to an acoustic duct

Highlights

• A structure preserving spatial discretization method for open port-Hamiltonian governed by the wave equation is presented.

• The method is extended in 2D for rectilinear and regular triangular meshes.

• Time integration is performed using implicit midpoint rule.

• Numerical results are presented in open and closed loop in the case of boundary control.

• The advantage of the regular triangular mesh over the rectilinear one regarding isotropy is discussed.

Abstract

A finite-difference spatial discretization scheme that preserves the port-Hamiltonian structure of infinite dimensional systems governed by the wave equation is proposed. The scheme is based on the use of staggered grids for the discretization of different variables of the system. The discretization is given in 2D for rectilinear and regular triangular meshes. The proposed method is completed with the midpoint rule for time integration and numerical results are provided, including considerations for interconnection and closed loop behaviors and isotropy comparison between the proposed meshes.

Introduction

Discretization of distributed parameter systems is a key issue for simulation and control purposes. Among all existing and standard methods, the ones aiming at preserving structural invariants and structural properties of the original system are of particular interest when control design is considered. The aim of this paper is to revisit the standard finite-difference scheme in the light of port-Hamiltonian formulations to deal with the structural reduction of 1D and 2D systems driven by the wave equation. Port-Hamiltonian systems (PHS) stem from the representation of internal energy exchanges of multi-physical systems in interaction with their environment. They are particularly well suited to describe complex open systems interconnected through power conjugated port variables. The power preserving interconnection of subsystems is thereby formalized by the notion of Dirac structures [1], [2], which represent a generalization of Kirchhoff's laws or Newton's laws and the velocity continuity in mechanics. For open systems, port-Hamiltonian systems cope with input and output flows and efforts [3]. By using this formalism, controlled systems can be interpreted in terms of the interconnection of physical systems exchanging energy among them, which from an engineering perspective provides a better understanding of the closed-loop system. This in turn permits to derive control strategies that are physically inspired [4]. Infinite-dimensional PHS governed by partial differential equations interconnected through their boundaries were introduced for modelling in [5], [6] and extensions for control in for instance [7], [2], [8], [9], [10]. In the infinite-dimensional case, the Dirac structure underlying the PHS is expressed through differential operators. Performing a discretization which preserves the geometric structure of a PHS allows to preserve conservation laws and thus to preserve properties such as energy conservation and passivity with respect to the natural inputs and outputs of the system, which are central for control purposes. It also respects the physical meaning of the boundary port variables that can be naturally used for the interconnection of open multi-physical systems.

Several recent works tackle this challenge [11], [12], [13], [14], [15]. In [12], [13], the structure of the system is preserved through the use of different (mixed) finite elements for the approximation of different variables. This approach found applications in modeling, reduction and control [16], [17], [18], [19] and was extended through high order polynomial approximations to pseudo spectral approximations in [14]. In [20] a structure preserving finite-volume discretization for PHS was performed for the 1D case. The mixed Galerkin structure-preserving discretization for systems of conservation laws was presented recently in [21].

Finite differences is another important numerical method whose main advantage is its simplicity. It is based on the discretization of differential operators through Taylor series which lead to various schemes with different advantages such as convenience for particular geometries or convergence orders (see [22], [23] for recent reviews). Among them, schemes presenting staggered grids [24], [25], [26] consider separately the state variables accordingly to their respective geometric nature. In [24], a generalized leapfrog structure is used in the time domain on the linear transport equation. It is shown that this method, if stable, preserves the conservation laws, which is an important property for the study of non-dissipative systems. Moreover, conditions for numerical stability are given therein. The proposed two steps methods apply only on closed systems and [20] can be considered as the extension to open systems from the finite volume point of view. For frequency domain modelling in 2D, the authors of [27], [28], [29] consider two staggered square grids, one of the grids being rotated. Changing the orientation of the grid allows to get high order approximations of the derivatives at the price of higher complexity, up to 25-point stencils. A simpler approach is proposed in [30] where finite differences are performed on simple staggered grids. The convenience of such a method for the wave propagation is illustrated on a seismologic example, where it allows to implement a velocity-stress formulation. Indeed, the use of staggered grids allows to impose the effort variables as boundary conditions, e.g. speed and pressure in acoustics, which is convenient for the study of open systems and interconnections, while a traditional finite-difference method would not. For the 2D and 3D cases, the meshing choice, for example the choice of rectilinear or triangular meshing, impacts the performance of the discretization. As mentioned in [31], if the use of rectangular grids leads straightforwardly to a good approximation of a considered PDE, a structured triangular mesh, also called hexagonal mesh [31], [32], can at the cost of a higher complexity provide a better computational efficiency and is a better choice to emulate the isotropy of the spatial propagation of a wave. Indeed, a convenient choice of the discretization scheme through a compact structured mesh leads to a more homogeneous distribution of the propagation in space. In [33], schemes for a regular triangular mesh are proposed, including an explicit second-order method which leads to a simple and intuitive discretization of the first order differential operator, convenient for the discretization of the conservation laws underlying the wave equation. Another scheme, first proposed in [34] and used in [31], relies on a balanced propagation of the wave along the three axes of a regular triangular mesh for the second order differential operator of the wave equation.

Regarding time integration, since a staggered-grid method for PHS results in the approximation of first order time derivatives, explicit Euler, which is generally used in finite-difference time domain (FDTD) methods, is unconditionally unstable. Implicit midpoint, a symplectic method [35], is shown to be suited for non-dissipative PHS in [36]. Moreover, recent work [37] shows that there is a coincidence between the midpoint rule and discrete gradient in the linear case, which respects the power balance of the system.

In this paper, we propose to combine the port-Hamiltonian framework with finite differences on staggered grids to derive control oriented reduced order systems for the 2D wave equation. The differential operators which define the infinite-dimensional Stokes Dirac structure are approximated with a consistency of order 2 [38] by matrices which define a Dirac structure. In that sense, a finite-dimensional system derived from the proposed method is a PHS endowed in a Dirac structure which approximates the original Stokes Dirac structure. This work is motivated by the study of a generic acoustic control problem proposed in [39], [40] and formulated as a port-Hamiltonian system in [41]. The considered system is an acoustic tube of rectangular section which is approximated using a 2D rectangular domain. Distributed actuation at the boundary (see Fig. 1) aims to dissipate part of the energy carried by the acoustic waves propagating inside the tube. Industrial applications can be found for example in aeronautics [42], air conditioning systems, motorised vehicles, etc., see [43] for a more exhaustive list. Even if the development of this work is motivated by this particular physical system, the discretization methods proposed in this paper apply to any linear 2D hyperbolic system.

The paper is organized as follows. In Section 2 the use of staggered grids for the spatial discretization through finite differences is motivated on the 1D wave equation. In Section 3 the 2D case is treated for a rectilinear and for a regular triangular mesh. It is shown that the method preserves the port-Hamiltonian structure. In Section 4 numerical results are presented for the 1D and the 2D case for the rectilinear and the regular triangular meshes. The study is performed under open and closed loop conditions and the isotropy of the two meshes is compared. Finally in Section 5 we present general conclusions and comments of lines of future research. The notation used in this paper is recapitulated in the Appendix. A short recall on the considered time integration methods is given in a second appendix.