Patent Application
20250131159
This application claims the benefit of European Patent Application No. EP 23204348.9, filed on Oct. 18, 2023, which is hereby incorporated by reference in its entirety.
The present embodiments relate to modeling acoustics in an unbounded domain to determine a field variable within a control volume of the domain.
The analysis of the sound characteristics of industrial machines is crucial in various applications (e.g., in the automotive, aerospace, and marine industries). On one side, in many instances, engineers are to design their products paying attention to the sound quality and sound levels, to make the products more attractive for potential customers. On the other side, the increasing awareness on the potential health risks linked to exposures to high sound levels has led to more stringent regulations on sound emissions in urban environment.
In addition, in the industrial world, there is a growing interest towards digital twins. In this context, high fidelity 3D simulation models are used not only during the design phase of a product cycle, but also during the operation phase, where the high fidelity 3D simulation models may be combined with measurement data, through virtual sensing algorithms, to gain further insights into the state of the product. High fidelity 3D simulation models are typically more complex and expensive to develop, but high fidelity 3D simulation models provide a more accurate and realistic simulation of the system being modeled than surrogate models.
Smart virtual sensing algorithms are known from [24].
In the context of this disclosure, a digital twin is a virtual replica of a physical object, system, or process that is used to simulate and analyze its behavior and performance. A digital twin is typically created using data from sensors, cameras, and other sources that capture information about the physical object or system, and combine the captured information with the simulation model of the said object. The digital twin may be used to test and optimize its performance under different conditions. The digital twin enables better monitoring and control of physical systems, allowing for predictive maintenance and improved performance.
The digital twin technology is instrumental in many scenarios (e.g., in structural health monitoring applications). Due to the nature of this problem, compact representations of high-fidelity models that may be used to run stable time domain simulations are to be provided.
Digital twin technology is known from, for example, the publication Hartmann, D., & Van der Auweraer, H. (2022) “The Executable Digital Twin: merging the digital and the physics worlds” arXiv preprint arXiv:2210.17402 or from De Miguel, A. G., Atak, O., & Blanco, M. A., “A Digital Twin for full-field monitoring in multi-channel control with applications to direct field acoustic testing.”
For many (vibro-) acoustic applications, sound waves are to be modeled travelling to infinity, with the sound waves travelling away and not being reflected back.
Some background knowledge known to a skilled person in the technical field of the invention has been provided by the following publications:
A few options are disclosed in the literature to numerically model transient (vibro-) acoustics wave propagation phenomena in unbounded domains; the most popular techniques include modelling the exterior domain using infinite elements and perfectly matched layers (PML).
PML is a technique used in computational acoustics to model the behavior of sound waves in unbounded regions. The PML was originally developed for use in electromagnetic simulations. The PML technique involves constructing a layer around the computational domain that gradually absorbs the waves as the waves propagate away from the domain. This layer is designed to mimic the properties of an infinite medium, effectively eliminating reflections, and preventing wave leakage into the surrounding space.
At continuous level, no reflections take place at the interface between the FEM and the exterior domain, for any angle of incidence [3]. Various PML formulations for time domain simulations exist [10; 7; 5].
However, there are a number of problems associated with such transient PML formulations: a) to make the system matrices frequency independent, new auxiliary variables are introduced in the PML layer. This provides that for transient calculations the number of degrees of freedom (DOFs) is considerably larger than corresponding frequency-domain simulations. This leads to large time and memory requirements to build the ROMs. As another problem, b) these formulations are usually quite complex and difficult to implement in a simulation software, and c) the stability analysis of many of the formulations presented in literature is not easy to carry out. Each formulation requires a separate inspection [5]. For the formulation presented in [10], some terms are to be suppressed to stabilize transient simulations, introducing new approximation errors. As another problem, d) from [26], it is known in the context of model order reduction (MOR) to guarantee that stability properties are preserved after the projection step (e.g., necessary to use a split-basis approach), which results in larger reduced order models. In a numerical experiment, analyzing a plate radiating in the free field, the author compared the efficiency of the PML and infinite elements for time-domain MOR: the author showed that infinite elements resulted in a smaller full order model (e.g., making the reduction procedure faster) and that the number of DOFs in the final ROM was three times smaller than the one obtained from the PML configuration. The added DOFs lead to larger calculation time for the PML formulation.
A number of model order reduction (MOR) techniques are known in the field of engineering and applied mathematics. Some of the most common techniques include 1. Proper Orthogonal Decomposition (POD). This technique reduces the dimensionality of a system by identifying the most important modes of variation in the system using a statistical method. This technique is often used in fluid dynamics, structural mechanics, and other fields. Another technique includes 2. Balanced Truncation. This technique reduces the dimensionality of a system by balancing the energy of the system between its inputs and outputs. This technique is often used in control systems and linear systems. Another technique includes 3. Krylov Subspace Methods. This technique reduces the dimensionality of a system by approximating the response of the system using a subspace of the Krylov sequence. This technique is often used in linear systems and control systems. Another technique includes 4. Model Reduction via Interpolation (MRI). This technique reduces the dimensionality of a system by approximating the transfer function of the system using a set of interpolation points. This technique is often used in control systems, circuit simulation, and other fields. Another technique includes 5. Reduced Basis Methods. This technique reduces the dimensionality of a system by constructing a basis for the solution space of the system using a set of snapshot solutions. This technique is often used in finite element analysis, structural mechanics, and other fields.
These techniques are used to reduce the complexity of mathematical models and simulations, making the mathematical models and simulations easier to solve and more computationally efficient. These techniques are widely used in engineering, physics, and other fields to analyze and optimize complex systems and processes.
The Infinite Element Method (IFEM) is a popular alternative to the PML, particularly competitive for time domain models. The Astley-Leis infinite element formulation [20] leads to frequency independent system matrices and may be easily used for transient simulations. Two types of Astley-Leis formulations exist: 1. the standard Astley-Leis formulation [39] is based on the truncated multipole expansion written in spherical coordinates (e.g., the Atkinson Wilcox expansion [40]); 2. the spheroidal formulation [28] on is based on the truncated multipole expansion written in spheroidal coordinates.
Hereby, the main difference between spherical and spheroidal coordinates is that spherical coordinates describe a point in terms of its distance from a fixed point and two angles, while spheroidal coordinates describe a point in terms of its distance from a fixed point and two angles that define its orientation relative to a more complex shape, such as an ellipsoid.
Both formulations rely on a singular mapping between the infinite physical elements and the corresponding finite parent elements. For each infinite element, the distance from, for example, the virtual source is mapped to a parent element radial coordinate; to form the shape functions, the usual polynomials written in parent coordinates are used. For the standard Astley-Leis formulation to accurately capture the correct solution, the virtual sources should coincide with the center of the smallest sphere circumscribing all the sound sources (e.g., center of radiation) in the problem analyzed. When this constraint is not met, the accuracy of the solution is not guaranteed.
In this context, a parent element is a reference element that is used as a template to construct the finite element model of a physical domain. The parent element may have a simple geometric shape, such as a triangle, quadrilateral, tetrahedron, or hexahedron.
The parent element is defined by a set of nodal points, or nodes, that are used to interpolate the solution within the element. The nodal points may be located at the vertices, edges, or faces of the parent element, and the solution is approximated by a polynomial or piecewise polynomial function that satisfies continuity and differentiability conditions. The parent element may be transformed and mapped onto the physical domain using a set of geometric transformations, such as affine transformations or isoparametric mappings.
In this context, shape functions are mathematical functions used to interpolate the values of a field variable within each finite element.
The use of the standard Astley-Leis formulation permits to use FEM domains of any convex shape for frequency-domain calculations. For transient simulations relying on such formulation, however, the infinite elements should be extruded orthogonally to the exterior FEM boundary envelope to provide the stability of the method [4]. This, together with the above-mentioned criterion on the positioning of the virtual sources, imposes to use a spherical truncation of the FEM domain. For the analysis of elongated scatterers in time domain, this provides that an unnecessarily large mesh is to be used, leading to prohibitive number of DOFs, as shown in
The spheroidal Astley-Leis formulation circumvents this issue by allowing to fit a smaller FEM domain around the scatter.
In [32, 33], it was demonstrated that for such formulation, the radial order of the shape functions often needs to be higher than for the standard Astley-Leis formulation, offsetting the benefits of using a smaller FEM domain.
It is therefore shown that the issue of automatically extruding the infinite elements has not yet been solved; it appears that the issue has not been addressed in prior literature.
A particularly challenging case is when only the FEM mesh is available and the geometry is not modelled with a CAD tool. In this case, the geometric data needed to extrude the infinite elements is to be computed directly from the mesh of the FEM region. A version of the PML (e.g., automatically matched layer (AML)) has been reported in [47] for which the user must only provide the scatterer geometry (e.g., either a CAD model or a FEM mesh); then, the acoustic interior convex domain is generated automatically using the, for example, Quickhull Algorithm, which is known from [36], and the AML is extruded automatically during run time. The AML, however, is only available for frequency-domain simulations.
In this context, the Quickhull algorithm is a computational geometry algorithm used to compute the convex hull of a set of points in a multi-dimensional space.
The convex hull of a set of points is the smallest convex polygon that encloses all the points in the set. In other words, the convex hull is the outermost boundary of the set of points. The Quickhull algorithm is a divide-and-conquer algorithm that recursively partitions the set of points into smaller subsets and finds the convex hull of each subset.
A divide-and-conquer algorithm is a well-known problem-solving strategy that involves breaking down a problem into smaller, more manageable sub-problems and solving each sub-problem separately. The solutions to the sub-problems are then combined to solve the original problem. Here, the strategy is applied to the geometry problem as known from [36].
Hence, the challenge to numerically model transient (vibro-) acoustics wave propagation phenomena in unbounded domains is still open in the prior art.
The scope of the present invention is defined solely by the appended claims and is not affected to any degree by the statements within this summary.
The present embodiments may obviate one or more of the drawbacks or limitations in the related art. For example, efficiency of numerically modelling transient (vibro-)acoustics wave propagation phenomena in unbounded domains may be improved.
A computer-implemented method includes the additional acts of: (c) constructing a convex shape around the sound source using Quickhull algorithm; (d) extruding infinite elements starting from the convex shape; (e) transforming the coordinate system into a field variable coordinate system by mapping a radial coordinate of the coordinates to an auxiliary coordinate extending from the convex shape to the infinite elements; (f1) constructing the shape functions using the auxiliary coordinate; and (g) determining the field variable by solving the numerical simulation model.
According to the present embodiments, the numerical simulation model concerns numerically modelling transient (vibro-)acoustics wave propagation phenomena in unbounded domains. Therefore, the model may be understood and termed as being a model for ‘numerically modelling transient (vibro-) acoustics wave propagation phenomena in unbounded domains.’
In one embodiment, the field variable is a pressure field. In one embodiment, the field variable may be a pressure field or any derivative of the parameter pressure field.
The basic idea of the present embodiments is an infinite element formulation, which may be referred to as the flexible infinite element, that may be used in combination with the finite element method (FEM) to efficiently model the time domain propagation of sound waves in unbounded media. For example, the new formulation addresses the challenge to make such transient unbounded vibroacoustic simulations very fast and compact, such that the transient unbounded vibroacoustic simulations may be used in Digital Twin scenarios, including applications where the simulations need to run close to real time.
The governing equations describing the physical behavior of sound waves may be written in frequency or in time domain. In frequency domain simulations, the input signal is represented as a sum of harmonic waves. By breaking down the signal into its individual frequency components, the steady-state pressure response at each frequency may be computed separately. Time domain simulations (e.g., transient simulations) model sound waves as variations of pressure over time. For digital twin applications, time domain descriptions are prominent.
The problem of efficiently and accurately simulating the propagation of sound waves in unbounded domains has been widely researched in the past decades. The non-reflection condition may be modelled in various ways with the finite element method (FEM) or, alternatively, with other numerical methods, such as the boundary element method (BEM). The time domain version of BEM does not lend itself to easy building of reduced order models (ROMs). Since such compact representations are needed in digital twin scenarios, the focus of the invention is on FEM.
Since FEM is a domain discretization method, modelling the infinite domains is not inherently possible. The FEM mesh is to be truncated, and non-reflecting boundary conditions are to be enforced on the exterior boundary to prevent the acoustic waves from travelling back into the FEM domain.
The present embodiments present a novel flexible infinite element formulation. This is a unique method that allows the use of general convex-shaped envelopes for infinite elements, while being time-stable; its geometrical flexibility makes it particularly attractive to simulate the sound scattered by elongated structures. The method is also compatible with MOR techniques, because the method has frequency independent system matrices. Therefore, the method also allows compact representations.
Unlike existing formulations that build the infinite elements shape functions based on a standard mapping between the distance from the virtual sources and the radial parent element coordinate, a new auxiliary mapping is introduced between the radial coordinate and a new auxiliary coordinate (see
Therefore, the numerical scheme may provide accurate results, irrespective of the positioning of the mapping nodes in the infinite domain. As such, the new formulation is geometrically flexible, allowing to attach the infinite elements to general convex-shaped envelopes, without compromising the accuracy for transient simulations.
To guarantee the time stability of the flexible infinite element formulation, an accuracy-preserving stabilization procedure has been successfully developed.
According to the present embodiments, the Quickhull algorithm is used to generate the convex acoustic domain around any given scatterer, while the mapping nodes inside the infinite elements are extruded orthogonally to exterior FEM boundary at run time, with the normal directions computed using the same rules as in [12]. Thanks to the adopted strategy, the user intervention is kept to a minimum for the construction of the acoustic and exterior domain, even when no CAD of the scatterer geometry is available.
In addition, as the proposed method results in frequency independent system matrices, the method may be easily combined with MOR methods such as Krylov based methods. This aspect, together with the geometric flexibility of the new formulation, unlocks new digital twin applications in the field of unbounded (vibro-) acoustics.
An embodiment provides that in act (e), the auxiliary coordinate for positions outside of the convex shape is defined as:
where:
This transformation enables a quick and accurate modelling and solving of the acoustics to determine the field variable within a reasonable time.
Another embodiment provides displaying the field variable and/or key figures calculated from the field variable or images illustrating the field variable to a user. By default, the displaying uses a monitor as a human machine interface informing the user about the acoustic behavior of the object.
Another embodiment includes an act of providing the field variable to an iterative object design process of the object for designing the object with the design goal of improving acoustic emission.
Another embodiment includes an act of generating the object according to a design as determined by the iterative object design process. This generation may be done by any manufacturing or production method.
The iterative object design may include as acts of: selecting at least one design parameter of the object (OBJ) design to be improved; providing an iteration abortion criterion as a threshold value for the field variable or a key figures calculated from the field variable; defining at least two variants (VRT) of the object (OBJ) differing in the selected design parameter; modelling the acoustics of the variants (VRT) according to the present embodiments; selecting the best variants from the defined variants based on the modelling results; repeating the acts of defining, modelling, selecting until the abortion criterion is met; selecting the design parameter of the variant that meets the abortion criterion; optionally select at least one different design parameter and processing the same iteration as defined above for the newly selected at least one different design parameter.
This process is robust and efficient. Other iterative object design strategies may be applicable, too.
Further, the present embodiments relate to a computer product arranged and configured to execute the acts of the computer-implemented method according to the present embodiments.
Further, the present embodiments relates to a computer-readable medium (e.g., a non-transitory computer-readable storage medium) encoded with executable instructions that, when executed, cause the computer product defined above to carry out a method according to the present embodiments.
Further, the present embodiments relate to a system for determining acoustic parameters of an emission of the system. The system includes at least one processor being prepared by upload of computer-executable code to perform a method according to the present embodiments.
The finite element method (FEM) is one of the most popular techniques to simulate the propagation of acoustic waves, as it allows to efficiently solve complicated problems, by discretizing the physical domain of interest. The simulation of wave fields in unbounded regions represents a challenge, because computers can only handle a finite number of elements. The simplest, but rather inaccurate, approach consists in applying an impedance boundary condition at the artificial boundary of the finite element domain [1]. High-order absorbing boundary conditions [2] and the perfectly matched layer (PML) [3] represent more accurate alternatives. In particular, the PML has been proven to be very efficient for the solution of wave problems in the frequency domain. The PML has also been extensively studied for transient applications (see [4]-[10] and references therein). However, for time-domain simulations, the PML requires the use of auxiliary variables; the resulting increase in the number of degrees of freedom (Dofs) makes the PML less attractive to simulate transient phenomena.
The infinite element method, introduced by Bettess and Zienkiewicz [11], avoids the issue of truncating the computational domain by employing special elements, which extend from the boundary of the FE domain to infinity. Over the years, different formulations have been proposed in the literature. In particular, two groups of formulations, which differ in the choice of the test function in the infinite domain, were developed: the unconjugated, where the test basis functions and the shape functions are the same, and the conjugated formulation, where the test basis functions are the complex conjugate of the shape functions. A detailed comparison of these different strategies is given in [12], where the advantages and disadvantages of both are analyzed in detail. The Astley-Leis formulation [13] differs from the classical conjugated formulation only by a weight factor that multiplies the test function, which is designed to cancel out the unbounded terms in the integrals that need to be computed to obtain the system matrices.
Alternatives exist also in terms of geometrical discretization of the infinite elements. The so-called mapped formulations rely on a mapping from the physical element to a parent element [13]-[16]: their main advantage is the geometrical flexibility, as they can be attached to arbitrary convex-shaped FE boundaries. By contrast, for Burnett elements [17]-[19], the geometry of the envelope is restricted to spheroids and no mapping is employed.
Astley-Leis elements have several attractive properties. First, the system matrices integrals can be evaluated using standard quadrature rules, thanks to the weight factor used in the test function. Moreover, they result in frequency-independent system matrices and maintain the sparsity of the FEM system, which makes them readily available for time-domain simulations [20]-[22] and suitable for Model Order Reduction (MOR) [23]. This is particularly interesting, as in recent years MOR for transient simulations has emerged as one of the key technologies in a number of applications [24]-[27]. The aim of MOR techniques is to build Reduced Order Models (ROMs) which capture the most relevant dynamics of the corresponding initial Full Order Models (FOM), using only a fraction of the Dofs. In [26] a numerical test is presented, for which the performance of the transient formulation from Kaltenbacher et al. [10] is compared with the Astley-Leis infinite element for mulation. In particular, it is shown that using the PML results in a larger number of Dofs, due to the presence of the auxiliary variables. Moreover, from a MOR perspective, it is necessary to use a split-basis approach, which results in tripling the number of Dofs of the final ROM, to guarantee time-domain stability preservation of the PML formulation after MOR. Since this performance comparison between the PML and the infinite element method is limited to a single numerical test, it is difficult to draw general conclusions. The study cited, however, seems to give a strong indication that using infinite elements for unbounded transient problems may be the best option, especially when MOR techniques are used.
Astley and Hamilton [28] showed that the Astley-Leis formulation is time-stable, provided that the infinite edges are aligned with the normal of the interface between the FEM and the exterior domain. However, the performance of traditional mapped infinite elements varies with the position of the virtual sources [14]; the mentioned stability requirements do not allow to freely choose the location of the virtual sources, and this often represents a problem for accuracy in transient analysis, when arbitrarily convex-shaped interior domains are analyzed. In practice, for most scenarios, this means that one should attached mapped Astley-Leis infinite elements only to spherical envelopes for transient simulations. Alternative mapped spheroidal infinite element formulations, based on the Holford expansion [29], have been proposed, both in the frequency and in the time domain [30], [31]. However, for such formulations, as well as for non-mapped spheroidal formulations [17]-[19], the benefits of a smaller interior domain may be nullified by the necessary increases in the radial order of the infinite element shape functions, as observed in [31], [32] and [33].
Li et al. [34]proposed an alternative conjugated mapped formulation, where radial shape function of the Lagrangian type are modified, to satisfy exactly the truncated multipole expansion along the infinite edge, irrespective of the location of the virtual sources. However, they only addressed frequency-domain problems, where the standard Astley-Leis formulation already provides satisfactory results; moreover, being limited to radial shape function derived from Lagrange polynomials, Li's elements have bad conditioning for high radial orders [35].
In this work, we propose a novel infinite element formulation for radial shape functions of any type, which satisfies exactly the truncated multipole expansion, irrespective of the orientation of the infinite edges. Consequently, we employ the new flexible infinite element to achieve a time-stable and accurate formulation, for arbitrarily convex shaped envelopes. The integration of the flexible infinite elements into an high-order software is discussed. A strategy is presented, to automatically extrude the infinite elements, when only a mesh of the interior domain is available. Additionally, we show how one can take advantage of the quickhull algorithm [36], implemented in the simulation package Simcenter 3D [37], developed by Siemens Industry Software, to automatically generate the envelope where the infinite elements are attached, around the scattering object. Finally, the combination of the proposed formulation with an existing MOR strategy is discussed.
We take the convention that the time-dependence of the fields is eiwt, where w is the angular frequency and t is the time. The behavior of the pressure field p in the unbounded domain Ω is described by the Helmholtz equation, in the frequency domain:
where k is the wavenumber. The domain Ω is subdivided into an interior region Ωi and an exterior region Ωe; we denote with Γ the interface between the two domains (see
At the interior boundary ΓN of the domain Ω a Neumann boundary condition is prescribed. Dirichlet and Robin boundary conditions can be easily handled as well, but won't be considered here for ease of notation. The Sommerfeld condition states that no waves can be reflected at infinity; in mathematical terms, in a d dimensional domain, this writes as:
where r is the spherical/cylindrical radius.
As in [38], approximating (2), we can write:
The corresponding boundary value problem reads as follows:
where ρ is the fluid density and vn is the prescribed normal velocity on ΓN.
Using the method of weighted residuals, we obtain the variational formulation:
where H1 is the Sobolev space where the solution is sought.
Eq. (5) is solved by using H1-conforming finite elements in Ωi and attaching infinite elements to the envelope Γ. The field variable p and the test function q are approximated by phϵVh⊂H1(Ω) and qhϵQh⊂H1(Ω) where Vh is the space of shape functions and Qh is the space of test functions. In the Astley-Leis formulation, thanks to the weight factor used in the test function, the Sommerfeld condition is automatically satisfied and does not need to be taken into account in the variational formulation [39]. In this case, the discrete variational formulation is:
which can be rewritten as a linear system of equations:
where p is the vector of Dofs, f is the forcing vector and M, C and Kϵn
Since, for the Astey-Leis formulation, the system matrices are frequency-independent, Eq. (7) can be readily transformed into the time domain:
where pτ and fτ are the inverse Fourier transforms of p and f, respectively.
The Atkinson-Wilcox theorem [40] states that at any point outside the smallest sphere circumscribing all the sound sources, the solution of problem (4) can be written as the multipole expansion:
where Gn is the directivity function, α and θ are the spherical angles and τ is the radial distance, measured from the center of radiation, that coincides with the center of the minimal sphere. For ease of notation, throughout the rest of the paper, we will assume that the origin O of the absolute Cartesian coordinate system coincides with the center of radiation.
The Atkinson-Wilcox theorem has been extended by Holford 1291, to include the spheroidal case. By placing the envelope Γ outside the minimal sphere (for the Atkinson-Wilcox expansion) or spheroid (for the Holford expansion), we ensure that the solution in the exterior domain Ωe is composed only of outward propagating waves: this consideration is used to select a proper trial solution in the infinite elements.
The so-called spheroidal formulations have received some criticism in the past [32], [33], since the smaller number of Dofs in the FEM region is often counterbalanced, in practical applications, by the need for higher radial orders in the exterior region. On the other hand, numerical experiments have been performed in the frequency domain, using mapped infinite elements based on expansion (9), where the conditions of the Atkinson-Wilcox theorem were explicitly violated [41], [35], [42]: namely, non-spherical envelopes were used, partially lying inside the minimal enclosing sphere. Although the convergence of this approach has not been formally proven [43], numerical experiments with sufficiently smooth envelopes gave satisfactory results (see [44] for more details). For all the test cases encountered during our research, our observations were in line with the cited literature. These claims will be further expanded in Section 6.
In this subsection, we describe the main features of the Astley-Leis formulation and we analyze under which conditions it is stable in the time domain.
The Astley-Leis formulation relies 313 a mapping from the infinite physical element to a finite parent element. In particular, in two dimensions (see
In
Similarly, we write the Cartesian coordinates of the mapping node Γie on Γ as:
The i-th infinite edge is the line originating from node Γie, passing through node Σie and extending to infinity (see
We denote with
where ∥⋅∥ is the Euclidean norm. The virtual source Ōie, with Cartesian coordinates:
is located (see
As a consequence the distance
The envelope Σ is the union of all points, in the physical space, whose parent element radial coordinate is
Notice that, with this mapping definition, we have
The interpolated distance
where n is the number of mapping nodes on the envelope Γ, in the e-th infinite element, and Nti are the standard tangential geometrical mapping functions (in our case, we use Lagrangian polynomials). Similarly, the interpolated distance
By substituting eq. (15) in (19), we obtain:
By then substituting (18) in (20), we obtain:
The interpolated distance between the mapping nodes in the element e and their corresponding virtual sources is defined as:
where N
By substituting (23) and (21) in (22), we finally obtain:
In other words, a mapping between the physical coordinate
This mapping is consistent with the notion that points with parent element radial coordinate
The vector of Cartesian coordinates for a point inside the c-th infinite element is found by interpolating the nodal Cartesian coordinates:
where xΓ
We emphasize that the notation adopted ere slightly different from what is usually found in the literature, because we want to distinguish between r as defined in Section 2.1.3 (e.g. the distance from the center of radiation) and
Additionally, we define rie=rie(
Finally, the so-called phase distance
When the infinite edges are extruded orthogonally to Γ, the phase distance
By analogy with the Atkinson-Wilcox expansion (9), the field shape functions are written as:
where ψ(s, t,
By inserting the radial mapping (24) in (31), it can be verified that the following relationship holds:
Extending the above argument to higher order polynomials in three dimensions, we obtain:
where m is the selected radial order in the infinite domain. In particular, along the i-th infinite edge, we have
In other words, for the i-th infinite edge, the pressure field generated by the m-th order virtual multipole source, placed at Ōi, belongs to the space of the shape functions. This is where the idea of virtual sources come from.
An interesting case to analyze is when all the virtual sources of the e-th infinite element coincide with the center of radiation O. In this case, along the i-th infinite edge, we obtain:
since
even when all the virtual sources coincide with the center of radiation. This comes from the fact that in (22)
In the Astley-Leis formulation, a Petrov-Galerkin approach is taken: the test function and the solution function belong to different function spaces; in particular, the test basis functions are the complex conjugate of the shape functions, multiplied by a weight factor w:
where the weight factor takes the form:
In contrast to the Bubnov-Galerkin method, where the test basis functions are equal to the shape functions, the Petrov-Galerkin method results in un-symmetric system matrices. However, this choice for the test basis functions is particularly convenient: firstly, being a conjugated formulation, the Astley-Leis formulation results in frequency-independent system matrices, which can be then readily used in the time domain; secondly, for pw>2, standard quadrature rules can be used for the evaluation of the integrals necessary to compute the system matrices. In [39] the choice of pw=2 is justified on physical grounds, although it is pointed out that any power pw>2 preserves the boundedness of the system matrices integrals. Irrespective of the choice of pw, the conjugated weighted formulation was shown to lead to a variational formulation in weighted Sobolev spaces, for which existence and uniqueness of the solution has been proven [12], [45].
In this subsection, we analyze the conditions under which the Astley-Leis formulation is stable. This analysis is an extension to the theory developed from Astley and Hamilton [28]. In particular, in Section 3.5 of [28] only the system matrices relative to the infinite elements were considered; here, by contrast, we include the case where the interior Doman is discretized with finite elements. Moreover, in Section 3.5 of [28] only the case where the inner surface of the infinite element is an are of a circle was inspected; in this subsection, we extend the analysis to arbitrary convex shapes. Additionally, to guarantee stability, we will show that it is not necessary to construct “zero-mass” elements, but it is sufficient to require that the element mass matrices relative to the infinite domain are positive semi-definite.
As stated in [46] the second-order system (8) is stable, provided that the damping, stiffness and mass matrices are positive semi-definite. A matrix Aϵn
For the sake of brevity, we will adopt A≥0 as a shorthand notation to indicate positive semi-definiteness of the matrix A. Because of the infinite elements' contribution, the mass matrix may not always be positive semi-definite, causing the system to become unstable.
First of all, we show that, if the mass matrix relative to the infinite domain, MIE ϵn
n
where the subscripts FE and IE denote the Dofs associated with the finite and the infinite elements, respectively; the last inequality in (40) comes from the fact that MIE≥0 by hypothesis and MFE≥0 when standard polynomial shape functions are used and the material properties of the medium in which sound propagates are physically realizable.
Next, we show that MIE≥0, provided that the element mass matrix MIEeϵn
n
where nE is the number of infinite elements and the inequality comes from the fact that MIEe≥0 by hypothesis, with e=1, . . . , nE.
We now analyze under which conditions MIEe≥0. In the Astley-Leis formulation, when the test basis functions are taken as the weighted complex conjugate of the shape functions MIEe takes the form (see [13]):
where e is the speed of sound. To make it easier to analyze its definiteness, we rewrite it to a much simpler form. In particular, when using Gauss quadrature formulas to perform the integration, expression (42) can be rewritten as
where Ψgϵn
ng is the number of Gauss points, gj=(sjg, tjg,
where gjw is the j-th Gauss weight and det Jj is the determinant of the Jacobian of the transformation between physical and parent element coordinates, evaluated at the j-th Gauss point.
Since De is a diagonal matrix, it is positive semi-definite if and only if all its components are not smaller than zero. This happens if and only if
We show that if De≥0 then MIE≥0. Let's choose a generic vector yϵn:
where the last inequality in (47) comes from the fact that De≥0 by hypothesis, with z=Ψgyϵn
Let's define the vector
where ne(s, t) is the unit vector normal to the envelope Γ at the point identified by the parent element coordinates (s, t,
Additionally, for the properties of the gradient operator, thanks to the relationship between the gradient of a function and its directional derivatives, we have:
Let's consider first the ideal scenario where the envelope Γ is represented exactly, and no error arises from the geometrical discretization, and let's choose the infinite edges in the e-th element to be extruded orthogonally to the envelope Γ. In this case,
Finally, inserting (50) and (49) in (51) and simplifying the resulting equation, we obtain
This means that, if the extrusion direction is orthogonal to the envelope Γ for all the infinite edges, we obtain ∥≡9x
Requiring that the infinite elements are extruded orthogonally to Γ means that, in general, the virtual sources cannot be located at the center of the smallest sphere enclosing all the sound sources. This limits the accuracy of all standard mapped formulations, especially for envelopes with an high aspect ratio, and explains the need to develop a new flexible infinite element formulation.
In this section, we describe the new flexible formulation, for which the geometric and field shape functions are constructed using two distinct mappings. This feature makes the accuracy of our formulation independent of the geometry of the infinite edges. As a consequence, the flexible infinite elements can be used to accurately simulate transient acoustic wave propagation problems, for arbitrarily convex-shaped envelopes.
We recall that we denote with r the distance from the center of radiation and with
everywhere inside the infinite elements, irrespective of the location of the virtual sources. To do that, we first redefine the phase distance as
where re is computed according to (28), and rΓe is the radial distance evaluated on Γ:
One may try to write the shape functions as
where ψr(s, t, re(s, t,
It should be noticed that this mapping largely resembles (25), with v, re and rΓe instead of
Exploiting the auxiliary mapping, we propose a new formulation for the shape functions:
where ψ are the usual polynomials, expressed with respect to the variables s, t, v instead of the variables s, t,
We emphasize that the only difference between (31) and (59) is that U gets substituted with the auxiliary variable v. By inserting the radial mapping (57) in (59), it can be verified that the following relationship holds:
Extending the above argument to higher order polynomials in three dimensions, we obtain:
where m is the selected radial order in the exterior domain. This is true not only for points along the infinite edges, but for any point inside the infinite element.
In other words, the multipole expansion (9) truncated at order m, belongs to the space spanned by the columns of the m-th order shape functions. The concept of virtual sources does not apply to our formulation, as relationship (61) is true, irrespective of the geometry of the infinite elements. This feature makes our formulation more flexible than traditional mapped formulations, since the accuracy does not depend on the positioning of the mapping nodes. As in the Astley-Leis formulation, the test functions are chosen as the weighted complex conjugate of the shape functions:
The new flexible infinite element requires a separate stability analysis, since the definition of μe in (54) differs from the definition of
First, we analyze the behavior of ∥∇xμ∥, showing that, when the infinite edges are orthogonal to the envelope Γ, ∥∇xμ∥≤1 in the near-field, while ∥∇xμ∥≥1 in the intermediate and far-field. Then, we propose a strategy to stabilize the time-domain formulation without compromising the accuracy.
3.2.1. Behavior of ∥∇xμ∥
In this subsection, we study the behavior of ∥∇xμ∥, when the infinite edges are extruded in the direction orthogonal to the envelope Γ, as in the time-stable infinite elements proposed in [28]. For simplicity, we consider the two-dimensional case, but all the considerations we make hold in three dimensions.
Let's consider a point on the envelope Γ, with parent element coordinates (t,
The envelope Γ is an iso-line of μ, with μ=0 by definition. It follows
since the gradient of μ is perpendicular to the iso-lines of μ. For the standard Astley-Leis formulation, we have already shown that relationship (52) holds. Let's define
The partial derivatives of
By triangle inequality (see
from which follows
As a consequence, we finally obtain
Once again, as in the classical Astley-Leis formulation, small errors may arise because of the geometrical approximation. In this Case, as in [28], we propose to neglect these small discrepancies, to ensure the stability of the method.
For
In the region between
Because of the behavior of ∥∇xμ∥, analyzed in the previous subsection a strategy to stabilize the flexible infinite element formulation is necessary. In the non-stabilized version, the mass matrix contribution at the element level writes:
To stabilize the new flexible formulation, we propose to substitute the non-stabilized mass matrix MIEe with {tilde over (M)}IEe, where
with
When computing (74) using Gauss quadrature, we can write
with the j-th component of the diagonal matrix {tilde over (D)}e
This ensures that {tilde over (M)}IEe≥0 and thus that the formulation is time-stable.
As an example, we consider a two-dimensional case, where the FE boundary is an elliptical envelope, with a width of 6 and a height of 2. In
In the previous subsection, we showed how to stabilize the flexible infinite element formulation. However, the stabilization procedure may potentially affect the accuracy of the method. This comes from the fact that some of the entries in {tilde over (D)}e may differ substantially from the corresponding entries in De, making {tilde over (M)}IEe a bad approximation of MIEe.
Exploiting relationship (71), which holds when the infinite edges are extruded orthogonally to the envelope Γ, we can design a strategy to ensure the accuracy of the flexible infinite elements after the stabilization procedure. We propose to tune the power pw in the weight factor w, to make sure that the stabilization procedure does not significantly affect the accuracy. In practice, this is achieved by controlling the relative difference between {tilde over (D)}e and De: the sa tiller this difference, the smaller the accuracy loss. We recall that the weight factor in the Astley-Leis formulation is (38), with pw=2. We observe that w=1 at the envelope and w→0 as
while the corresponding entry in {tilde over (D)}e is
We emphasize that, since relationship (71) holds for orthogonally extruded infinite edges, gκ cannot lie on the envelope Γ. Fort this reason, it is guaranteed that w(gκ)<1 and we can impose
by choosing
where ϵ is a tolerance value, set by the user. The parameter pw is chosen constant across the whole domain, to guarantee the continuity of the test functions across the whole domain, such that (81) is satisfied for all the infinite elements. A smaller tolerance t implies a smaller difference between De and {tilde over (D)}e and, as a consequence, a smaller accuracy loss after the stabilization step. We emphasize that stability is guaranteed by (75), while the accuracy of the stabilized formulation is guaranteed by (81). In Appendix A, we present a bound for the accuracy loss due to the stabilization procedure, linking it to the chosen tolerance ϵ. The oily cost of selecting pw>2 is that an higher number of Gauss points are needed to carry out the computation, since the order of the polynomials to be integrated increases with pw. Typically, for numerical acoustic simulations using the FEMI, the most computationally intensive operation is the factorization of the fully assembled matrix; for this reason, despite an increase in the computational cost of the assembly step, the total computational cost of the flexible formulation is comparable to the Astley-Leis formulation, for the same number of Dofs. Additionally, in many cases, as we will document in Section 5 and 6, the flexible infinite elements result in a drastic decrease in the number of Dofs, thanks to their geometric flexibility, making the simulations less computationally intensive.
As already explained in previous sections, the infinite edges need to be extruded orthogonally to the envelope T to ensure stability in the time domain. Unfortunately, it is not trivial to compute the normal of the envelope Γ for general convex shapes. In some cases, the information about the normal vectors could be provided by CAD tools, with the drawback of requiring an additional preprocessing step. In other cases, only meshed domains are available, without a reference geometry. We solve this issue by following the workflow presented in [47]. With this approach, only a mesh of the internal domain is necessary. The required normal vectors at the mapping nodes on Γ are computed as follows:
The process is repeated for all the mapping nodes on the envelope Γ, for all the infinite elements.
In this subsection, we discuss the implementation of the infinite elements into an high-order FEM code. The considerations made here hold both for the Astley-Leis formulation and the flexible infinite element. We refer to the adaptive order p-FEM approach described in [48], where Lobatto shape functions (also called integrated Legendre polynomials) are used, and to [49], where the method was generalized to include anisotropic orders. For a number of Helmholtz problems, tis approach has been shown to drastically improve the efficiency of the numerical models, when compared with conventional FEM, exploiting the fact that high-order polynomial approximations are more effective at controlling the pollution effect. An added benefit of using a hierarchical set of shape functions, such as the Lobatto polynomials, is found when solutions in the frequency domain are required over a range of frequencies. It is possible to compute the element matrices only once for the highest order, and then, at each frequency, one can extract only the required portion of these matrices to assemble the global matrix. A key feature of the proposed method is the use of a simple local error indicator to select a priori the polynomial order in each element. For calculations over a range of frequencies this approach requires only a single mesh, and the error estimator is able to select the suitable polynomial order in each element to maintain the efficiency and accuracy as the frequency is changed, by identifying the lowest order required for the accuracy target.
Infinite elements based on Lagrange polynomials have been shown to present conditioning problems for high radial orders [35]. The use of Jacobi [35] and Berstein [50] polynomials has been proposed to improve the conditioning and optimize the performance of Krylov-based iterative solvers. As pointed out in [35] Lobatto polynomials, even though sub-optimal from a conditioning perspective, outperform Lagrange polynomials in this respect.
By adopting the same type of polynomials for finite and infinite elements, the software implementation becomes particularly easy. In this case, the infinite element subroutine can be developed by reusing most of the routines called by its finite element counterpart. The only differences are that the new radial geometric mapping (23) needs to be implemented, and different formulas nee to be used for the computation of the system matrices, which does not require any additional implementation effort. In other words, in a software implementation, the infinite elements can be simply seen as specialized anisotropic high-order elements. This is an important aspect, as the implementation burden of the infinite element method may hinder its broad Use in practical applications. An additional benefit of using Lobatto polynomials also in the exterior domain is that the hierarchical structure of the system matrices is preserved for unbounded-domain simulations, and can be exploited for frequency-domain calculations.
The shape functions orders inside the FE domain are selected based on the a priori error estimator described in [48]: the radial order on the infinite edges is selected by the user, while the orders of the shape functions on the envelope Γ are designed to ensure continuity with the FE domain; finally, vertex, edge and face Dofs are set to zero at infinity (
In this section, the flexible infinite element is compared to the Astley-Leis formulation, with numerical examples in two and three dimensions, both in the frequency and in the time domain. The FE domain is meshed with the mesh generator Gmsh [51], while the infinite elements are extruded at run time, using one of the following strategies:
In the following examples, we will often refer to the relative L2-error between the numerical and the reference solution, unum and uref respectively, computed on the envelope Γ:
and in the interior domain (the FE domain):
When using the flexible formulation for transient problems, we will mea-sure the error induced in the frequency domain by the stabilization procedure (see discussion in Section 3.2.3), for the largest frequency of interest, as:
where unumst is obtained using the stabilized formulation, while unumunst is obtained using the non-stabilized formulation.
First, we analyze the case of a monopole source in a 2D domain, for which the analytical solution is known:
where H20 is the Hankel function of the second kind and the monopole source location coincides with the center of radiation. The source is enclosed by a convex envelope and the two meshes generated using the confocal-rays and the normal-rays approach are shown in
For the following frequency-domain results, the order of the shape functions is fixed (both in the interior and the exterior domain). Numerical convergence studies are done by investigating the effect of the order chosen on the accuracy of the methods. First of all, we compare the performance of the new flexible infinite element against the Astley-Leis formulation, for the confocal-rays configuration. Even though, for the configuration considered, all the virtual sources coincide with the center of radiation, the Astley-Leis formulation and the flexible formulation are still different. This is because, for the Astley-Leis formulation, the radial coordinates are found, as in (18) and (22) using the standard interpolating functions, while, for the flexible infinite element, they are computed as a function of the Cartesian coordinates using (28) and (55). In
For all four frequencies, the L2-error of the standard Astley-Leis formulation reaches a plateau, while the flexible infinite element approach displays better p-convergence properties. This results are explained by the fact that (53) is satisfied everywhere for the flexible infinite elements, while for the Astley-Leis formulation it is only valid along the infinite edges, when the confocal-rays extrusion approach is adopted (see discussion in 2.2.2). In other words, the truncated multipole expansion (9) belongs to the space spanned by the shape functions in the whole exterior domain for the flexible formulation, but only along the infinite edges for the Astley-Leis formulation. The resulting accuracy difference is higher for coarser meshes, which are typically usec in the context of high-order adaptive FEM, while it is negligible when fine meshes are used.
Next, we compare the flexible infinite elements with the Astley-Leis infinite elements, when the normal-rays extrusion strategy is adopted. This serves as a preliminary to the time-domain comparison, where the normal rays strategy is the only one capable of delivering stable results. In
It is well known that the confocal-rays extrusion approach leads, in general, to unstable formulations in the time domain, as discussed in Section 2.3. For this reason, we focus only on the normal-rays approach. For this example, the monopole source is placed it (0.5, 0.4), to show that the flexible infinite elements work well even when the center of radiation does not coincide with the physical source. The input signal used for the simulation is a sine-wave, with k=2, filtered using the Hamming window. So, the forcing vector fτ(τ), with τϵ[0, T], is
where f is the forcing vector in the frequency domain. The order of the shape functions in the FEM domain is selected adaptively, using the a priori error indicator from Bériot et al. [48], with target L2-error E2=0.01, while the radial order in the infinite domain is set to m=6 for both cases. The power in (38) for the flexible formulation is set to pw=6 and the error induced by the stabilization procedure is E∞st=0.0007, computed according to (85) for K=20, the largest wavenumber of interest for the problem analyzed. In other words, it is clear that the error induced by the stabilization procedure is negligible, since it is largely below, the target error used to compute the shape functions orders in the FEM domain.
A reference numerical solution is computed on a larger circular FEM domain, centered at the monopole source position, with radius 2, using the Astley-Leis infinite elements with the infinite edges extruded in the radial direction, with extrusion length
In
In this subsection, we consider the scattering of a plane wave by a sphere and we test the performance of the flexible infinite elements and of the Astley-Leis infinite elements. Several polyhedral shapes, with rounded edges, were used as enclosing surfaces. Rounding the corners is necessary to guarantee the accuracy of the method, regardless of the infinite element formulation used, since the normal vector has to be continuous throughout the whole surface, for the extrusion direction to be orthogonal to the envelope at any point.
The scattering sphere has unit radius, while the pseudo-polyhedral domains have midradius 2 aid fillet radius 0.2. The mesh used in the internal domain has elements with a characteristic length h=0.3, with a refinement corresponding to the rounded corners, to resolve the geometric details, and the order in the FEM region is once again set adaptively, using the a priori error indicator from Bériot et al. [48], with target L2-error E2=0.01
The reference analytical solution for this problem is:
where asca is the radius of the scattering sphere, jm is the m-th order spherical Bessel function, hm(1) is the m-th order first-kind spherical Hankel function, and Pm is the mth-order Legendre polynomial.
The best-interpolation error is computed as the relative L-error between the analytical solution and its L2-projection onto the finite element space:
This error corresponds to the best numerical solution that can be reached, given the FE mesh and the order used, irrespective of the infinite element formulation used. In
In this subsection, out of all the the pseudo-polyhedral shapes presented in Section 6.3.1, only the most challenging one, the pseudo-tetrahedron, was considered. We now consider the results from transient simulations, where the input signal is a sine-wave with k=10, filtered by the Hamming function, as in the 2D case. In
As in the two-dimensional case, the reference numerical solution is computed on a larger FEM domain, consisting of a sphere with radius 2 surrounding the scatterer, using the Astley-Leis infinite elements with the infinite edges extruded in the radial direction, such that all the virtual sources coincide with the center of the scattering sphere. The radial order in the infinite domain for the reference solution is set to 10 and, once again, the radial quantities are computed according to expressions (28) and (55). The reference numerical solution in the the domain is computed constructing “zero-mass” infinite elements, following the procedure from Astley and Hamilton [28]. As in the two-dimensional case, good accuracy of the reference numerical model is guaranteed by the use of a large FEM domain, high radial order m and the modified interpolation of radial quantities, combined with the fact that all the virtual sources coincide with the center of radiation.
Also for this case, it can be noticed that the flexible formulation outperforms the traditional Astley-Leis formulation.
In this section, we present a workflow to tackle unbounded problems of industrial complexity, using the flexible infinite element formulation. In particular, to showcase the effectiveness of the method, we will consider the shark submarine, shown in
In this subsection, we describe the automatic process to get a mesh for the computational domain. First, a surface mesh made of second-order triangular elements, with characteristic length h=0.02, is created on the submarine. Secondly, an approximate convex hull is generated automatically using the quickhull algorithm [36], implemented in the simulation package Simcenter 3D [37] developed by Siemens Industry Software. The minimal distance from the submarine surface to the convex envelope is d=0.0. Once again, quadratic triangular elements with characteristic length h=0.02 are used to mesh the external surface. Then, a volume mesh, that comprises tetrahedral quadratic elements, is created between the two surface meshes. The infinite elements are extruded automatically along the normal direction, so that the transient formulation is time-stable, following the process described in Section 3.3. The final mesh, shown in
Since the flexible infinite elements, as the traditional Astley-Leis infinite elements, result in frequency-independent and sparse system matrices, their use in combination with Krylov-based MOR is particularly attractive. To accelerate transient simulations, we use the AKSA [2], [26] algorithm, which aims at building automatically a reduced order model, whose frequency response function (FRF) matches the FRF of the initial full order model (FOM), within a user-specified tolerance, in the frequency range of interest, by only using a few system solves and enriching the projection basis with higher-order moments. To obtain a time-domain reduced model of the shark submarine, for example, valid between k=5 and k=50, with a relative error tolerance tol=0.01, only 5 full system solves were necessary.
We now test the performance of the flexible infinite element for the scattering of plane waves in the frequency domain. The simulations were performed for four plane waves, traveling at different angles of incidence a and θ. As a reference, we use the results obtained using the automatic perfectly matched layer (AML) implementation [47], with 5 layers. For the configuration using infinite elements, the order in the FEM domain is assigned adaptively, using the a priori error estimator front Bériot et al. [48], with target L2-error E2=0.01, while the radial order in the exterior domain is user-defined; for the AML case, the order is assigned adaptively in the whole computational domain. For all cases analyzed, the mesh used is the one shown in
Finally, we consider the scattering of a monopole source placed in front of the shark, at (0, 0, 0.5).
The target L2-error, used to compute the shape functions orders in the FEM domain, is set once again to E2=0.01. The radial order for the infinite elements is set to m=8, fixed throughout the exterior domain. The power in (38) for the flexible formulation is set to pw=8 and the error induced by the stabilization procedure is E∞st=0.0002, computed according to (85) for k=50, the largest wavenumber of interest for the problem analyzed. Once again, it is clear that the error induced by the stabilization procedure is negligible, since it is largely below the target error used to compute the shape functions orders in the FEM domain.
The input signal consists of a sine wave with k=45, filtered using the Hamming window. A ROM, valid between k=5 and k=50, with a relative error tolerance tol=0.01, is generated using the AKSA. The FRFs of the FOM and the ROM are compared in
Dofs. We ran all our numerical experiments on a Dell desktop computer with 384 Gb of RAM and 3.00 GHz clock speed. The offline MOR procedure took 5 minutes and 45 seconds, and a single full-system solve, an operation performed several tines during the AKSA, required 15.5 Gb of memory. As a reference, we compute the transient response to the same excitation, using the standard Astley-Leis formulation, with a spherical envelope surrounding the submarine. In this case, the model has 1 039 249 nodes, 801 405 elements and 2 515 697 Dofs, 5 times more than the convex-hull model. The offline MOR procedure took 95 minutes, considerably more than the previous case, and for a full-system solve 151.9 Gb of memory are necessary; additionally, the projection basis, that needs to be stored throughout the whole MOR process, is 5 times larger than before, as the number of rows corresponds to the Dofs of the FOM. As in the previous case, the final ROM has 90 Dofs. In Table 1 we summarize the data on the computational requirement to build a ROM using the AKSA for the two models.
In
In this work, we introduced a novel infinite element formulation, called flexible infinite element, for the simulation of unbounded wave propagation problems. One of the main disadvantages of traditional mapped formulations is that the accuracy depends on the geometrical configuration of the infinite elements. For transient simulators, where the infinite elements have to be extruded orthogonally to the FE boundary, the user is not free to locate the virtual sources where desired. This limits the accuracy of such formulations, especially in the case in which envelopes of a high aspect ratio are used. In the flexible formulation, a new set of shape functions is used, that makes the accuracy independent of the location of the mapping nodes. The newly developed element can be used for transient simulations, with arbitrarily convex-shaped FE boundaries. A strategy to guarantee the stability of the proposed method in the time domain, and at the same time preserve the accuracy, was discussed in detail. The accuracy of the flexible formulation was tested with several examples, both in the frequency and in the time domain. It was shown how the new formulation can be easily implemented in a software based on high-order adaptive FEM. Finally, a workflow including the automatic extrusion of the infinite elements and the autocratic order reduction of the resulting models was presented. Since typically higher order polynomials are used for the test functions, to guarantee the accuracy of the stabilized transient flexible formulation, the cost of the assembly step is higher than for standard Astley-Leis infinite elements (see discussion in Section 3.2.3). However, for most cases of interest, the most computationally expensive operation is the factorization step. In particular, the computational savings made possible by the use of the flexible infinite elements were demonstrated, especially for high-aspect ratio scatterers. For future research, it would be interesting to develop a new error estimator, to automatically determine the radial order of the shape functions in the infinite domain. Finally, it would be interesting to research new strategies to automate the choice of the power pw in the test functions for transient simulations, based on accuracy requirements from the user.
In this appendix, we analyze the effects of the stabilization procedure on the accuracy of the solution P of the linear system of equation (7). In particular, we develop an error bound, which links the accuracy loss to the chosen tolerance c in equation (80). Let's call Afull the fully assembled matrix
and let the perturbed system be
where ΔAfull is the perturbation in Afull, introduced by the stabilization procedure, and Δp(ω) is the subsequent perturbation in the solution vector As stated in [52] we have
where the ∞-norm of a matrix Aϵa×b is
and the condition number κ∞ is defined as
Since the stabilization procedure consists in substituting the infinite elements' mass matrix contribution MIE with the stabilized version {tilde over (M)}IE, we have
It is interesting to note that, since ∥ΔAfull∥∞ is proportional to ω2, the accuracy loss is larger for the highest frequency of interest. At the element level, we can write
The maximum entry of the matrix De-{tilde over (D)}e is at most ϵ, which is imposed through (80) and (81). As a consequence, the maximum entry of the matrix MIEe-{tilde over (M)}IEe is smaller than ng2ϵ, since Ψgϵn
where nΓ is the number of mapping nodes on the envelope Γ and nelemi is the number of elements which share the i-th node. So, we have that the maximum entry in the matrix MIE-{tilde over (M)}IE is at most nmaxng2ϵ. Using the ∞-norm definition, we can write
with
where nnzi is defined as the number of non-zero entries in the i-th row of the matrix MIE. By substituting inequality (A.9) in (A.6), we obtain
which can be substituted in (A.3), to obtain
Expression (A.12) represents the link between the tolerance ϵ, set by the user at the element level, and the accuracy loss coning from the stabilization procedure.
Additional possible implementations or alternative solutions of the present embodiments also encompass combinations of features described above or below with regard to the embodiments that are not explicitly mentioned herein. The person skilled in the art may also add individual or isolated aspects and features to the most basic form of the invention.
It is part of the invention that not all steps of the method necessarily have to be performed on the same component or computer instance but may also be performed on different computer instances.
The properties, features, and advantages of the present embodiments described above, as well as the manner they are achieved, become clearer and more understandable in the light of the following description and embodiments, which will be described in more detail in the context of the drawings. This following description does not limit the invention on the contained embodiments. Same components or parts may be labeled with the same reference signs in different figures. In general, the figures are not for scale. An embodiment may also be any combination of the dependent claims or above embodiments with the respective independent claim.
These and other aspects of the invention will be apparent from and elucidated with reference to the embodiments de-scribed hereinafter.
Embodiments are now described, by way of example only, with reference to the accompanying drawings, of which:
The illustration in the drawings is in schematic form.
It is noted that in different figures, similar or same elements may be provided with the same reference signs.
The method may be implemented as or stored on or downloadable as a computer product CMP arranged and configured to execute the acts of the computer-implemented method according to the present embodiments. The computer product CMP may be computer hardware including executable instructions stored in a data storage. An additional or alternative option is to provide the method according to the present embodiments via a computer-readable medium CRM encoded with executable instructions that, when executed, cause a computer or computer product to carry out a method according to the present embodiments or one of its embodiments. Another option is providing the present embodiments as a system SYS for determining acoustic parameters of an emission of the object. The system SYS includes at least one processor CPU being prepared by upload of computer-executable code to perform a method according to the present embodiments.
Although the present invention has been described in detail with reference to embodiments, it is to be understood that the present invention is not limited by the disclosed examples, and numerous additional modifications and variations may be made thereto by a person skilled in the art without departing from the scope of the invention.
The use of “a” or “an” throughout this application does not exclude a plurality, and “comprising” does not exclude other steps or elements. Also, elements described in association with different embodiments may be combined. Reference signs in the claims should not be construed as limiting the scope of the claims.
Independent of the grammatical term usage, individuals with male, female, or other gender identities are included within the term.
The elements and features recited in the appended claims may be combined in different ways to produce new claims that likewise fall within the scope of the present invention. Thus, whereas the dependent claims appended below depend from only a single independent or dependent claim, it is to be understood that these dependent claims may, alternatively, be made to depend in the alternative from any preceding or following claim, whether independent or dependent. Such new combinations are to be understood as forming a part of the present specification.
While the present invention has been described above by reference to various embodiments, it should be understood that many changes and modifications can be made to the described embodiments. It is therefore intended that the foregoing description be regarded as illustrative rather than limiting, and that it be understood that all equivalents and/or combinations of embodiments are intended to be included in this description.
| Number | Date | Country | Kind |
|---|---|---|---|
| 23204348.9 | Oct 2023 | EP | regional |
