## Abstract

Energy distributions of high-frequency linear wave fields are often modelled in terms of flow or transport equations with ray dynamics given by a Hamiltonian vector field in phase space. Applications arise in underwater and room acoustics, vibroacoustics, seismology, electromagnetics and quantum mechanics. Related flow problems based on general conservation laws are used, for example, in weather forecasting or in molecular dynamics simulations. Solutions to these flow equations are often large-scale, complex and high-dimensional, leading to formidable challenges for numerical approximation methods. This paper presents an efficient and widely applicable method, called *discrete flow mapping*, for solving such problems on triangulated surfaces. An application in structural dynamics, determining the vibroacoustic response of a cast aluminium car body component, is presented.

## 1. Introduction

Wave energy transport in the high-frequency limit can be described in terms of the underlying ray dynamics, neglecting interference and other wave effects. The problem is thus reduced to tracking ray densities in phase space and becomes part of a wider class of mass, particle or energy transport problems driven by an underlying deterministic velocity field. Applications are found, for example, in fluid dynamics [1], weather forecasting [2], linear wave dynamics or in general in describing the evolution of phase space densities of a dynamical system. The governing equations for the particle or ray densities can often be written in terms of conservation laws [3], an example is the Liouville equation describing the evolution of ray densities in phase space for Hamiltonian flows. This class of problems is particularly interesting in the context of approximating linear wave equations in terms of their underlying ray dynamics in the short wavelength limit.

Numerical approaches for solving transport problems of this kind are typically based on the method of characteristics, that is, the solutions are found along trajectories or rays determined by the underlying vector field. The flow equations can be formulated in terms of a linear propagator, the so-called Frobenius–Perron (FP) operator [4], which describes the evolution of phase space densities in time. The FP operator acting on a phase space density *ρ* may be expressed as an integral operator of the form
1.1Here, the flow map *φ*^{τ} describes the propagation of trajectories over time *τ*. Efficient numerical methods for solving flow problems in more than one-dimension for a wide class of physically relevant systems are still non-existent. A variety of techniques have been developed based on an FP-operator approach, however, all with a fairly limited range of applicability. Difficulties arise due to the high-dimensionality of the phase space and the singular nature of the operator describing the underlying deterministic dynamics. One approach for dealing with such problems is Ulam's method [5], which is based on subdividing the phase space into distinct cells and considering transition rates between these phase space regions. Other methods include wavelet and spectral methods for the infinitesimal FP-operator [6,7] and periodic orbit expansion techniques [4,8]. The modelling of many-particle dynamics, such as protein folding, has been approached using short trajectories of the full, high-dimensional molecular dynamics simulation to construct reduced Markov models [9]. For a discussion of convergence properties of the Ulam method in one- and several dimensions, see Bose & Murray [10] and Blank *et al.* [11], respectively.

More direct methods are based on tracking swarms of trajectories in phase, space often referred to as *ray tracing* [12]. Methods related to ray tracing but tracking the time-dynamics of interfaces in phase space, such as moment methods and level set methods, have been developed by Osher & Fedkiw [13], Engquist & Runborg [14], Ying & Candès [15] and Boon *et al.* [16] among others. They find applications in acoustics, seismology and computer imaging, albeit restricted to problems with few reflections; for an excellent overview, see Runborg [17]. In the following, we focus on ray-tracing approximations of linear wave problems, although the methodology developed here can be used in a more general context.

Ray-tracing and tracking methods often become inefficient when considering stationary (wave) problems in bounded domains, or in general for ray-tracing problems including multiple scattering trajectories and chaotic dynamics. An example is the wave field in a finite cavity driven by a continuous monochromatic excitation. Here, multiple reflections of the rays and complicated folding patterns of the associated level-surfaces often lead to an exponential increase in the number of branches that need to be considered. It is thus necessary to use approximation methods for the associated ray-tracing solutions, often based on ergodicity and mixing assumptions of the underlying ray dynamics. A popular tool among the mechanical engineering community in the context of vibroacoustic modelling is statistical energy analysis (SEA) [18–20]. A related method for electromagnetic fields is the random coupling model, which makes use of random field assumptions [21]. In SEA, the structure is subdivided into a set of subsystems and ergodicity of the underlying ray dynamics as well as quasi-equilibrium conditions are postulated. The result is that the density in each subsystem is taken to be approximately constant leading to greatly simplified equations based only on coupling constants between subsystems. The disadvantage of these methods is that the underlying assumptions are often hard to verify *a priori* or are only justified when an additional averaging over ‘equivalent’ subsystems is considered. The shortcomings of SEA have been addressed by Langley [22], Langley & Bercin [23] and more recently in a series of papers by Le Bot [24–26].

A computational method called *dynamical energy analysis* (DEA), which is based on an FP-operator approach, has been introduced by Tanner [27] and further developed in Chappell *et al*. [28]. DEA is a Ulam-type method subdividing configuration space into smaller subsystems, but using a basis approximation (spectral method) to obtain an improved resolution of the phase space density. By increasing the resolution in both the position and momentum variable, one systematically interpolates between SEA and full ray tracing, thus relaxing the underlying ergodicity and quasi-equilibrium assumptions in SEA. A more computationally efficient approach using a boundary element method for the spatial approximation has been applied to both two- and three-dimensional problems in Chappell and co-workers [29,30]. A major advantage of DEA is that by removing the SEA requirements of diffusive wave fields (equivalent to the ergodicity assumption) and quasi-equilibrium conditions, the choice of subsystem division is no longer critical. The resulting increase in flexibility in this choice leads to a much more widely applicable method.

In this work, we exploit the freedom in choosing the subsystems by extending DEA towards solving the Liouville equation (or other hyperbolic equations) on meshes. Meshed surfaces are replete in numerical simulation problems, thanks largely to the huge popularity of finite-element methods (FEMs). Considering the elements of a mesh as our basic domains therefore automatically renders our techniques applicable to a wide class of problems including the ability to handle complex geometries; an example of which is given in figure 1*a*. A highly efficient solution procedure can be developed for triangulated surfaces since the problem is simplified to considering the local ray dynamics in flat planar regions. In addition, solving phase space flow equations on meshes makes it possible to consider geodesic dynamics on curved surfaces or ray dynamics in non-homogeneous media. This opens up the range of possible applications enormously, for example, in modelling high-frequency vibrations of thin elastic shells [31,32] or in underwater acoustics [33].

The main outcome of this paper will be the introduction and development of the *discrete flow mapping* (DFM) technique, a new and efficient method for solving stationary phase space flow equations in complex domains. DFM has similarities to both finite-volume methods and boundary integral methods, combining the advantages of each. One achieves the reduction in dimensionality to the boundary of each subdomain characteristic of boundary integral methods, which in phase space is a reduction by two. In addition, one can treat non-homogeneous domains as in finite-element and finite-volume methods by applying an approximation of local homogeneity and refining the mesh to improve this approximation. A great advantage of DFM is that it can be applied directly to existing finite-element models with relatively coarse meshes, and provides a solution for the high-frequency case that automatically contains the geometrical details absent from SEA-type methods. In addition, the flow directivity can be resolved, unlike in an SEA model where ergodicity is assumed.

The paper is structured as follows: we will give an integral equation formulation for the stationary Liouville equation on triangulated surfaces in §2. We then detail the implementation of DFM in §3; using the linearity of the integral operator, we approximate the solution using a mixture of boundary element and spectral methods. We also discuss the treatment of reflection and transmission on interfaces, such as at edges or due to abrupt changes in material parameters. Finally, in §4, DFM is applied to model the phase space flow on a sphere and the geodesic flow on an irregularly shaped car body part, demonstrating its power and efficiency in practice.

## 2. Integral equation formulation

We will focus on triangulated surfaces consisting of *N*_{Ω} triangles *Ω*_{j}, *j*=1,…,*N*_{Ω} such as depicted in figure 1*a*. We consider piecewise constant Hamiltonians of the form *H*_{j}(*r*,*p*)=*c*_{j}|*p*|=1 in *Ω*_{j} describing the energy of a flow, where *c*_{j} is the flow velocity for *r*∈*Ω*_{j}, and the momentum coordinate *p* lies on a circle of radius . This Hamiltonian is associated to the Helmholtz equation with inhomogeneous wave velocity *c*(*r*) [17]. In §3*b*, we also consider vectorial wave equations, such as plate equations, which include different wave modes. In this case, the wave propagation needs to be characterized by more than one Hamiltonian per triangle *j*, that is, we need to consider Hamiltonians , where *l* refers to the mode type. We restrict our discussion to the scalar case in the following for simplicity of notation.

Let us denote the phase space on the boundary of the triangle *Ω*_{j} as . The associated coordinates are *X*_{j}=(*s*_{j},*p*_{j})∈*Q*_{j} with *s*_{j} parametrizing ∂*Ω*_{j}, the boundary of the *j*th triangle, and parametrizing the component of the inward momentum (or slowness) vector tangential to ∂*Ω*_{j}. We denote the boundary flow map as which takes a vector in *Q*_{j} and maps it under the flow given through *H*_{j} to a vector in *Q*_{i} (figure 1*b*). Note that *φ*_{ij} is generally only defined on a subset of *Q*_{j}, namely the preimage of *Q*_{i}, that is . This preimage is empty if *Ω*_{j} and *Ω*_{i} are not adjacent (figure 1).

The stationary density *ρ*(*X*_{i}) on *Q*_{i}, *i*=1,…,*N*_{Ω}, due to an initial boundary distribution *ρ*^{(0)} on *Q*_{j}, *j*=1,…,*N*_{Ω}, may be computed using the following boundary integral equation [27,29,30],
2.1where
2.2and *K*_{Γ} describes the propagation of the flow. Here, we consider purely deterministic flows with
Note that the transfer operator is a modified version of the FP operator (1.1) defined over a triangulated surface. The time-dependent flow map *φ*^{τ} has been replaced with the boundary flow map *φ*_{ij}. In addition, a weight function *w* has been included to describe reflection/transmission probabilities at boundaries and to incorporate dissipation. In this work, we consider a dissipative factor of the form , where *L* is the length of the flow trajectory and *μ* is a damping coefficient. A diffusion component may be added by replacing the *δ*-distribution in *K*_{Γ} with a finite-width kernel.

For the case *w*=1, the transfer operator is of FP type (1.1) with a maximum eigenvalue of 1. Here the problem (equation (2.1)) becomes ill-posed, and thus we consider the case *μ*>0 henceforth. One can also obtain a well-posed problem using other forms of dissipation, such as an absorbing or open boundary region [29]. Note that because the flow map only maps to neighbouring triangles, a matrix representation of over the whole of is, in general, sparse. Equation (2.1) is a boundary integral formulation for the stationary Liouville equation as detailed in Chappell & Tanner [30]. The derivation of *ρ*^{(0)} for a high-frequency point source is also given in Chappell & Tanner [30].

## 3. Implementation of discrete flow mapping

### (a) Discretization

Typically, it is assumed that *Ω* consists of a large number of triangles *N*_{Ω} describing a geometrically complex domain. A discrete approximation of *ρ* is sought in phase space coordinates on the triangle boundaries. The spatial approximation is given by taking polynomial basis functions on each triangle edge *α* with *α*=1,2,3. That is, we split the approximation at corners, see Chappell *et al.* [28] for further discussion on the benefits of this choice. In the following, we use piecewise constant functions on the triangle edges for the sake of simplicity. In contrast to the position coordinate, the momentum coordinate has support on the interval only. It is therefore proposed to use a Legendre polynomial basis approximation in this coordinate [29]. A key advantage of these choices is that the integrand in the operator (2.2) remains a very simple function of the position argument and the corresponding integral can be performed analytically. This dramatically reduces the costs of evaluating (2.2) compared with the implementation used in Chappell *et al.* [29].

The overall approximation on *Q*_{i} for *i*=1,…,*N*_{Ω} is then of the form
3.1where *N*_{p} is the order of the momentum basis expansion. The momentum basis functions are given by
3.2where *P*_{β} is the Legendre polynomial of order *β*. The piecewise constant spatial basis functions are given by for *s*_{i} on the edge *α*, and zero elsewhere. Here, *A*_{α} is the length of the edge *α*∈{1,2,3} of *Ω*_{i}. Imposing a weak form of the integral equation (2.1) using the orthonormal inner product for Legendre polynomials 〈⋅,⋅〉 yields
3.3where *I* is the identity matrix,
3.4and the vectors , have entries given by *ρ*_{(i,α,β)}, , respectively. Expanding the inner product in (3.4) and using the definition of the transfer operator (2.2), the discretized transfer operator *B* acting on *Q*_{i} for all *i*=1,…,*N*_{Ω} may be expressed as
3.5Here, we write to denote the splitting of the position and momentum parts of the boundary map. Using the properties of the weight function and the spatial basis functions, equation (3.5) simplifies to
3.6

Here, *w*_{rt} is equal to the transmission probability *w*_{t} when *i*≠*j*, and is equal to the reflection probability *w*_{r}=1−*w*_{t} otherwise. Higher-dimensional transmission/reflection matrices arise in the case of multi-mode dynamics as discussed in the next section. Also, *L*_{ij} is the length of the trajectory from the point represented by *s*_{j} to , and is the admissible range of values for *s*_{j}. Restriction to this range is necessary to ensure that a ray starting on edge *l* of triangle *j* with direction coordinate *p*_{j} will be transferred to a particular adjacent edge *α* of triangle *i* as shown in figure 2. Note that for flat polygonal elements such as triangles, in fact only depends on *p*_{j}, and hence only the damping term in equation (3.6) retains dependence on *s*_{j}. The inner integral thus has a simple form and can be computed analytically, leaving only a single integral to be evaluated numerically. It is this step that leads to vast improvement in the computing times, and enables us to consider many thousands of elements with high-order approximations in direction space. Two key points are:

— Triangulation is used to ensure that the elements of a piecewise linear mesh are flat. This enables us to deal with complex enclosures/surfaces while keeping the system locally simple.

— Owing to the semi-analytic phase space integration, working with many thousands of subsystems and high-order direction space approximations becomes tractable.

It is worth pointing out that the analytic integration method is not restricted to triangular elements, but works for any flat polygons. The general idea of the DFM is related to the work in Chappell *et al.* [28] and Chappell & Tanner [30]. However, in these papers, the DEA method was developed for general subsystems of arbitrary shape and the double integrals in the equations equivalent to equation (3.6) needed to be computed entirely numerically. Computing the matrix *B* is then time consuming, and the computations were restricted to a small number of subsystems (typically 10 or less).

### (b) Transmission at interfaces and ray tracing on curved surfaces

For modelling high-frequency wave propagation, it will be necessary to take into account transmission (and consequently reflection) at interfaces with abruptly changing material parameters or owing to edges and branch lines. The latter will arise, for example, when plates are connected along a common junction such as depicted in figure 3. These effects are included in our approach through transmission/reflection probabilities at triangle boundaries coinciding with the interface. The reflection/transmission coefficients are obtained from local wave solutions at the interface, incorporating the dependence on the direction of the incoming and scattered waves. A number of scenarios important in the context of both wave propagation and ray tracing are outlined below.

*Scalar wave equations*. For scalar wave equations and simple planar interfaces, the transmission coefficients are derived from the wave equation using continuity of the wave function and its normal derivative for an incident plane wave [34]. One thus obtains
3.7Here, *k*_{j} and *k*_{i} are the wave numbers in the elements containing the incoming and outgoing rays, respectively, with *k*_{i}=*ω*/*c*_{i} and *ω* is the angular frequency. A change of wave vector at the interface may be due to a change in the material parameters for example. The angle of incidence *θ*_{j} of the incoming ray with respect to the normal at the element boundary is simply . The direction of the outgoing ray *θ*_{i} is determined using Snell's Law.

*Elastic waves in coupled plates*. A more general case is represented by interfaces forming junctions between plates of varying thickness and meeting at arbitrary angles. The transmission/reflection coefficients for a set of plates coupled at a common interface (such as depicted in figure 3) may be computed using the methods presented in Langley & Heron [35] and Craik *et al.* [36]. In particular, we consider the connection between plates as line junctions, that is, the interior properties of the junction are not modelled and, the mass and moment of inertia are neglected. Let us consider a line junction that couples *n* different (for simplicity) semi-infinite plates. The boundary conditions at the line junction correspond to dynamic conditions involving stresses, moments and kinematic conditions for the displacement and rotation of the *n* plates. To construct the transmission coefficients, we calculate the response of the system with respect to excitation by an incoming plane wave. The incoming wave has a fixed wavenumber and a characteristic mode, that is, it is of bending (*b*), pressure (*p*) or shear (*s*) type. The outgoing waves typically have components in all *n* plates and are a mixture of all mode types. Evanescent modes may be included to complete the description. Possible material differences between the plates can lead to different wavenumbers in different plates. For a given forcing with a particular incoming mode in a particular plate, we can solve for the unknown modal coefficients in all plates. In practice, we find the transmission probabilities directly by calculating the ratio of outgoing to incoming normal power fluxes. A detailed description can be found in Chappell *et al.* [37].

*Curved surfaces*. In the case of a geodesic flow on a smooth curved surface, it is necessary to mimic the geodesic paths on the corresponding triangulated surface. We use the theory outlined in Kimmel & Sethian [38] and Martinez *et al.* [39], making use of the fact that for ray paths not intersecting vertices on triangulated surfaces, the notions of shortest and straightest (discrete) geodesic are equivalent [39]. Hence, the straightest geodesic choice of *θ*_{j}=*θ*_{i} will approximate the direction of the geodesic flow on homogeneous regions of the triangulated surface, and abrupt changes of surface thickness or material will result in a Snell Law effect [38]. A suitable choice of quadrature method here ensures that the rays never pass through vertices of the triangulation (although a point source may lie on a vertex), for example Gaussian quadrature rules where endpoints are never used as abscissae.

*Curved shells-curvature corrections*. For modelling the vibration of thin shells, we need to consider ray tracing on curved surfaces, where the dynamics are derived from thin shell theory [31]. In a lowest order approximation, curved rays again follow the geodesics of the surface [32,40]. This approximation derives from applying thin shell theory in the high-frequency limit as in Norris & Rebinsky [41] and Norris [31], which is valid for wavelengths shorter than the radii of curvature of the shell, but larger than the thickness. Corrections to the geodesic ray approximations need to be considered if the local radii of curvature are of the same order as the wavelength. It is possible to construct modified ray paths from the dispersion curves given by thin shell theory [41], but this requires a detailed knowledge of the local curvature. In the interests of keeping the model as simple as possible, we will follow a different approach here. We treat the meshed structure as a set of plate-like elements and estimate reflection/transmission properties owing to the finite angle between mesh elements using the plate-junction theory sketched earlier and detailed in Langley & Heron [35] and Craik *et al.* [36]. This enables us to mimic wave barriers due to regions of high curvature as demonstrated in §4.

In each of the cases considered earlier, reflection/transmission coefficients are incorporated as part of the weight function *w* in the boundary integral kernel *K*_{Γ} in equation (2.2), or its finite dimensional approximation. We assume that the transmission coefficients depend only on the incoming momentum *p*_{s} via the angle of incidence; the outgoing angle is given by Snell's law taking into account refraction owing to differences in the wave speed across the interface.

## 4. Numerical examples

We will demonstrate the efficiency and flexibility of DFM with the help of two examples. First, we determine the ray density produced by a point source on a sphere, where an analytic solution is available for verification. Second, we compute the energy density distribution on the curved surface of a cast aluminium Range Rover body part and compare with solutions of the corresponding wave equation obtained using the FEM.

### (a) Flow on a sphere

The ray density *ρ* generated by a geodesic flow emanating continuously from a point source on a sphere can be determined analytically. As a function of the polar angle *ϕ*, with source point *ϕ*=0, the ray density is given by:
4.1Here, *μ* is the damping coefficient introduced in §2 and *C* is a constant depending of the strength of the source. The derivation of equation (4.1) follows simply from the fact that the ray density on the sphere is inversely proportional to the element of surface area. The exponential terms result from damping contributions of the form ; summing over the trajectory lengths *L* at any point on the sphere results in a geometric series with a contribution each time the ray orbits the sphere.

The sphere is thus a good candidate for verifying our approach on an approximately spherical triangulated surface. The dynamics on the sphere are integrable and the exact solution for the density contains a singularity at both poles (*ϕ*=0 and *π*) as the rays are unidirectional along great circles passing through the source point. This example is therefore slightly atypical because ergodic or mixing ray dynamics in complex geometries generally lead to more smoothly distributed ray densities. The sphere is thus a true challenge for DFM, which due to the finite-order basis approximation will always incorporate diffusive behaviour and consequently a smoothing of singularities.

Figure 4 demonstrates that DFM can deal even with such a singular case. Higher-order implementations of the basis approximation in direction space lead to an improved resolution of the refocussed singularity. Table 1 gives the mean relative error averaged over the upper hemisphere containing the source point shown in figure 4*a*. The results are given for Delaunay triangulations with differing numbers of triangles and for different orders of direction space basis. The results are computed at the centroid of each triangle and compared against the exact solution for the same value of *ϕ*. For the results presented here, we have taken *μ*=1 and *C*=(800*π*^{2})^{−1}, which is the scaling for a unit excitation of the Helmholtz equation with *k*=100*π*. The factor is derived from the high-frequency asymptotics of the two-dimensional fundamental solution for the Helmholtz equation [30] and matching the asymptotics with equation (4.1) as the distance from the source becomes small.

### (b) An application in vibroacoustics

In this section, we consider the transport of high-frequency flexural or bending wave energy through curved structures using thin shell and high-frequency ray models. At present, there is wide interest in modelling aluminium shells in the automotive industry owing to a desire for lighter and hence lower emission vehicles. In particular, large moulded aluminium components are replacing more traditional multi-component beam-plate constructions, which has the additional advantage of eliminating problems due to fatigue at joints. High-frequency vibroacoustic models based on an SEA treatment will be unsuitable in these circumstances, because complex geometrical features are not included in SEA. In addition, a subdivision of the model into subsystems is not clear cut for such castings as all components are well connected, see, for example, the structure in figure 1*a* representing an aluminium shock tower from a Range Rover.

DFM can overcome these problems, since it can be easily applied in the framework of existing grids for finite-element models, requires no choice of subsystem division and incorporates the full geometry and directionality of the energy flow. The response of a thin moulded aluminium car component (range rover shock tower) to a point force applied normally to the surface is modelled and compared with an FEM approximation for the full wave model performed using Nastran. Figure 1*a* shows the problem set-up, including the mesh and indicating the local DFM map *φ*_{ij}. In order to maintain a tractable model size for the finite-element calculation and to study frequency ranges of industrial interest, the computation is performed at frequencies between 8 and 10 kHz. This approximately corresponds to a one-third-octave band centred at 9 kHz. A key assumption of the thin shell theory is that the radius of curvature is large compared with the wavelength [41], which is only partly satisfied for this structure in the frequency band considered. The wavelength is typically around 3–6 cm depending on the shell thickness. We therefore use a modified geodesic-ray-tracing technique by incorporating reflection/transmission owing to finite angles between mesh elements as described in §3*b*. This leads to geodesic trajectories without reflections in relatively flat areas, but incorporates wave barrier behaviour in regions of high curvature.

Figure 5 shows the DFM result (bending mode only) compared with the Nastran solution for the kinetic energy averaged over 41 evenly spaced frequencies spanning the prescribed range and with a hysteretic damping level of 3 per cent, which is typical for such a structure. The point source is positioned at the front of the structure. The ray computation is performed using a triangulated surface consisting of 11 623 triangles and with a sixth-order Legendre polynomial basis in direction space. The Nastran grid contains 40 670 elements comprising a mix of piecewise linear triangles and quadrilaterals. As would be expected, one sees more oscillation in the full wave model. The DFM prediction of the overall energy flow in regions of both high and low curvature matches well with the Nastran results. The flow of energy along the side walls of the central mount, as well as features such as shadow regions owing to holes in the structure and channelling effects owing to variations in the shell thickness can be observed as common in both plots. Such geometrically dependent features would be entirely absent from SEA-type models and represent a major advance for high-frequency vibroacoustic simulation methods.

Figure 6 gives the kinetic energy averaged over four different receiver regions, which are displayed as the blue regions on the inset structure plots. The response is now shown for a range of damping levels between 1 and 6 per cent, and for receiver regions with differing levels of proximity to the source point. The correspondence between the mean of the 41 Nastran solutions and DFM is remarkably good in figure 6*a–c*. The regions are spread across the whole structure and provide strong verification of our approach. In figure 6*d*, the DFM results deviate slightly for high damping. Note that this region is close to the region in figure 6*c* and the deviations thus reflect local variations in the FEM solution. The results clearly demonstrate that curvature related wave-barrier effects are correctly predicted by DFM using local reflection/transmission matrices.

The computations in this section were performed in parallel on four cores of a desktop PC in less than 6 min. The code was written in C++ with openMP. The large sparse non-symmetric linear system of dimension 732 249 that arises has been solved efficiently and accurately using a stabilized biconjugate gradient iterative solver.

## 5. Conclusions

We present DFM as a highly efficient method for solving stationary phase space flow equations in complex domains and apply the method to two illustrative examples. In particular, we highlight the versatility of the method, its capability to handle large-scale models efficiently and emphasize the fact that both geometrical details and flow directivity are fully included in the calculation at a moderate computational cost. In the special case of ray focusing, the method clearly captures the foci, albeit not reproducing the actual singularity. In the vibroacoustic example, the full complexity of the model is taken into account via the mesh functionality of DFM. Deviations from a pure geodesic flow owing to regions of high curvature have been reproduced without compromising the efficiency of the method. DFM can thus serve as a practical tool for solving phase space transport problems with applications in engineering ranging from acoustics to structural mechanics. Extensions of the method to electrical engineering, as well as to fluid mechanics and other flow problems are evident.

## Acknowledgements

Support from the EPSRC (grant no. EP/F069391/1 and a follow-on knowledge transfer grant) and the EU (FP7 IAPP grant MIDEA) is gratefully acknowledged. The authors also gratefully acknowledge the support, guidance, data and Nastran access provided by Stephen Fisher and Jaguar Land Rover.

- Received March 6, 2013.
- Accepted April 16, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.