Method for modeling acoustics in an underware

Abstract

A method for modeling acoustics in a delimited area, i.e., a computer-implemented method of modeling acoustics of an object as a numerical simulation model to determine field variables within a delimited control volume, the method comprising: (a) providing a coordinate system mapping positions in the control volume (b) providing a geometric specification of a sound source in the control volume based on the coordinate system. The objective of the invention is to improve the efficiency of numerical modeling of a transient (vibration) acoustic wave propagation phenomenon in a boundary-free domain. Additional steps are proposed: (c) constructing a convex shape around a sound source using a fast convex hull algorithm (d) stretching an infinite element starting from the convex shape (e) constructing a field shape function (g) using auxiliary coordinates by mapping radial coordinates of a coordinate system to auxiliary coordinates extending from the convex shape to the infinite element to transform the coordinate system into a field variable coordinate system (f1) Determining field variables by solving a numerical simulation model provides a computer system for simulating a system by applying a computer-implemented method.

Description

Technical Field

[0001] The present invention relates to a computer-implemented method for modeling acoustics in an unbounded domain to determine field variables within a control volume of the domain, the method comprising:

[0002] (a) Provide a coordinate system for mapping the control volume,

[0003] (b) providing a geometrical specification of an object that acts as a sound source and emits sound in the control volume, wherein the geometrical specification is based on the coordinate system.

[0004] Furthermore, the present invention also relates to a computer product, a computer-readable medium and a system.

[0005] The analysis of the sound characteristics of industrial machines is crucial in various applications (e.g. in the automotive, aerospace and marine industries). On the one hand, in many cases engineers need to pay attention to the sound quality and sound levels when designing their products to make them more attractive to potential customers. On the other hand, the increasing awareness of the potential health risks associated with exposure to high sound levels has led to stricter regulations on sound emissions in urban environments.

[0006] Additionally, in the industrial world, there is a growing interest in digital twins. In this context, high-fidelity 3D simulation models are used not only during the design phase of the product cycle, but also during the operational phase, where they are combined with measurement data through virtual sensing algorithms to gain a deeper understanding of the product's state. High-fidelity 3D simulation models are typically more complex and more expensive to develop, but provide a more accurate and realistic simulation of the modeled system compared to proxy models.

[0007] The intelligent virtual sensing algorithm is known from

[24] .

[0008] In the context of this disclosure, a digital twin is a virtual replica of a physical object, system, or process used to simulate and analyze its behavior and performance. Digital twins are typically created using data from sensors, cameras, and other sources that capture information about a physical object or system and combine it with a simulation model of the object. Digital twins can be used to test and optimize their performance under different conditions. They can better monitor and control physical systems, enabling predictive maintenance and performance improvements.

[0009] Digital twin technology is very useful in many scenarios, for example, in structural health monitoring applications. Due to the nature of the problem, a compact representation of a high-fidelity model is required, which can be used to run stable time-domain simulations. For example, digital twin technology is known from the following publications: "The Executable Digital Twin: merging the digital and the physics worlds" (2022) arXiv preprint arXiv:2210.17402 by D. Hartmann and H. Van der Auweraer, or "A Digital Twin for full-field monitoring in multi-channel control with applications to directfield acoustic testing" by AG De Miguel, O. Atak and MA Blanco.

[0010] In many (vibro)acoustic applications, a model of sound waves travelling to infinity is desired, meaning that the sound waves are travelling far away and are not reflected back. Background Art

[0011] The following publications provide some background information known to those skilled in the art:

[0012] [1] D. Givoli, “Computational absorbing boundaries”, Computational acoustics of noise propagation in fluids-finite and boundary element methods, Springer, 2008, pp. 145–166.

[0013] [2] A. Bayliss and E. Turkel, “Radiation boundary conditions for wave-like equations”, Communications on Pure and applied Mathematics, vol. 33 (1980), pp. 707–725.

[0014] [3] J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves”, Journal of computational physics, vol. 114 (1994), pp. 185–200.

[0015] [4] U. Basu and AK Chopra, “Perfectly matched layers for transient elastodynamics of unbounded domains”, International journal for numerical methods in engineering, vol. 59 (2004), pp. 1039–1074.

[0016] [5] J. Diaz and P. Joly, “A time domain analysis of pml models in acoustics”, Computer methods in applied mechanics and engineering, vol. 195 (2006), pp. 3820–3853.

[0017] [6] Z. Chen and X. Wu, “Long-time stability and convergence of the uniaxial perfectly matched layer method for time-domain acoustic scattering problems”, SIAM Journal on Numerical Analysis, vol. 50 (2012), pp. 2632–2655.

[0018] [7] A. Modave, A. Kameni, J. Lambrechts, E. Delhez, L. Pichon, C. Geuzaine, “An optimal pml for scattering problems in the time domain”, The European Physical Journal—Applied Physics, vol. 64 (2013).

[0019] [8] A. Modave, J. Lambrechts, and C. Geuzaine, “Perfectly matched layers for convex truncated domains with discontinuous galerkin time domain simulations,” Computers & Mathematics with Applications, vol. 73 (2017), pp. 684–700.

[0020] [9] S. Francois, H. Goh, L. F. Kallivokas, “Non-convolutional second-order complex-frequency-shifted perfectly matched layers for transient elastic wave propagation”, Computer Methods in Applied Mechanics and Engineering, 377 (2021), 113704.

[0021]

[10] B. Kaltenbacher, M. Kaltenbacher, and I. Sim, “A modified and stable version of a perfectly matched layer technique for the 3-d second order wave equation in time domain with an application to aeroacoustics,” Journal of Computational Physics, vol. 235 (2013), pp. 407–422.

[0022]

[11] P. Bettess and O. Zienkiewicz, “Diffraction and refraction of surface waves using finite and infinite elements”, International Journal of Numerical Methods in Engineering, vol. 11 (1977), pp. 1271–1290.

[0023]

[12] K. Gerdes, “The conjugated vs. the unconjugated infinite element method for the Helmetholtz equation in exterior domains”, Computer Methods in Applied Mechanics and Engineering, vol. 152 (1998), pp. 125–145.

[0024]

[13] R. Astley, G. Macaulay, J.-P. Coyette, and L. Cremers, “Three-dimensional wave-envelope elements of variable order for acoustic radiation and scattering. Part i. formulation in the frequency domain”, The Journal of the Acoustical Society of America, vol. 103 (1998), pp. 49–63.

[0025]

[14] L. Cremers, K. Fyfe, and J. Coyette, “A variable order infinite acoustic wave envelope element”, Journal of Sound and Vibration, vol. 171 (1994), pp. 483–508.

[0026]

[15] O. Zienkiewicz, K. Bando, P. Bettess, C. Emson, and T. Chiam, “Mapped infinite elements for exterior wave problems”, International Journal of Numerical Methods in Engineering, vol. 21 (1985), pp. 1229–1251.

[0027]

[16] W. Eversman, “Mapped infinite wave envelope elements for acoustic radiation in a uniformly moving medium”, Journal of Sound and Vibration, vol. 224 (1999), pp. 665–687.

[0028]

[17] DS Burnett, “A three-dimensional acoustic infinite element based on a prolate spheroidal multipole expansion”, Journal of the Acoustical Society of America, vol. 96 (1994), pp. 2798–2816.

[0029]

[18] DS Burnett and RL Holford, “Prolate and oblate spheroidal acoustic infinite elements”, Computer Methods in Applied Mechanics and Engineering, vol. 158 (1998), pp. 117–141.

[0030]

[19] DS Burnett and RL Holford, “An ellipsoidal acoustic infinite element”, Computer Methods in Applied Mechanics and Engineering, vol. 164 (1998), pp. 49–76.

[0031]

[20] R. Astley, “Transient wave envelope elements for wave problems”, Journal of Sound and Vibration, vol. 192 (1996), pp. 245–261.

[0032]

[21] R. Astley and J. Hamilton, “Infinite elements for transient flowacoustics,” in 7th AIAA / CEAS Aeroacoustics Conference and Exhibit, 2001, p. 2273.

[0033]

[22] R. Astley, J.-P. Coyette, and L. Cremers, “Three-dimensional wave-envelope elements of variable order for acoustic radiation and scattering. Part ii. formulation in the time domain”, Journal of the Acoustical Society of America, vol. 103 (1998), pp. 64–72.

[0034]

[23] S. van Ophem, O. Atak, E. Deckers, and W. Desmet, “Stable model order reduction for time-domain exterior vibro-acoustic finite elements simulations”, Computer Methods in Applied Mechanics and Engineering, vol. 325 (2017), pp. 240–264.

[0035]

[24] "Modelorder reduction and smart virtual sensing for unbounded vibroacoustics using high order fem and infinite elements" by D.Bizzarri, O.Atak, S.van Ophem, H.Bériot, T.Tamarozzi, P.Jiranek, L.Scurria, A.Garcia de Miguel, M.Alvarez Blanco, K.Janssens, W.Desmet, "ISMA (International Conference on Noiseand Vibration Engineering) / USD Conference Proceedings (Proceedings of ISMA / USD)》2022(2022).

[0036]

[25] A. van de Walle, “The power of model order reduction invibroacoustics: And its applications in model-based sensing,” PhD thesis, University of Leuven, 2018.

[0037]

[26] "Novel reduction techniques for exterior vibro-acoustic models and their use in model-based sensing and identification" by S. van Ophem (2019).

[0038]

[27] S. van Ophem, B. Forrier, D. Bizzarri, O. Atak, A. de Miguel, M. Elkafafy, M. Dal Borgo, T. Tamarozzi, M. Alvarez Blanco, K. Janssens, et al., “Physics-based virtual acoustic sensing for enhanced direct field acoustic excitation testing,” ISMA / USD Proceedings, 2022 (2022).

[0039]

[28] R. Astley and J. Hamilton, “The stability of infinite element schemes for transient wave problems”, Computer Methods in Applied Mechanics and Engineering, vol. 195 (2006), pp. 3553–3571.

[0040]

[29] R. Holford, “A multipole expansion for the acoustic field exterior to a prolate or oblate spheroid,” Journal of the Acoustical Society of America, submitted (1998).

[0041]

[30] R. Astley, “Mapped spheroidal wave-envelope elements for unbounded wave problems”, International Journal of Numerical Methods in Engineering, vol. 41 (1998), pp. 1235–1254.

[0042]

[31] R. Astley, “Transient spheroidal elements for unbounded wave problems”, Computer Methods in Applied Mechanics and Engineering, vol. 164 (1998), pp. 3–15.

[0043]

[32] R. Astley, “Infinite elements for wave problems: a review of current formulations and an assessment of accuracy”, International Journal of Numerical Methods in Engineering, vol. 49 (2000), pp. 951–976.

[0044]

[33] R. Astley and J.-P. Coyette, “The performance of spheroidal infinite elements”, International Journal of Numerical Methods in Engineering, vol. 52 (2001), pp. 1379–1396.

[0045]

[34] L.-X. Li, S. Kunimatsu, J.-S. Sun, and H. Sakamoto, “A new conjugated mapped infinite element,” Journal of Computational Acoustics, vol. 12 (2004), pp. 543–570.

[0046]

[35] D. Dreyer and O. von Estorff, “Improved conditioning of infinite elements for exterior acoustics”, International Journal of Numerical Methods in Engineering, 58 (2003), 933–953.

[0047]

[36] CB Barber, DP Dubkin, and H. Huhdanpaa, “The quickhull algorithm for convex hulls,” ACM Transactions on Mathematical Software (TOMS), vol. 22 (1996), pp. 469–483.

[0048]

[37] Siemens Industrial Software Simcenter Nastran User Guide (Version 2022.1), https: / / docs.sw.siemens.com / en-US / doc / 289054037 /

[0049] PL20201105151514625.xid1853788 / xid1853789, 2022.

[0050]

[38] R. J. Astley, “Infinite elements,” in Computational Acoustics of Noise Propagation in Fluids-Finite and Boundary Element Methods, Springer, 2008, pp. 197–230.

[0051]

[39] R. Astley, G. Macaulay, and J. Coyette, “Mapped wave envelope elements for acoustical radiation and scattering”, Journal of Sound and Vibration, vol. 170 (1994), pp. 97–118.

[0052]

[40] CH Wilcox, “A generalization of theorems of Rellich and Atkinson”, Proceedings of the American Mathematical Society, vol. 7 (1956), pp. 271–276.

[0053]

[41] "Zur modellierung mitboundary elementen und finiten elementen bei schallabstrahlungsberechnungen.1" by O.Zaleski, W.-C.von Karstedt, O.von Estorff, in Deutschsprachige Anwender-Konferenz SYSNOISE, 1999.

[0054]

[42] D. Dreyer, S. Petersen, and O. von Estorff, “Effectiveness and robustness of improved infinite elements for exterior acoustics”, Computer Methods in Applied Mechanics and Engineering, vol. 195 (2006), pp. 3591–3607.

[0055]

[43] R. Astley and J. Hamilton, “Numerical studies of conjugated infinite elements for acoustical radiation”, Journal of Computational Acoustics, vol. 8 (2000), pp. 1–24.

[0056]

[44] D. Dreyer, Efficient infinite elements for exterior acoustics, Shaker Press, 2004.

[0057]

[45] K. Gerdes, “A review of infinite element methods for exteriorhelmholtz problems”, Journal of Computational Acoustics, vol. 8 (2000), pp. 43–62.

[0058]

[46] J. Cipolla, “Acoustic infinite elements with improved robustness,” in ISMA Proceedings, 2002, pp. 2181–2187.

[0059]

[47] H. Bériot and A. Modave, “An automatic perfectly matched layer for acoustic finite element simulations in convex domains of general shape”, International Journal of Numerical Methods in Engineering, vol. 122 (2021), pp. 1239–1261.

[0060]

[48] ​​H. Bériot, A. Prinn, and G. Gabard, “Efficient implementation of high-order finite elements for helmholtz problems”, International Journal of Numerical Methods in Engineering, vol. 106 (2016), pp. 213–240.

[0061]

[49] H. Bériot and G. Gabard, “Anisotropic adaptivity of the p-fem for time-harmonic acoustic wave propagation”, Journal of Computational Physics, vol. 378 (2019), pp. 234–256.

[0062]

[50] J. Biermann, O. von Estorff, S. Petersen, and C. Wenterodt, “Higher order finite and infinite elements for the solution of helmholtz problems”, Computer Methods in Applied Mechanics and Engineering, vol. 198 (2009), pp. 1171–1188.

[0063]

[51] C. Geuzaine and J.-F. Remacle, “Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities”, International Journal of Numerical Methods in Engineering, vol. 79 (2009), pp. 1309–1331.

[0064]

[52] Math 3795: Introduction to computational mathematics by D. Leykekhman, 2008.

[0065] Several options for numerically modeling transient (vibro)acoustic wave transport phenomena in unbounded domains are disclosed in the literature; the most popular techniques include:

[0066] - modeling the external domain using infinite elements, and

[0067] -Perfectly matched layer [PML].

[0068] PML is a technique used in computational acoustics to model the behavior of acoustic waves in unbounded regions. PML was originally developed for use in electromagnetic simulations. The PML technique involves building a layer around the computational domain that gradually absorbs the waves as they travel away from the domain. The layer is designed to simulate the properties of an infinite medium, effectively eliminating reflections and preventing the waves from leaking into the surrounding space.

[0069] At the continuous level, no reflection occurs at the interface between the FEM and the external domain for any angle of incidence [3]. There are various PML formulations for time domain simulations [10; 7; 5].

[0070] However, there are several problems associated with such transient PML concepts:

[0071] a) To make the system matrix frequency independent, new auxiliary variables are introduced in the PML layer. This means that for transient calculations, the number of DOFs (degrees of freedom) is quite large compared to the corresponding frequency domain simulations. This will result in a lot of time and memory required to build the ROM.

[0072] b) These ideas are often quite complex and difficult to implement in simulation software.

[0073] c) Stability analysis of many of the proposals presented in the literature is not easy to perform and each proposal needs to be examined individually [5]. For the proposals presented in

[10] , some terms need to be suppressed to stabilize the transient simulation, which introduces new approximation errors.

[0074] d) It is known from the literature

[26] that in the case of model order reduction (MOR), in order to guarantee that stability properties are preserved after the projection step, it is necessary to use a split basis method, which results in a larger reduced order model. In numerical experiments analyzing a plate radiating in a free field, the authors compared the efficiency of PML and infinite elements for time-domain MOR: the authors showed that infinite elements result in a smaller full-order model (making the reduction procedure faster) and that the number of DOFs in the final ROM is 3 times less than that obtained from the PML configuration. The increased DOFs lead to longer computational times for the PML concept.

[0075] Several Model Order Reduction (MOR) techniques are known in the field of engineering and applied mathematics. Here are some of the most common ones:

[0076] 1. Proper Orthogonal Decomposition (POD): This technique reduces the dimensionality of a system by using statistical methods to identify the most important patterns of variation in the system. It is often used in fluid dynamics, structural mechanics, and other fields.

[0077] 2. Balanced Truncation: This technique reduces the dimensionality of a system by balancing the energy of the system between its input and output. It is often used in control systems and linear systems.

[0078] 3. Krylov subspace method: This technique reduces the dimensionality of the system by approximating the response of the system using a subspace of the Krylov sequence. It is often used in linear systems and control systems.

[0079] 4. Model Order 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. It is often used in control systems, circuit simulation, and other fields.

[0080] 5. Reduced Basis Method: This technique reduces the dimensionality of a system by building a basis for the solution space of the system using snapshot solutions. It is often used in infinite element analysis, structural mechanics, and other fields.

[0081] These techniques are used to reduce the complexity of mathematical models and simulations, making them easier to solve and more computationally efficient. They are widely used in engineering, physics, and other fields to analyze and optimize complex systems and processes.

[0082] The infinite element method (IFEM) is a popular alternative to PML and is particularly competitive for time-domain models.

[0083] The Astley-Leis infinite element formulation

[20] produces a frequency-independent system matrix and can be easily used for transient simulations. There are two types of Astley-Leis formulations:

[0084] 1. The standard Astley-Leis conception

[39] is based on a truncated multipole expansion written in spherical coordinates (the so-called Atkinson Wilcox expansion

[40] );

[0085] 2. On the other hand, the spherical conception

[28] is based on a truncated multipole expansion written in spherical coordinates.

[0086] So the main difference between spherical and spherical coordinates is that spherical coordinates describe a point in terms of its distance from a fixed point and two angles, whereas spherical 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.

[0087] Both formulations rely on a singular mapping between infinite physical elements and the corresponding finite parent elements. For each infinite element, the distance to the so-called virtual source is mapped to the parent element radial coordinate; to form the shape function, a common polynomial written in terms of the parent coordinates is used. In order for the standard Astley-Leis formulation to accurately capture the correct solution, the virtual source should coincide with the center of the smallest sphere (also called the radiation center) that encloses all the sound sources in the problem being analyzed. When this constraint is not met, the accuracy of the solution cannot be guaranteed.

[0088] In this case, the parent element is a reference element used as a template for building the finite element model of the physical domain. It may have a simple geometric shape such as a triangle, quadrilateral, tetrahedron, or hexahedron.

[0089] The parent element is defined by a set of nodes or knots that are used to interpolate the solution within the element. The knots can 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 the continuity and differentiability conditions. The parent element can be transformed and mapped onto the physical domain using a set of geometric transformations such as affine transformations or isoparametric mappings.

[0090] In this case, the shape function is a mathematical function used to interpolate the values ​​of the field variables within each finite element.

[0091] The use of the standard Astley-Leis formulation allows the use of any convex FEM domain for frequency domain calculations. However, for transient simulations relying on such formulations, the infinite elements should be stretched orthogonally to the outer FEM boundary envelope to guarantee the stability of the method [4]. This, together with the above-mentioned criteria for the positioning of virtual sources, forces the use of a spherical cutoff of the FEM domain. For the analysis of elongated scatterers in the time domain, this means that when using the standard Astley-Leis infinite elements, an unnecessarily large mesh must be used, resulting in a prohibitive number of DOFs, such as Fig.18 shown.

[0092] The spherical Astley-Leis concept can circumvent this problem by allowing a smaller FEM domain to be fitted around the scatterer.

[0093] In [32, 33] it was shown that for such formulations the radial order of the shape function often needs to be higher than for the standard Astley-Leis formulation, offsetting the benefits of using a smaller FEM domain.

[0094] It thus appears that the problem of automatically stretching infinite elements has not been solved, which seems to be a problem that has not been addressed in previous literature.

[0095] A particularly challenging situation is when only a FEM mesh is available and the geometry was not modeled with a CAD tool. In this case, the geometric data required to extrude the infinite elements needs to be computed directly from the mesh of the FEM domain. A version of PML called Automatic Matching Layer (AML) was reported in

[47] for which the user has to provide only the scatterer geometry (CAD model or FEM mesh); the acoustic interior convex domain is then automatically generated using the so-called Quickhull algorithm (known from

[36] ) and the AML is automatically extruded during runtime. However, AML is only available for frequency domain simulations.

[0096] In this context, the fast convex hull algorithm is a computational geometry algorithm for computing the convex hull of a set of points in multidimensional 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, it is the outermost boundary of the set of points. The fast convex hull algorithm is a divide and conquer algorithm that recursively splits the set of points into smaller subsets and finds the convex hull of each subset.

[0097] The divide and conquer algorithm is a well-known problem solving strategy that involves breaking down a problem into smaller, more manageable subproblems and solving each subproblem individually. The solutions to the subproblems are then combined to solve the original problem. Here, this strategy is applied to a geometric problem as known from

[36] .

[0098] Therefore, the challenge of numerically modeling transient (vibro)acoustic wave propagation phenomena in unbounded domains remains open in the prior art. Summary of the invention

[0099] It is an object of the present invention to improve the efficiency of numerical modeling of transient (vibro)acoustic wave propagation phenomena in unbounded domains.

[0100] In order to solve this problem, the present invention presents a method according to the present invention. A computer-implemented method of the type mentioned earlier is proposed, the method comprising the additional steps of:

[0101] (c) constructing a convex shape around the sound source using a fast convex hull algorithm,

[0102] (d) Stretching the infinite element starting from the convex shape,

[0103] (e) transforming the coordinate system into a field variable coordinate system by mapping the radial coordinates of the coordinates to auxiliary coordinates extending from the convex shape to the infinite element,

[0104] (f1) Use auxiliary coordinates to construct shape functions,

[0105] (g) Determine the field variables by solving the numerical simulation model.

[0106] According to the present invention, the numerical simulation model involves numerical modeling of transient (vibration) acoustic wave propagation phenomena in an unbounded domain. Therefore, the model can be understood and referred to as a model for "numerical modeling of transient (vibration) acoustic wave propagation phenomena in an unbounded domain".

[0107] Preferably, the field variable is a pressure field. Most preferably, the field variable can be a pressure field or any derivative of a parametric pressure field.

[0108] The basic idea of ​​the present invention is an infinite element concept, which can be called flexible infinite elements, which can be used in conjunction with the finite element method (FEM) to effectively model the time domain propagation of acoustic waves in unbounded media. In particular, the new concept addresses the challenge of making such transient unbounded vibroacoustic simulations very fast and compact so that they can be used in digital twin scenarios, including applications where the simulations need to run in near real time.

[0109] The governing equations describing the physical behavior of acoustic waves can be written in the frequency domain or in the time domain. In frequency domain simulation, the input signal is represented as a sum of harmonics. By decomposing the signal into its individual frequency components, the steady-state pressure response at each frequency can be calculated separately. On the other hand, time domain simulation (also called transient simulation) models acoustic waves as pressure variations over time. In digital twin applications, the time domain description is very important.

[0110] The problem of efficiently and accurately simulating the propagation of acoustic waves in unbounded domains has been extensively studied over the past few decades. Non-reflective conditions can be modeled in various ways using the finite element method (FEM) or alternatively other numerical methods such as the boundary element method (BEM). The time domain version of the BEM does not lend itself to the construction of reduced order models (ROMs). Since such compact representations are needed in digital twin scenarios, the focus of the present invention is on the FEM.

[0111] Since FEM is a domain discretization method, it is inherently infeasible to model infinite domains. The FEM mesh must be truncated and non-reflecting boundary conditions must be implemented on the external boundaries to prevent acoustic waves from traveling back into the FEM domain.

[0112] The present invention presents a novel flexible infinite element concept. It is a unique approach that allows the use of generally convex envelopes for infinite elements while being time stable; its geometric flexibility makes it particularly attractive when modeling sound scattered by slender structures. It is also compatible with the MOR technique, since it has a frequency-independent system matrix. Therefore, it also allows a compact representation.

[0113] Unlike the existing concept, which constructs infinite element shape functions based on a standard mapping between the distance to a virtual source and radial parent coordinates, a new auxiliary mapping is introduced between radial coordinates and new auxiliary coordinates. The usual parent coordinates are used for the definition of geometric interpolation, while the polynomial of the shape function within the infinite element is expressed in terms of the new auxiliary parent coordinates, which consist of the standard parent coordinates in the tangential direction and the new auxiliary coordinates in the radial direction, as shown in the figure. In this concept, the shape function closely mimics the truncated spherical multipole expansion in any geometric configuration.

[0114] Therefore, the numerical scheme is able to provide accurate results regardless of the positioning of the mapping nodes in the infinite domain. The new concept is therefore geometrically flexible, allowing infinite elements to be attached to generally convex envelopes without compromising the accuracy of transient simulations.

[0115] In order to ensure the time stability of the flexible infinite element concept, an accuracy-maintaining stabilization procedure was successfully developed.

[0116] According to the present invention, a fast convex hull algorithm is used to generate a convex acoustic domain around any given scatterer, while the mapping nodes within the infinite element are extruded at runtime with normal directions orthogonal to the outer FEM boundaries calculated using the same rules as in

[12] . Due to the adopted strategy, user intervention is maintained to minimize the construction of the acoustic and exterior domains, even when no CAD of the scatterer geometry is available.

[0117] Additionally, since the proposed method produces a frequency-independent system matrix, it can be easily combined with MOR methods such as Krylov-based methods. This aspect and the geometric flexibility of the new concept open new digital twin applications in the field of unbounded (vibro)acoustics.

[0118] A preferred embodiment of the present invention provides that, in step (e), the auxiliary coordinates for positions outside the convex shape are defined as:

[0119]

[0120] in:

[0121] s, t, is the 3D coordinate of the position outside the parameterized convex shape (s, t are tangential coordinates, is the radial coordinate);

[0122] Identify the point on the convex boundary Γ that is closest to the external position, is the distance between the radiation center and the specific point on Γ;

[0123] defining a distance between a center of radiation and said location;

[0124] are auxiliary coordinates.

[0125] This transformation enables the acoustics to be modeled quickly and accurately and solved to determine the field variables in a reasonable time.

[0126] Further preferred embodiments provide for displaying to the user said field variables and / or key figures calculated from said field variables or images illustrating said field variables.By default, the display informs the user about the acoustic behavior of the object using the display as a human-machine interface.

[0127] A further preferred embodiment comprises the step of providing said field variables to an iterative object design process of said object for designing said object with improved acoustic emissions as a design goal.

[0128] A further preferred embodiment comprises the step of generating said object according to the design as determined by said iterative object design process.The generating can be done by any manufacturing or production method.

[0129] The iterative object design can include the following steps:

[0130] selecting at least one design parameter of the object (OBJ) design to be improved;

[0131] providing an iteration termination criterion as a threshold value of the field variable or a key number calculated based on the field variable;

[0132] defining at least two variants (VRT) of the object (OBJ), the two variants differing in selected design parameters;

[0133] Acoustic modeling of said variant (VRT) according to the invention;

[0134] selecting the best variant from the defined variants based on the modeling results;

[0135] Repeat the steps of definition, modeling, and selection until the termination criteria are met;

[0136] Selecting design parameters of a variant that satisfies a termination criterion;

[0137] Optionally, at least one different design parameter is selected and the same iterations as defined above for the newly selected at least one different design parameter are processed.

[0138] The process is robust and efficient.Other iterative object design strategies may also be applicable.

[0139] Furthermore, the invention also relates to a computer product arranged and configured to perform the steps of the computer-implemented method according to the invention.

[0140] Furthermore, the invention also relates to a computer-readable medium encoded with executable instructions which, when executed, cause a computer product as defined above to implement the method according to the invention.

[0141] Furthermore, the invention relates to a system for determining acoustic parameters of an emission of an object, the system comprising at least one processor which is prepared to execute the method according to the invention by uploading a computer executable code.

[0142] 1. Introduction

[0143] The finite element method (FEM) is one of the most popular techniques for simulating acoustic wave propagation, as it allows complex problems to be solved efficiently by discretizing the physical domain of interest. Wave field simulations in unbounded domains present challenges, as computers can only handle a finite number of elements. The simplest, but rather inaccurate approach consists in applying impedance boundary conditions at the artificial boundaries of the finite element domain [1]. Higher-order absorbing boundary conditions [2] and perfectly matched layers (PMLs) [3] represent more accurate alternatives. In particular, PML has been shown to be very effective for solving frequency-domain wave problems. PML has also been extensively studied for transient applications (see [4] to

[10] and references therein). However, for time-domain simulations, PML requires the use of auxiliary variables; the resulting increase in the number of degrees of freedom (Dof) makes PML less attractive for simulating transient phenomena.

[0144] The infinite element method proposed by Bettess and Zienkiewicz

[11] avoids the problem of truncating the computational domain by employing special elements that extend from the boundary of the FE domain to infinity. Different formulations have been proposed in the literature over the years. In particular, two groups of formulations have been developed that differ in the choice of the test function in the infinite domain: non-conjugate formulations, where the test basis function and the shape function are the same, and conjugate formulations, where the test basis function is the complex conjugate of the shape function. A detailed comparison of these different strategies is presented in

[12] , where the advantages and disadvantages of both are analyzed in detail. The Astley-Leis formulation

[13] differs from the classical conjugate formulation only in the weighting factor by which the test function is multiplied, which is designed to cancel the unbounded terms in the integral that needs to be calculated to obtain the system matrix.

[0145] There are also alternative geometric discretizations of infinite elements. So-called mapping proposals rely on a mapping from physical elements to mother elements

[13] –

[16] : their main advantage is geometric flexibility, since they can be attached to arbitrary convex FE boundaries. In contrast, for Burnett elements

[17] –

[19] , the geometry of the envelope is restricted to a spheroid and no mapping is employed.

[0146] Astley-Leis elements have several attractive properties. First, due to the weighting factors used in the test function, the system matrix integral can be evaluated using the standard quadrature rule. In addition, they produce frequency-independent system matrices and preserve the sparsity of the FEM system, which makes them readily available for time-domain simulations

[20] to

[22] and suitable for model order reduction (MOR)

[23] . This is particularly interesting because in recent years MOR for transient simulations has become one of the key techniques in many applications

[24] to

[27] . The goal of the MOR technique is to construct a reduced-order model (ROM) that captures the most relevant dynamics of the corresponding initial full-order model (FOM) using only a small fraction of the Dof. Numerical tests were presented in

[26] , where the performance of the transient concept from Kaltenbacher et al.

[10] was compared with the Astley-Leis infinite element concept. Specifically, the results show that the use of PML results in a higher Dof due to the presence of auxiliary variables. Moreover, from the MOR point of view, it is necessary to use a split basis method, which results in a three-fold increase in the number of Dofs of the final ROM, in order to ensure that the time-domain stability of the PML concept after MOR is preserved. Since the performance comparison between the PML and infinite element methods is limited to a single numerical test, it is difficult to draw general conclusions. However, the cited studies seem to give strong indications that using infinite elements for unbounded transient problems may be the best choice, especially when using the MOR technique.

[0147] Astley and Hamilton

[28] showed that the Astley-Leis formulation is time-stable as long as the infinite edges are aligned with the normals of the interface between the FEM and the external domain. However, the performance of conventional mapped infinite elements varies with the location of the virtual sources

[14] ; the stability requirements mentioned do not allow the location of the virtual sources to be chosen freely, and this often represents an accuracy problem in transient analysis when analyzing arbitrary convex interior domains. In practice, for most scenarios this means that the mapped Astley-Leis infinite elements should be attached only to the spherical envelope for transient simulations. Based on the Holford expansion

[29] , alternative mapped sphere infinite element formulations have been proposed in both the frequency and time domains

[30] ,

[31] . However, for such formulations as well as non-mapped sphere formulations

[17] to

[19] , the benefits of a smaller interior domain may be cancelled by the necessary increase in the radial order of the infinite element shape function, as observed in

[31] ,

[32] , and

[33] .

[0148] Li et al.

[34] proposed an alternative conjugate mapping formulation in which the Lagrangian-type radial shape function is modified to exactly satisfy the truncated multipole expansion along the infinite edge, regardless of the location of the virtual source. However, they only solved the problem in the frequency domain, where the standard Astley-Leis formulation already provides satisfactory results; moreover, Li's elements are poorly conditioned for high radial orders due to the restriction that the radial shape function is derived from Lagrangian polynomials

[35] .

[0149] In this work, we propose a novel infinite element formulation for any type of radial shape function that exactly satisfies the truncated multipole expansion regardless of the orientation of the infinite edges. Thus, we employ new flexible infinite elements to achieve a time-stable and accurate formulation of arbitrary convex envelopes. The integration of the flexible infinite elements into high-level software is discussed. Strategies are proposed to automatically stretch the infinite elements when only the mesh of the interior domain is available. In addition, we show how to exploit the fast convex hull algorithm

[36] implemented in Simcenter 3D

[37] , a simulation package developed by Siemens Industry Software, to automatically generate the envelopes attached to the infinite elements around the scatterer. Finally, the combination of the proposed formulation with existing MOR strategies is discussed. The paper is organized as follows. In Section 2, the Helmholtz problem and the resulting system of time-domain equations are introduced, along with a description of the standard Astley-Leis formulation. In Section 3, the novel flexible infinite element is proposed and its stability properties are discussed. In Section 5, we demonstrate the effectiveness of the proposed approach through two academic examples, while in Section 6 we consider industrial complexity problems. Finally, we present the conclusions in Section 7.

[0150] 2. Finite element and infinite element in acoustics

[0151] 2.1. Problem Statement

[0152] 2.1.1. Governing equations

[0153] We define the time dependence of the field as e iωT , where ω is the angular frequency and τ is time.

[0154] The behavior of the pressure field p in the unbounded domain Ω is described by the Helmholtz equation in the frequency domain:

[0155]

[0156] where k is the wave number. The domain Ω is subdivided into the inner domain Ω i and external domain Ω e ; We use Γ to denote the interface between two domains (see Figure 1 ).

[0157] Figure 1 Schematic diagram showing the external acoustic problem. Ω i and Ω e are the inner domain and the outer domain, respectively. Γ is the interface between the two domains, and Γs is the envelope where the Sommerfeld condition is imposed.

[0158] In the domain Ω i The inner boundary of Γ N At , Neumann boundary conditions are imposed. Dirichlet and Robin boundary conditions can also be easily handled, but for ease of presentation they will not be considered here. The Sommerfeld condition states that no waves can be reflected at infinity; in mathematical terms, in the d-dimensional domain, this can be written as:

[0159] Since r→∞, (2)

[0160] Where r is the radius of the sphere / cylinder.

[0161] As shown in

[38] , approximating (2), we can rewrite:

[0162] In Γ S Up, (3)

[0163] in, is the gradient operator in Cartesian coordinates, and Γs is the envelope with outward unit normal n, placed at a radial distance R from the envelope Γ, where R→∞.

[0164] The corresponding boundary value problem is as follows:

[0165]

[0166] where ρ is the fluid density and v n is in Γ N The normal velocity specified above.

[0167] 2.1.2. Variational Conception and Discretization

[0168] Using the weighted residual method, we obtain the variational formulation:

[0169]

[0170] Among them, H 1 is the Sobolev space in which the solution is sought.

[0171] Equation (5) is obtained by i Use H 1The solution is obtained by conforming the finite element and attaching the infinite element to the envelope Γ. The field variable p and the test function q are given by and Approximately, where V h is the shape function space, and Q h is the test function space. In the Astley-Leis formulation, the Sommerfeld condition is automatically satisfied due to the weighting factors used in the test function and does not need to be considered in the variational formulation

[39] . In this case, the discrete variational formulation is:

[0172]

[0173] It can be rewritten as a system of linear equations:

[0174] [-ω 2 M+iωC+K]p(ω)=f(ω) (7)

[0175] Where p is the Dof vector, f is the force vector, and M, C, and K∈R nDOF×nDOF are the mass matrix, damping matrix and stiffness matrix respectively, and n DOF is the number of Dof. The exact formula for the system matrix can be found in

[13] .

[0176] Since, for the Astley-Leis formulation, the system matrix is ​​frequency independent, equation (7) can be easily transformed into the time domain:

[0177]

[0178] Among them, p τ and f τ are the inverse Fourier transforms of p and f respectively.

[0179] 2.1.3. Multipole expansion

[0180] The Atkinson-Wilcox theorem

[40] states that at any point outside the smallest sphere enclosing all sound sources, the solution to problem (4) can be written as a multipole expansion:

[0181]

[0182] Among them, G n is the directivity function, α and θ are spherical angles, and r is the radial distance, measured from the center of the radiation coinciding with the center of the smallest sphere. For ease of notation, in the rest of this article we will assume that the origin O of the absolute Cartesian coordinate system coincides with the center of the radiation.

[0183] Holford

[29] has extended the Atkinson-Wilcox theorem to include the spherical case. By placing the envelope Γ outside the minimal sphere (for the Atkinson-Wilcox expansion) or the spheroid (for the Holford expansion), we ensure that the outer domain Ω e The solution in consists only of outward propagating waves: this consideration is used to select appropriate trial solutions in the infinite element.

[0184] The so-called sphere conception has been criticized in the past

[32] ,

[33] , since in practical applications the lower number of Dof in the FEM region is often offset by the need for a higher radial order in the outer domain. On the other hand, numerical experiments have been performed in the frequency domain using mapped infinite elements based on the expansion (9), where the conditions of the Atkinson-Wilcox theorem are explicitly violated

[41] ,

[35] ,

[42] : namely, non-spherical envelopes are used, which partially lie inside the smallest enclosed sphere. Although the convergence of this approach has not been formally proven

[43] , numerical experiments with sufficiently smooth envelopes give promising results (see

[44] for more details). For all test cases encountered during our study, our observations are consistent with those in the cited literature. These are further expanded in Section 6.

[0185] 2.2. Overview of Astley-Leis infinite element

[0186] In this subsection, we describe the main features of the Astley-Leis formulation, and we analyze the conditions under which it is stable in the time domain.

[0187] 2.2.1. Geometric interpolation

[0188] The Astley-Leis conception relies on a mapping from infinite physical elements to finite parent elements. Specifically, in two dimensions (see Figure 2 ), the physical element is mapped to the secondary mother element In three dimensions, the parent element can be a triangular prism or a cube, depending on the shape of the corresponding physical element, and the parent element coordinates are named s, t, and symbol is introduced for all variables related to the concept of virtual source (see how to define virtual source later in this section). Specifically, after introducing radial mappings (16) and (24), the radial coordinates of the mother element are The connection to the distance from the virtual source in physical space will become clear.

[0189] Figure 2 A schematic diagram of the mapping infinite element is shown. In our convention, the radiation center coincides with the origin O of the Cartesian coordinate system.

[0190] exist Figure 2 In , we schematically represent the infinite elements mapped in two dimensions, showing both the physical elements and the parent elements. For the i-th node Γi i e , whose position is defined during the meshing process of the inner domain, on the envelope Γ (at the interface between the outer and inner domains), and in the eth infinite element, the user defines the mapping node in the outer domain The Cartesian coordinates of the position are:

[0191]

[0192] Similarly, we will map the node Γi i e The Cartesian coordinates on Γ are written as:

[0193]

[0194] The i-th infinite edge is from node Γi i e Issued, passing through node and extending to infinity (see Figure 2 ). For all nodes on Γ, infinite edges are defined.

[0195] We use Represents a mapping node With Γi i e Distance between:

[0196]

[0197] Where ‖·‖ is the Euclidean norm. Virtual source The Cartesian coordinates are:

[0198]

[0199] It is located at a distance of Γ along the direction of the i-th infinite edge toward the inner domain i e distance (See Figure 2 ):

[0200]

[0201] Therefore, the virtual source With external domain The distance between the mapping nodes in for:

[0202]

[0203] The envelope Σ is the union of all points in the physical space, and its parent radial coordinate is For any point on the infinite edge, the distance from the virtual source is Mapping to the mother element radial coordinates

[0204]

[0205] Note that with this mapping definition, we have

[0206]

[0207] The interpolated distance between the mapped nodes on Γ and their corresponding virtual sources Defined as:

[0208]

[0209] where n is the number of mapping nodes on the envelope Γ in the e-th infinite element, and is the standard tangent geometry mapping function (in our case, we use Lagrange polynomials). Similarly, the interpolated distances between the mapped nodes on the envelope Σ and their corresponding virtual sources are Defined as:

[0210]

[0211] By substituting equation (15) into (19), we obtain:

[0212]

[0213] Then by substituting equation (18) into equation (20), we obtain:

[0214]

[0215] The interpolation distance between a mapping node in element e and its corresponding virtual source is defined as:

[0216]

[0217] in, and is the radial geometric mapping function, defined as:

[0218]

[0219] Substituting equation (23) and equation (21) into equation (22), we finally get:

[0220]

[0221] In other words, for any point inside the infinite element, the physical coordinates can be obtained With radial coordinates The mapping between them is , and Equation (16) is only valid along the i-th infinite edge. The interpolation distance of the virtual source mapped to the e-th infinite element explains The relationship between the concept of virtual source and thus the use of the symbol It is reasonable to use it as the radial coordinate of the parent element. Note that by this mapping, we obtain:

[0222]

[0223] This mapping is consistent with the concept of having a parent radial coordinate The points are located on the envelope Γ in physical space (see Figure 2 ),have The points are located on Σ and have The point is at infinity in the radial direction from the interior domain. In addition, along the i-th infinite edge, we have:

[0224]

[0225] The Cartesian coordinate vector inside the e-th infinite element is obtained by interpolating the Cartesian coordinates of the nodes:

[0226]

[0227] in, and is the Cartesian coordinate vector of the mapping node on the i-th infinite edge.

[0228] We emphasize that the notation adopted here is slightly different from that usually found in the literature, because we want to distinguish between r defined in Section 2.1.3 (i.e., the distance from the center of radiation) and The interpolation distance between the mapping node and its corresponding virtual source is defined in Eq. (22). Specifically, since we agree that the origin O of the Cartesian coordinate system and the radiation center r e Coincidence means that the distance between the point in the e-th infinite element and the radiation center can be written as:

[0229]

[0230] In addition, we will It is defined as the distance between any point on the i-th infinite edge and the radiation center O.

[0231] Finally, the so-called phase distance Defined as:

[0232]

[0233] When the infinite edge is stretched orthogonally to Γ, the phase distance The physical coordinates are The distance between the point and the envelope Γ, such as Figure 2 shown.

[0234] 2.2.2. Shape functions and test functions

[0235] By analogy with the Atkinson-Wilcox expansion (9), the field shape function is written as:

[0236]

[0237] in, is the usual vector of polynomial shape functions used in the FE domain. As an example, for ease of representation, let us consider a first-order Lagrangian polynomial in two dimensions:

[0238]

[0239] Inserting the radial mapping (24) into equation (31), we can verify that the following relationship holds:

[0240]

[0241] Expanding the above argument to a higher-order polynomial in three dimensions, we obtain:

[0242]

[0243] where m is the radial order chosen in the infinite domain. Specifically, along the i-th infinite edge, we have

[0244]

[0245] In other words, for the i-th infinite edge, the pressure field generated by the m-th order virtual multipole source is placed at This is where the idea of ​​a virtual source comes from.

[0246] The interesting case is when all virtual sources of the e-th infinite element coincide with the radiation center O. In this case, along the i-th infinite edge, we obtain:

[0247]

[0248] Because when When, along the i-th infinite edge In other words, when all virtual sources coincide with the radiation center, the approximate field variable ph It is possible to simulate the multipole expansion (9), truncated along the infinite edge with radial order m, ensuring good accuracy of the approximation. In all other cases, relation (35) is not necessarily true. This explains why for the traditional mapping infinite element conception, the accuracy depends on the location of the virtual source. Moreover, inside the e-th infinite element, in general, we have:

[0249]

[0250] This comes from the fact that in equation (22), the same Cartesian coordinate vector x is used to obtain the Cartesian coordinate vector x in equation (27). e The same interpolation function is obtained However, the radial distance r e is not a linear combination of Cartesian coordinates and therefore cannot be found by interpolating the node values ​​using the same interpolation function as used for Cartesian coordinates. As an example, consider an infinite element whose interior surface is a circular arc, discretized using a linear element, such as Figure 3 As shown. With respect to nodes Γ1 and Γ2, both virtual sources coincide with the center of the arc. This means that between Γ1 and the virtual source The radial distance between coincides with the radius of the arc; The same is true. When the node value and Using the standard interpolation of formula (18), we can get However, for points inside the infinite element, the discrete approximation Γ along the envelope Γ h , the distance from the center of the arc is not This is equivalent to, for Γ h The point on the envelope is defined as the distance Γ from the envelope in the radial direction. The interpolated virtual source position at Figure 3 ). In other words, when using the interpolation formulas (18) and (22) to calculate and , the interpolated position of the virtual source does not coincide with the arc center O. This means that, in general, even when the virtual source coincides with the radiation center relative to all infinite edges, the relation (35) is valid only along the infinite edges. The impact of this fact on the quality of the approximation will be analyzed in Section 5.1.1.

[0251] Figure 3 Schematic diagram of infinite element is shown. The envelope Γ is a circle and is represented by a piecewise linear envelope Γ h Discretize. The virtual source relative to nodes Γ1 and Γ2 coincides with the arc center O. hThe interpolated virtual sources of other points on the arc do not coincide with the arc center O.

[0252] In the Astley-Leis formulation, the Petrov-Galerkin approach is used: the test function and the solution function belong to different function spaces; specifically, the test basis function is the complex conjugate of the shape function, multiplied by a weight factor w:

[0253]

[0254] Where the weight factor takes the following form:

[0255]

[0256] In contrast to the Bubnov-Galerkin method, where the test basis functions are equal to the shape functions, the Petrov-Galerkin method results in an asymmetric system matrix. However, this choice of test basis functions is particularly convenient: firstly, as a conjugate formulation, the Astley-Leis formulation results in a frequency-independent system matrix, which can then be easily used in the time domain; secondly, for p w >2, the standard quadrature rule can be used to evaluate the integral required to compute the system matrix. In

[39] , p is chosen w = 2 is reasonable in the physical field, although it is pointed out that any power p w >2 preserves the boundedness of the system matrix integral. w Regardless of the choice of , the conjugate weighted formulation has been shown to lead to a variational formulation in a weighted Sobolev space, whose solutions have been shown to exist and be unique

[12] ,

[45] .

[0257] 2.3 Overview of stability conditions

[0258] In this subsection, we analyze the conditions under which the Astley-Leis conception is stable. The analysis is an extension of the theory proposed by Astley and Hamilton

[28] . Specifically, in section 3.5 of

[28] , only the system matrix relative to the infinite element was considered; in contrast, here we include the case where the interior domain is discretized using finite elements. Furthermore, in section 3.5 of

[28] , only the case where the interior surface of the infinite element is a circular arc was considered; in this subsection, we extend the analysis to arbitrary convex shapes. In addition, to ensure stability, we will argue that it is not necessary to construct "zero mass" elements, but it is sufficient to require that the element mass matrix relative to the infinite domain is semi-positive definite.

[0259] As described in

[46] , the second-order system (8) is stable as long as the damping matrix, stiffness matrix and mass matrix are positive semidefinite. is defined to be positive semidefinite:

[0260]

[0261] For simplicity, we will use A ≥ 0 as a shorthand notation to indicate the positive semi-definiteness of the matrix A. Due to the contribution of infinite elements, the mass matrix cannot always be positive semi-definite, causing the system to become unstable.

[0262] First, we show that if, with respect to the mass matrix of an infinite domain, (where n IE is the number of Dof associated with the infinite element) is semi-positive definite, then the global mass matrix M is also semi-positive definite. Let us consider the general vector

[0263]

[0264] Among them, the subscript FE and IE denote the association with the Dof of the finite element and the infinite element respectively; the last inequality in equation (40) comes from: Assuming M IE ≥ 0, and when using the standard polynomial shape function, M FE ≥0, and the material properties of the medium in which the sound propagates are physically realizable.

[0265] Next, we show that as long as the element mass matrix For all infinite elements, M is positive semidefinite. IE ≥ 0. Let us consider the general vector

[0266] Among them, n E is the number of infinite elements, and the inequality comes from the assumption that Where e=1,...,n E .

[0267] We now analyze In the Astley-Leis formulation, when the test basis functions are considered as weighted complex conjugates of the shape functions, It takes the following form (see

[13] ):

[0268]

[0269] Where c is the speed of sound. To make it easier to analyze its determinism, we rewrite it into a simpler form. Specifically, when the Gaussian quadrature formula is used to perform the integration, expression (42) can be rewritten as

[0270]

[0271] in, is a matrix of polynomial functions, evaluated at Gaussian points:

[0272]

[0273] n g is the number of Gaussian points, is the mother element coordinate vector of the jth Gaussian point, and D e is a diagonal matrix whose jth component along the diagonal is

[0274]

[0275] in, is the jth Gaussian weight, detJ j is the Jacobian of the transformation between physical coordinates and matrix coordinates, evaluated at the jth Gaussian point.

[0276] Due to D e is a diagonal matrix, so it is positive semidefinite if and only if all its components are not less than zero. This happens if and only if:

[0277]

[0278] We prove that if D e ≥0, then M IE ≥ 0. Let us choose a general vector y∈R n :

[0279]

[0280] The last inequality in (47) comes from the fact that, assuming D e ≥0, where It can be seen that if for any point in the infinite element Then the second-order system (8) is guaranteed to be time stable.

[0281] Let's transform the vector is defined as a vector pointing in a direction orthogonal to the envelope Γ, and its magnitude is

[0282]

[0283] Among them, n e (s,t) is in the mother coordinates The unit vector perpendicular to the envelope Γ at the point marked. For the classic Astley-Leis conception, the phase distance Defined in formula (29); After differentiation, we get

[0284]

[0285] In addition, for the characteristics of the gradient operator, due to the relationship between the gradient of a function and its directional derivative, we have:

[0286]

[0287] Let us first consider the ideal scenario, where the envelope Γ is exactly represented and the geometric discretization does not introduce any errors, and let us choose the infinite edge in the e-th element to stretch orthogonally to the envelope Γ. In this case, The physical coordinates are The distance between the point and Γ. In this way, there is a constant The line is perpendicular to n e (s,t). This means that, for the characteristics of the gradient operator, Perpendicular to Contour lines, parallel to n e (s,t), and thus with Therefore, we can write

[0288]

[0289] Finally, inserting equations (50) and (49) into equation (51) and simplifying the resulting equation, we obtain

[0290]

[0291] This means that if for all infinite edges the stretching direction is orthogonal to the envelope Γ, then we obtain anywhere in the external domain However, due to geometric approximations, small numerical errors can occur, which may cause some points to appear Astley and Hamilton

[28] suggested that by forcing M IE = 0 (which is equivalent to imposing ) to compensate for this effect, arguing that as long as the lateral node spacing is fine enough to fully resolve the curvature of Γ, the resulting change in the frequency domain solution can be neglected. As shown earlier in this section, in fact, forcing M IE ≥0 is sufficient, which is equivalent to forcing This can be done by only The point forces And The point To achieve unchanged.

[0292] The requirement that the infinite elements be stretched orthogonally to Γ means that, in general, the virtual source cannot be located at the center of the smallest sphere that encloses all sound sources. This limits the accuracy of all standard mapping concepts, especially for envelopes with high aspect ratios, and explains the need to develop new flexible infinite element concepts.

[0293] 3. Flexible and infinite concept

[0294] In this section, we describe a new flexible formulation where the geometry and field shape functions are constructed using two different mappings. This feature makes the accuracy of our formulation independent of the geometry of the infinite edges. Therefore, the flexible infinite elements can be used to accurately model transient acoustic wave propagation problems with arbitrary convex envelopes.

[0295] 3.1. Modified shape function

[0296] We remember that we used r to represent the distance from the center of the radiation and represents the interpolated distance from the virtual source. Our goal is to modify the shape function concept to obtain

[0297]

[0298] at any location inside the infinite element, regardless of the location of the virtual source. To do this, we first redefine the phase distance as

[0299]

[0300] Among them, r e According to formula (28), is the radial distance evaluated on Γ:

[0301]

[0302] One could try writing the shape function as

[0303]

[0304] in, will be obtained as the tensor product between a usual polynomial function in the tangential direction and a rational function in the radial direction, specifically designed to satisfy Eq. (53). However, constructing these rational functions is not trivial, and the effort to implement such formulas in existing FEM software would be considerable. Instead, we introduce a new auxiliary mapping between the radial distance r and the auxiliary coordinate v:

[0305]

[0306] It should be noted that this mapping is largely similar to Eq. (25), using v, r eand replace and The schematic diagram of auxiliary mapping is as follows Figure 4 shown.

[0307] Using the auxiliary mapping, we propose a new conception of shape function:

[0308]

[0309] Here, ψ is a common polynomial expressed with respect to variables s, t, and v, rather than the variables s, t, and v in equation (30) as conceived by Astley-Leis.

[0310] Figure 4 Schematic diagram showing auxiliary mapping in flexible infinite elements.

[0311] As an example, for ease of representation, let us consider a first-order Lagrange polynomial in two dimensions:

[0312]

[0313] We emphasize that the only difference between (31) and (59) is is replaced by the auxiliary variable v. Inserting the radial mapping (57) into (59), it can be verified that the following relationship holds:

[0314]

[0315] Expanding the above argument to higher order polynomials in three dimensions, we obtain:

[0316]

[0317] where m is the radial order chosen in the exterior domain. This is true not only for points along the infinite edge, but also for any point in the interior of the infinite element.

[0318] In other words, the multipole expansion (9) truncated at the mth order belongs to the space spanned by the columns of the mth order shape function. The concept of virtual sources does not apply to our formulation, since the relation (61) is true regardless of the geometry of the infinite element. This feature makes our formulation more flexible than the conventional mapping formulation, since the accuracy does not depend on the positioning of the mapping nodes. As in the Astley-Leis formulation, the test function is chosen to be the weighted complex conjugate of the shape function:

[0319]

[0320] 3.2. Time-domain stability of the modified infinite element

[0321] Because μ in formula (54) e The definition of is the same as in formula (29) The definition of is different, so the new flexible infinite element requires a separate stability analysis for the standard mapping concept. In particular, even when the infinite element is stretched orthogonally to the envelope Γ and the geometric discretization is not approximated, it is not guaranteed to be This poses a challenge to the stability of the method.

[0322] First, we analyzed The behavior of , shows that when the infinite edge is orthogonal to the envelope Γ, in the near field In the middle and far fields We then propose strategies to stabilize temporal constructs without compromising accuracy.

[0323] 3.2.1. Behavior

[0324] In this section, we study the case where the infinite edge is stretched in a direction orthogonal to the envelope Γ. The behavior of is similar to the time-stable infinite element proposed in

[28] . For simplicity, we consider the two-dimensional case, but all our considerations apply to three dimensions.

[0325] Figure 5 Schematic diagram showing the internal region of a flexible infinite element. The infinite edge is stretched orthogonally to the envelope Γ.

[0326] Let us consider the points on the envelope Γ whose coordinates are First, we show that as long as the infinite edges in the element are aligned with the normal vectors of the envelope Γ at the mapped nodes on Γ, then For the characteristics of the gradient operator, due to the relationship between the gradient of a function and its directional derivative, we have:

[0327]

[0328] The envelope Γ is the isovalue of μ, and μ is defined as 0. It follows:

[0329]

[0330] Because the gradient of μ is perpendicular to the contour of μ. For the standard Astley-Leis conception, we have already proved that relation (52) holds. Let us define

[0331]

[0332] as well as

[0333]

[0334] and μ about The partial derivatives of can be written as

[0335]

[0336] as well as

[0337]

[0338] According to the triangle inequality (see Figure 5 ), we have

[0339]

[0340] It follows from this

[0341]

[0342] As a result, we finally get

[0343]

[0344] Again, as in the classic Astley-Leis formulation, small errors may occur due to geometric approximations. In this case, as in

[28] , we propose to ignore these small differences to ensure the stability of the method.

[0345] for It is easy to conclude This means

[0346]

[0347] exist and In the area between The exact behavior of depends on the geometry of the infinite element. In practice, we often encounter The geometry of Figure 6 ).

[0348] 3.2.2. Strategies for stabilizing the system

[0349] Due to the analysis in the previous section The behavior of , requires a strategy of stable and flexible infinite element conception. In the unstable version, the element-level mass matrix contribution is written as:

[0350]

[0351] To stabilize new flexible ideas, we recommend using Replace the unstable mass matrix in

[0352]

[0353] in

[0354]

[0355] When using Gaussian quadrature to evaluate equation (74), we can write

[0356]

[0357] The jth component of the diagonal matrix is

[0358]

[0359] This ensures And thus ensure that the concept is time-stable.

[0360] As an example, we consider the two-dimensional case where the FE boundary is the envelope of an ellipse with width 6 and height 2. Figure 6 For the case where the envelope is stretched orthogonally and from the coordinate The mapping point We will and Plotted as a function of v. Figure 6 Show and An infinite edge is stretched orthogonal to the envelope of the ellipse (with an aspect ratio of 3) at the mapping node (0.5, 0.986).

[0361] 3.2.3. Accuracy of Stable Concept

[0362] In the previous section, we demonstrated how to stabilize the flexible infinite element concept. However, the stabilization procedure may potentially affect the accuracy of the method. This is because Some of the items in D e The corresponding entries in are different, so become A bad approximation of .

[0363] Using the relation (71) that holds when the infinite edge is stretched orthogonally to the envelope r, we can devise a strategy to ensure the accuracy of the flexible infinite element after the stabilization procedure. We propose to adjust the power p in the weight factor w w , to ensure that the stabilization procedure does not significantly affect accuracy. In practice, this is achieved by controlling With D eThe smaller the difference, the smaller the loss in accuracy. We recall that the weight factor in the Astley-Leis conception is (38), where p w = 2. We observe that w = 1 at the envelope, and w → 0 as By increasing p w , we make The weighting factor is even smaller. Let us consider the jth Gaussian point g κ ,in (We have already mentioned that this usually happens ). Matrix D e The corresponding diagonal entries in are

[0364]

[0365] and The corresponding entry in is

[0366]

[0367] We emphasize that since relation (71) holds for orthogonal stretching of infinite edges, g κ cannot lie on the envelope Γ. For this reason, it is ensured that w(g κ )<1 and we can force

[0368]

[0369] By choosing

[0370]

[0371] Where ∈ is the tolerance value, which is set by the user. w The constant is chosen over the entire domain to ensure the continuity of the test function over the entire domain, so that Eq. (81) is satisfied for all infinite elements. The smaller the tolerance ∈, the greater the error D. e and The smaller the difference between , and hence the smaller the accuracy loss after the stabilization step. We emphasize that stability is guaranteed by (75), while the accuracy of the stabilization concept is guaranteed by (81). In Appendix A, we present bounds on the accuracy loss due to the stabilization procedure, relating it to the chosen tolerance ∈. Choosing p w The only cost of > 2 is that a larger number of Gaussian points are required to perform the calculations, since the order of the polynomial to be integrated increases with p. wIncrease. Typically, for numerical acoustic simulations using FEM, the most computationally intensive operation is the factorization of the fully assembled matrix; for this reason, despite the increased computational cost of the assembly step, the total computational cost of the flexible conception is comparable to that of the Astley-Leis conception for the same number of Dofs. Additionally, in many cases, as we will document in Sections 5 and 6, flexible infinite elements, due to their geometric flexibility, lead to a drastic reduction in the number of Dofs, making the simulation less computationally intensive.

[0372] 3.3. Mesh stretching in the normal direction

[0373] As already explained in the previous section, the infinite edges need to be stretched orthogonally to the envelope Γ to ensure stability in the time domain. Unfortunately, computing the normals of the envelope Γ for a generally convex shape is not trivial. In some cases, information about the normal vectors can be provided by the CAD tool, with the drawback that an additional preprocessing step is required. In other cases, only a meshed domain is available and no reference geometry is available. We address this issue by following the workflow proposed in

[47] . With this approach, only the mesh of the inner domain is required. The required normal vectors at the mapping nodes on Γ are calculated as follows:

[0374] For a vertex, the normal vector n e (s,t) is taken as the average of the element normals touching the vertex.

[0375] ·Inside a quadratic element, for nodes defined in the middle of an edge (or face), n e (s,t) is calculated as the average of the normal vectors at the vertices belonging to this edge (or face).

[0376] The user must specify the stretch length, which for the i-th mapping node on Γ is the distance defined as (12) Inside the e-th element, its parent coordinate is (s i ,t i ,-1). Then, the mapping node inside the external domain The location is determined by:

[0377]

[0378] This process is repeated for all mapped nodes on the envelope Γ and for all infinite elements.

[0379] 4. Software implementation of high-order finite element and infinite element methods

[0380] In this subsection, we discuss the implementation of infinite elements as high-order FEM codes. The considerations made here apply both to the Astley-Leis conception and to flexible infinite elements. We refer to the adaptive order p-FEM method described in

[48] , where Lobatto shape functions (also called integrated Legendre polynomials) are used, and to

[49] , where the method is generalized to include anisotropic orders. For many Helmholtz problems, this approach has been shown to greatly improve the efficiency of the numerical model when compared to conventional FEM, exploiting the fact that high-order polynomial approximations are more effective in controlling contamination effects. When a solution in the frequency domain is required over a certain frequency range, an additional benefit is found in the use of a hierarchical set of shape functions (such as Lobatto polynomials). For the highest order, the element matrices can be computed only once, and then, at each frequency, one can extract only the required parts of these matrices to assemble the global matrix. The key feature of the proposed method is the use of a simple local error indicator to a priori select the polynomial order in each element. This approach requires only a single grid for calculations over a range of frequencies, and the error estimator is able to select the appropriate polynomial order in each element to maintain efficiency and accuracy as frequency varies by identifying the lowest order required for the accuracy target.

[0381] Infinite elements based on Lagrange polynomials have been shown to have high radial order regulation problems

[35] . Jacobi

[35] and Berstein

[50] polynomials have been proposed to improve regulation and optimize the performance of Krylov based iterative solvers. As pointed out in

[35] , Lobatto polynomials, although suboptimal from a regulation point of view, outperform Lagrange polynomials in this respect.

[0382] By adopting the same type of polynomials for both finite and infinite elements, software implementation becomes particularly easy. In this case, infinite element subroutines can be developed by reusing most of the routines called by their finite element counterparts. The only differences are the need to implement a new radial geometry mapping (23) and the need to use a different formula to calculate the system matrix, which does not require any additional implementation effort. In other words, in the software implementation, infinite elements can simply be regarded as specialized anisotropic high-order elements. This is an important aspect because the implementation burden of the infinite element method may hinder its widespread use in practical applications. An additional benefit of using Lobatto polynomials in the external domain as well is that the hierarchical structure of the system matrix is ​​retained for unbounded domain simulations and can be used for frequency domain calculations.

[0383] The order of the shape function inside the FE domain is chosen based on the a priori error estimator described in

[48] ; the radial order on infinite edges is chosen by the user, while the order of the shape function on the envelope Γ is designed to ensure continuity with the FE domain; finally, the Dof of vertices, edges, and faces is set to zero at infinity (the ) to force the decay of the pressure field.

[0384] 5. Numerical Test

[0385] In this section, flexible infinite elements are compared with the Astley-Leis formulation and numerical examples are provided in two and three dimensions in both the frequency and time domains. The mesh is generated in the FE domain using the mesh generator Gmsh

[51] , while the infinite elements are stretched at runtime using one of the following strategies:

[0386] Confocal-ray method, where the infinite edges stretch in radial direction and converge towards the radiation center. In the Astley-Leis conception, this means that all virtual sources coincide with the radiation center.

[0387] • Normal-ray method, where the infinite edge is aligned with vector n, perpendicular to the FE boundary envelope Γ.

[0388] Curved triangular or tetrahedral elements are used in the inner domain, while the outer mesh consists of quadrilateral or prismatic elements. The results are obtained using an in-house Matlab code.

[0389] In the following examples, we will often refer to the numerical solutions u respectively calculated on the envelope Γ num With reference solution u ref The relative L 2 error:

[0390]

[0391] And in the internal domain (FE domain):

[0392]

[0393] When using the flexible formulation for transient problems, we measure the error caused by the stabilization procedure in the frequency domain for the maximum frequency of interest (see discussion in Section 3.2.3) as follows:

[0394]

[0395] in, is obtained using a stable concept, while is obtained using an unstable concept.

[0396] 5.1.2D Monopole Source

[0397] First, we analyze the case of a monopole source in the 2D domain, for which the analytical solution is known:

[0398]

[0399] in, is a second-class Hankel function, and the monopole source position coincides with the radiation center. The source is surrounded by a convex envelope, and the two meshes generated using the confocal ray method and the normal ray method are as follows Figure 7 As shown. Near the source, the mesh is refined due to the singularity of the analytical solution at r = 0, while larger elements are used in the rest of the domain to fully exploit the efficiency of the higher-order FEM. The element characteristic length at the interface between the inner and outer domains is h = 0.25.

[0400] Figure 7 The 2D mesh used for the numerical experiments is shown. The FEM mesh, the envelopes Γ and Σ (in blue), and the portion of the infinite edges between Γ and Σ are shown. In the confocal ray method, all virtual sources coincide with the radiation center. In the normal ray method, all infinite edges are locally orthogonal to the envelope Γ. In both cases, the mesh is refined near the monopole source location, which coincides with the radiation center.

[0401] 5.1.1. Frequency domain results

[0402] For the following frequency domain results, the order of the shape function is fixed (in both the inner and outer domains). Numerical convergence studies are done by investigating the influence of the chosen order on the accuracy of the method. First, we compare the performance of the new flexible infinite element with the Astley-Leis formulation for a confocal ray configuration. Even though for the considered configuration all virtual sources coincide with the radiation center, the Astley-Leis formulation and the flexible formulation still differ. This is because, for the Astley-Leis formulation, the radial coordinates are obtained using standard interpolation functions, as shown in Eqs. (18) and (22), while for the flexible infinite element they are computed as functions of the Cartesian coordinates using Eqs. (28) and (55). In Figure 8 In this paper, we compare the L on the envelope Γ of the two methods at several frequencies. 2 The error is plotted as a function of the order of the shape function.

[0403] For all four frequencies, the L2 error of the standard Astely-Leis conception reaches a plateau, while the flexible infinite element method shows better p-convergence properties. This result can be explained by the fact that for the flexible infinite element, Eq. (53) is satisfied everywhere, while for the Astley-Leis conception it is only valid along the infinite edges when the confocal ray extrusion method is employed (see discussion in 2.2.2). In other words, the truncated multipole expansion (9) belongs to the space spanned by the shape function in the entire external domain for the flexible conception, but only along the infinite edges for the Astley-Leis conception. The resulting accuracy difference is higher for coarser meshes (which are usually used in the context of high-order adaptive FEMs), but can be neglected when fine meshes are used.

[0404] Figure 8 Results for the confocal ray configuration are shown for the flexible infinite element (black) and the standard Astley-Leis element (grey). All virtual sources coincide with the locations of the radiation centers and physical sources.

[0405] Next, we compare the flexible infinite element with the Astley-Leis infinite element when the normal ray stretching strategy is adopted. This is a preliminary comparison in the time domain, where the normal ray strategy is the only one that can provide stable results. Fig. 9 In , we documented a significant improvement in accuracy due to the flexible conception. In this case, the results can be explained by the fact that for the flexible conception, (53) is again satisfied everywhere, since the space spanned by the shape function does not depend on the orientation of the infinite element; in contrast, for the Astley-Leis infinite element, (53) is not satisfied, since the virtual source does not coincide with the radiation center.

[0406] Fig. 9 Results are shown for the flexible infinite element (black) and the standard Astley-Leis element (grey), normal ray configuration. The mesh stretching is performed as described in Section 3.3.

[0407] 5.1.2. Time domain

[0408] It is well known that confocal ray stretching methods often lead to unstable formulations in the time domain, as discussed in Section 2.3. For this reason, we focus only on the normal ray method. For this example, the monopole source is placed at (0.5, 0.4) so ​​that the flexible infinite element works 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 = 20, filtered using a Hamming window. Therefore, the force vector f τ (τ)(T∈[0,T]) is

[0409]

[0410] where f is the force vector in the frequency domain. The order of the shape function in the FEM domain is adaptively selected using the a priori error indicator of Bériot et al.

[48] , where the target L 2 The error E2 = 0.01, and for both cases, the radial order in the infinite domain is set to m = 6. The flexibly conceived power in (38) is set to p w =6, and the error caused by the stabilization procedure is According to equation (85), k = 20, i.e. the maximum wave number of interest for the analyzed problem. In other words, it is clear that the error introduced by the stabilization procedure can be neglected since it is significantly lower than the target error for the order of the shape function calculated in the FEM domain.

[0411] The reference numerical solution is calculated on a larger circular FEM domain, centered at the monopole source location, with a radius of 2, using Astley-Leis infinite elements, with infinite edges stretched in the radial direction, with a stretch length of All virtual sources are made to coincide with the positions of the monopole sources. The radial order in the infinite domain of the reference solution is set to 10. The radial quantities are calculated according to Eqs. (28) and (55) to improve the accuracy of the reference solution (see the discussion in Section 5.1.1). In other words, the radial quantities are not calculated by interpolating the radial distances relative to the node values, but rather taking the distances of the interpolated points in the infinite element from the radiation center. This is equivalent to defining the interpolated positions of the virtual sources to coincide with the radiation center, thus obtaining better accuracy (see the discussion in Section 2.2.2). Due to this modification, in this specific configuration, the Astley-Leis conception is equivalent to the flexible conception. The accuracy improvements due to this modification have been recorded in Section 5.1.1. The reference numerical solution in the time domain is calculated by constructing “zero mass” infinite elements according to the procedure of Astley and Hamilton

[28] , which does not significantly affect the accuracy on the circular domain. The use of a large FEM domain, the implementation of the modified interpolation of the radial quantities, and the fact that all virtual sources coincide with the positions of the physical sources ensure good accuracy of the reference model.

[0412] exist Fig.10 In Figure 1, the pressure histories of two points on the convex envelope (point A = (-1.3, 0.2) and point B = (-0.78, -0.24)) are plotted as a function of time for the Astley-Leis conception and the flexible conception together with the reference solution. The numerical solution obtained using the flexible conception closely matches the reference solution, while the Astley-Leis conception does not provide satisfactory results. Fig.10The performance of the flexible conception of arbitrary convex envelopes in the time domain and the Astley-Leis conception is shown using the normal ray stretching strategy (see Figure 7 ). Pressure history at points A and B on the convex envelope. Comparison between the reference (black line), the flexible infinite element (grey dashed line), and the traditional Astley-Leis conception (red line). The radial order in the outer domain is set to m=6, while the order in the inner domain is determined by an a priori error estimator.

[0413] 5.2. Scattering of plane waves by a sphere

[0414] In this subsection, we consider the scattering of plane waves by a sphere, and we test the performance of flexible and Astley-Leis infinite elements. Several polyhedral shapes with rounded edges are used as closed surfaces. Regardless of the infinite element concept used, the accuracy of the method requires rounding of the corners, since the normal vector must be continuous over the entire surface so that the stretching direction is orthogonal to the envelope at any point.

[0415] The scattering sphere has unit radius, while the pseudopolyhedral domain has a midradius of 2 and a fillet radius of 0.2. The mesh used in the interior domain has elements of characteristic length h = 0.3, the fillets are refined to resolve geometric details, and the order in the FEM domain is again set adaptively using the a priori error indicator of Bériot et al.

[48] , where the target L 2 Error E2=0.01.

[0416] 5.2.1. Frequency domain results

[0417] The reference analytical solution to this problem is:

[0418]

[0419] Among them, a sca is the radius of the scattering sphere, j m is the mth-order spherical Bessel function, is a spherical Hankel function of the first kind of order m, and P m is the mth order Legendre polynomial.

[0420] The optimal interpolation error is calculated as the difference between the analytical solution and its L 2 The relative L between the projections onto the finite element space 2 error:

[0421]

[0422] This error corresponds to the best numerical solution that can be achieved given the FE mesh and the order used, regardless of the infinite element concept used. Fig.11 In this paper, we demonstrate numerical results obtained using flexible infinite elements stretched in the normal direction. We perform several simulations, varying the radial order m in the outer domain from 1 to 10; the shape function order in the FEM domain is chosen using an a priori error estimator, as already explained at the beginning of this section, and remains constant between the different simulations. It can be observed that for both frequencies, the inner domain E2(Ω i ) in the relative L 2 The error is close to the optimal interpolation error for increasing radial order m. It can be seen that the pseudo-tetrahedron is the most challenging shape because it is the closest to the scattering surface.

[0423] Fig.11 Each colored line is shown to correspond to a pseudo-polyhedron domain. Darker colors correspond to pseudo-polyhedrons with more edges. Dashed line: optimal interpolation error in the FE domain. The order of the shape function in the inner domain is adaptively set using an a priori error estimator and remains constant between different simulations. The radial order m is consistent throughout the outer domain and varies between different simulations.

[0424] 5.2.2. Time domain results

[0425] In this subsection, of all the pseudo-polyhedral shapes proposed in Section 6.3.1, only the most challenging one, the pseudo-tetrahedron, is considered. We now consider the results of transient simulations where the input signal is a sine wave with k = 10, filtered by a Hamming function, as in the 2D case. Fig.12 In , we show the pressure results at envelope points A and B for the flexible infinite element and the Astley-Leis infinite element; for both cases, the radial order in the infinite element is set to 5, while the order in the interior domain is determined by the a priori error estimator, as explained at the beginning of this section. The power of the flexible concept in (38) is set to p w =6, and the error caused by the stabilization procedure is According to (85) k = 10, i.e. the maximum wave number of interest for the analyzed problem. Again, it is clear that the error introduced by the stabilization procedure can be neglected since it is significantly lower than the target error used to calculate the order of the shape functions in the FEM domain.

[0426] As in the 2D case, the reference numerical solution is computed using Astley-Leis infinite elements stretched in the radial direction with infinite edges on a larger FEM domain consisting of a sphere of radius 2 surrounding the scatterer so that all virtual sources coincide with the center of the scattering sphere. The radial order in the infinite domain of the reference solution is set to 10 and the radial quantities are again computed according to expressions (28) and (55). The reference numerical solution in the time domain is computed by constructing “zero mass” infinite elements following the procedure of Astley and Hamilton

[28] . As in the 2D case, the good accuracy of the reference numerical model is ensured by the use of a large FEM domain, a high radial order m and a modified interpolation of the radial quantities, combined with the fact that all virtual sources coincide with the center of radiation.

[0427] Also for this case, it can be noted that the flexible conception outperforms the traditional Astley-Leis conception.

[0428] Fig.12 The performance of the flexible and Astley-Leis formulations of pseudo-tetrahedrons enclosing a scattering sphere in the time domain using a normal ray stretching strategy is shown. Pressure history at points A and B on the convex envelope. Comparison between the reference (black line), the flexible infinite element (dashed grey line), and the traditional Astley-Leis formulation (red line). The radial order in the outer domain is set to m=6, while the order in the inner domain is determined by an a priori error estimator.

[0429] 6. Application in industrial complexity problems

[0430] In this section, we propose a workflow for solving infinite problems of industrial complexity using the flexible infinite element concept. Specifically, to demonstrate the effectiveness of the approach, we will consider a shark submarine, such as Fig.13 As shown. For frequency domain validation, we will simulate the scattering of plane waves, while for the time domain, we will consider monopole sources located near the FE boundary. An automatic algorithm is used to create a convex envelope around the scatterer and to automatically stretch the infinite edges in the normal direction. To speed up the transient simulations, we use the MOR procedure proposed by van de Walle

[25] and van Ophem

[26] , called the Automatic Krylov Subspace Algorithm (AKSA). The high degree of automation of the proposed workflow, as well as the computational efficiency achieved by using a convex envelope of generic shape, makes it very attractive for realistic transient applications.

[0431] Fig.13 A shark model is shown.

[0432] 6.1. Automatic generation of convex domains

[0433] In this subsection, we describe the automatic procedure for obtaining the computational domain mesh. First, a surface mesh consisting of second-order triangular elements with a characteristic length of h = 0.02 is created on the submarine. Second, an approximate convex hull is automatically generated using the quickhull algorithm

[36] , which is implemented in Simcenter 3D

[37] , a simulation package developed by Siemens Industry Software. The minimum distance from the submarine surface to the convex hull is d = 0.05. Third, the outer surface is meshed using secondary triangular elements with a characteristic length of h = 0.02. Then, a volume mesh consisting of tetrahedral secondary elements is created between the two surface meshes. The infinite elements are automatically stretched along the normal direction so that the transient concept is time-stable, following the procedure described in Section 3.3. The final mesh is as follows Fig.13 As shown, it consists of 183519 quadratic elements and 232959 nodes. Since the convex envelope lies inside the smallest sphere enclosing the scatterer, theoretically the Atkinson-Wilcox theorem cannot guarantee convergence. In Sections 6.3.1 and 6.3.2, we will show that in practice, when flexible infinite elements are attached to the convex hull, good accuracy can be achieved in both the time and frequency domains.

[0434] 6.2. Model Reduction Procedure

[0435] Since flexible infinite elements, like traditional Astley-Leis infinite elements, produce 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

[25] ,

[26] algorithm, which aims to automatically construct a reduced-order model whose frequency response function (FRF) matches that of the initial full-order model (FOM) within a user-specified tolerance over the frequency range of interest by using only a few system solutions and enriching the projection basis with high-order moments. To obtain a time-domain simplified model of the shark submarine, for example, valid between k = 5 and k = 50 with a relative error tolerance of tol = 0.01, only five full system solutions are required.

[0436] Results

[0437] 6.3.1. Frequency domain results

[0438] We now test the performance of flexible infinite elements for plane wave scattering in the frequency domain. Simulations are performed for four plane waves propagating at different angles of incidence α and θ. As a reference, we use the results obtained using an automatic perfectly matched layer (AML) implementation

[47] with 5 layers. For the configuration using infinite elements, the order in the FEM domain is adaptively assigned using the a priori error estimator of Bériot et al.

[48] , where the target L 2The error E2 = 0.01, while the radial order in the outer domain is user defined; for the AML case, the order is adaptively distributed throughout the computational domain. For all analyzed cases, the mesh used is Fig.13 For k = 10 and k = 100, the relative L on the convex envelope is 2 Error Fig.14 In all cases, the infinite element shows good convergence, and the relative L of radial order 8 2 The error is always below 0.03, close to the target error. Note that in some cases, the relative L 2 The error is even lower than the target error; this can happen because the adaptive rule for determining the order in the inner domain is based on an inaccurate estimate of the error. Fig.14 Frequency domain results are shown for a plane wave hitting a shark submarine at different spherical angles α and θ (in radians). Performance of the flexible infinite element concept with infinite edges stretched in the normal direction. The AML implementation is used to compute the reference solution.

[0439] 6.3.2. Time domain results

[0440] Finally, we consider the scattering from a monopole source at (0,0,0.5) placed in front of the shark.

[0441] The objective L for calculating the order of the shape function in the FEM domain 2 The error is again set to E2 = 0.01. The radial order of the infinite element is set to m = 8, fixed throughout the outer domain. The flexibly conceived power in (38) is set to p w =8, and the error caused by the stabilization procedure is According to (85) k = 50, i.e. the maximum wave number of interest for the analyzed problem. Again, it is clear that the error introduced by the stabilization procedure can be neglected since it is significantly lower than the target error used to calculate the order of the shape functions in the FEM domain.

[0442] The input signal consists of a sine wave with k=45, filtered using a Hamming window. AKSA is used to generate a valid ROM between k=5 and k=50, with a relative error tolerance of tol=0.01. Fig.15 Compare the FOM and FRF of the ROM for a random entry in the solution vector. The final ROM has a Dof of 90, while the original FOM has a Dof of 465905.

[0443] Fig.15 The FRF is shown comparing the ROM (grey line) and the FOM (black stars) for random entries in the solution vector. Results obtained by appending flexible infinite elements to the convex hull.

[0444] We ran all our numerical experiments on a Dell desktop computer with 384Gb of RAM and a 3.00GHz clock frequency. The offline MOR procedure took 5 minutes and 45 seconds, and a single full system solution (an operation performed multiple times during AKSA) required 15.5Gb of storage. For reference, we used the standard Astley-Leis formulation to compute the transient response to the same excitation, with a spherical envelope around the submarine. In this case, the model had 1039249 nodes, 801405 elements, and 2515697 Dof, which is 5 times that of the convex hull model. The offline MOR procedure took 95 minutes, more than in the previous case, and required 151.9Gb of storage for the full system solution; in addition, the projection basis that needs to be stored during the entire MOR process is 5 times larger than before, because the number of rows corresponds to the Dof of the FOM. As in the previous case, the final ROM has 90 Dof. In Table 1, we summarize the data on the computational requirements of building ROMs for both models using AKSA.

[0445]

[0446] Table 1: Computational requirements of AKSA, comparison between the model with a convex hull and a spherical envelope surrounding the submarine.

[0447] Fig.16 The pressure history at points A and B on the submarine surface is shown using Astley-Leis infinite elements attached to the spherical envelope (black line) and flexible infinite elements attached to the convex hull (grey dashed line).

[0448] exist Fig.16 In , we plotted the pressure at different points on the submarine surface as a function of time for both the convex hull and spherical envelope models. The use of the convex hull model, enabled by the newly developed flexible infinite element concept, is clearly more efficient, as the cost of the MOR process is much lower, while the accuracy is similar to that obtained using the spherical envelope. Fig.17 In Figure 1 we show the visualization of transient acoustic waves at different snapshots using Simcenter 3D (version 2022.1).

[0449] Fig.17 Shows the transient results visualized in Simcenter 3D. The monopole source is placed in front of the shark.

[0450] 7. Conclusion

[0451] In this work, we introduce a novel infinite element formulation, called flexible infinite elements, for modeling unbounded wave propagation problems. One of the main drawbacks of the conventional mapping formulation is that the accuracy depends on the geometric configuration of the infinite elements. For transient simulations where the infinite elements must be stretched orthogonally to the FE boundaries, the user cannot freely position the virtual sources at the desired locations. This situation limits the accuracy of such formulations, especially when using envelopes with high aspect ratios. In the flexible formulation, a new set of shape functions is used, which makes the accuracy independent of the location of the mapping nodes. The newly developed elements are able to be used for transient simulations with arbitrary convex FE boundaries. The strategies to guarantee the stability of the proposed method in the time domain while maintaining the accuracy are discussed in detail. The accuracy of the flexible formulation in both the frequency and time domains is tested through several examples. It is demonstrated how the new formulation can be easily implemented in software based on high-order adaptive FEM. Finally, a workflow including the automatic stretching of infinite elements and the automatic order reduction of the resulting model is proposed. Since high-order polynomials are often used for test functions, the assembly step is more expensive than the standard Astley-Leis infinite elements (see discussion in Section 3.2.3) in order to guarantee accuracy for the stable transient flexible concept. However, for most cases of interest, the most computationally expensive operation is the factorization step. In particular, it is demonstrated that computational savings can be achieved by using flexible infinite elements, especially for high-aspect-ratio scatterers. For future research, it would be of interest to develop new error estimators to automatically determine the radial order of the shape function in the infinite domain. Finally, the automatic selection of the p in the test function for transient simulations, based on the user's accuracy requirements, could be investigated. w The new strategy will be of great interest.

[0452] Appendix A. Error bounds for stable and flexible concepts

[0453] In this appendix, we analyze the impact of the stabilization procedure on the accuracy of the solution p to the linear system of equations (7). Specifically, we propose error bounds that relate the accuracy loss to the tolerance ∈ chosen in equation (80). Let A full Fully Assembled Matrix

[0454] A full : = [-ω 2 M+iωC+K] (A.1)

[0455] and let the disturbed system be

[0456] (A full +ΔA full )(p(ω)+Δp(ω))=f(ω), (A.2)

[0457] Among them, △A full A is introduced by the stable program fullis the perturbation in , and ∆p(ω) is the subsequent perturbation in the solution vector. As described in

[52] , we have

[0458]

[0459] Among them, the matrix A∈R a×b The ∞ norm of

[0460]

[0461] And the condition number κ ∞ Defined as

[0462] κ ∞ (A full ):=||A full || ∞ ||A full -1 || ∞ (A.5)

[0463] Since the stable program includes the stable version Replace the infinite element mass matrix contribution M IE , so we have

[0464]

[0465] It is interesting to note that due to ||ΔA full ||∞ and ω 2 is proportional to , so the loss of accuracy is greater for the highest frequencies of interest. At the element level, we can write

[0466]

[0467] matrix The maximum entry of is at most ∈, which is imposed by equations (80) and (81). Therefore, the matrix The largest entry is less than Because when the Lobatto polynomial is chosen, and all its entries are less than or equal to 1. Constructing the global matrix M IE The element matrix of the entries in The maximum number is n max ,in

[0468]

[0469] Among them, n Γ is the number of mapping nodes on the envelope Γ, and is the number of elements that share the i-th node. Therefore, we have the matrix The maximum number of entries in is at most Using the ∞ norm definition, we can write

[0470]

[0471] in

[0472]

[0473] in, Defined as the matrix M IE The number of nonzero entries in the i-th row of . By substituting inequality (A.9) into (A.6), we obtain

[0474] ||ΔA full || ∞ ≤(wn g ) 2 n nz n max ∈, (A.11)

[0475] Substituting it into (A.3), we get

[0476]

[0477] Expression (A.12) represents the connection between the tolerance ∈ set by the user at the element level and the loss of accuracy from the stabilization procedure.

[0478] Additional feasible implementations or alternative solutions of the present invention also include combinations of features described above or below with respect to the embodiments (not explicitly mentioned herein). Those skilled in the art can also add individual or isolated aspects and features to the most basic form of the present invention.

[0479] It is part of the invention that not all steps of a method must be executed on the same component or computer instance, but can also be executed on different computer instances.

[0480] In view of the description and embodiments below that will be described in more detail in the context of the accompanying drawings, the attributes, features and advantages of the present invention described above and the modes for realizing them will become clearer and more easily understood. The following description is not intended to limit the present invention to the embodiments included. In different figures, the same parts or parts can be marked with the same reference numerals. Generally, these figures are not drawn to scale. It should be understood that the preferred embodiments of the present invention can also be any combination of the present invention or the above-described embodiments with the present invention.

[0481] These and other aspects of the invention will be apparent from and elucidated with reference to the embodiments described hereinafter. BRIEF DESCRIPTION OF THE DRAWINGS

[0482] Embodiments of the present invention will now be described, by way of example only, with reference to the accompanying drawings, in which:

[0483] Figure 1 Schematic diagram showing exterior acoustic problems;

[0484] Figure 2 A schematic diagram showing a mapped infinite element is shown;

[0485] Figure 3 A schematic diagram showing an infinite element;

[0486] Figure 4 A schematic diagram showing auxiliary mapping in a flexible infinite element;

[0487] Figure 5 A schematic diagram showing the interior region of a flexible infinite element;

[0488] Figure 6 (a) shows an unstable concept;

[0489] Figure 6 (b) shows a stable conception;

[0490] Figure 7 (a) shows the 2D grid and confocal rays used for numerical experiments;

[0491] Figure 7 (b) shows the 2D grid used for numerical experiments, normal rays;

[0492] Figure 8 (a) shows the results of the confocal ray configuration when k=1, the flexible infinite element (black) and the standard Astley-Leis element (grey);

[0493] Figure 8 (b) shows the results of the confocal ray configuration when k=5, the flexible infinite element (black) and the standard Astley-Leis element (grey);

[0494] Figure 8 (c) shows the results of the confocal ray configuration when k=10, the flexible infinite element (black) and the standard Astley-Leis element (grey);

[0495] Figure 8 (d) shows the results of the confocal ray configuration when k=25, the flexible infinite element (black) and the standard Astley-Leis element (grey);

[0496] Fig. 9 (a) shows the results of the normal ray configuration for the flexible infinite element (black) and the standard Astley-Leis element (grey) when k=1;

[0497] Fig. 9 (b) shows the results of the normal ray configuration for the flexible infinite element (black) and the standard Astley-Leis element (grey) when k=5;

[0498] Fig. 9 (c) shows the results of the normal ray configuration for the flexible infinite element (black) and the standard Astley-Leis element (grey) when k=10;

[0499] Fig. 9 (d) shows the results of the normal ray configuration for the flexible infinite element (black) and the standard Astley-Leis element (grey) when k=25;

[0500] Fig.10 (a) shows the use of the normal ray stretching strategy, point A = (-1.3, 0.2);

[0501] Fig.10 (b) shows the use of the normal ray stretching strategy, point B = (-0.78, -0.24);

[0502] Fig.11 (a) shows that when k = 1, each colored line corresponds to a pseudo-polyhedron domain;

[0503] Fig.11 (b) shows that when k = 10, each colored line corresponds to a pseudo-polyhedron domain;

[0504] Fig.12 (a) shows the use of the normal ray stretching strategy, point A = (1.137, 1.222, 0.560);

[0505] Fig.12 (b) shows the use of the normal ray stretching strategy, (b) point B = (0.522, 1.426, 0.642);

[0506] Fig.13 (a) shows a shark model, a CAD model;

[0507] Fig.13 (b) shows a shark model, surface mesh, characteristic length h = 0.02;

[0508] Fig.13 (c) shows the shark model, the convex envelope automatically generated at the minimum distance d=0.05;

[0509] Fig.13 (d) shows a cross-section of the final tetrahedral mesh of the shark model;

[0510] Fig.14(a) shows the frequency domain results of a plane wave hitting a shark submarine at different spherical angles α and θ (expressed in radians), k = 10;

[0511] Fig.14 (b) shows the frequency domain results of a plane wave hitting a shark submarine at different spherical angles α and θ (expressed in radians), k = 100;

[0512] Fig.15 The FRF shows the ROM (grey line) and FOM (black star) for random entries in the solution vector.

[0513] Fig.16 (a) shows the pressure history at point A on the submarine surface;

[0514] Fig.16 (b) shows the pressure history at point B on the submarine surface;

[0515] Fig.17 (a) shows the transient results visualized in Simcenter 3D, τ = 5;

[0516] Fig.17 (b) shows the transient results visualized in Simcenter 3D, τ = 205;

[0517] Fig.17 (c) shows the transient results visualized in Simcenter 3D, τ = 405;

[0518] Fig.17 (d) shows the transient results visualized in Simcenter 3D, τ = 605;

[0519] Fig.17 (e) shows the visualization of transient results in Simcenter 3D, legend (dimensionless pressure);

[0520] Fig.18 showing a representation of the underlying acoustic problem;

[0521] Fig.19 A flow chart showing the steps of a method according to the invention.

[0522] The illustrations in the drawings are in schematic form.

[0523] It is noted that in different drawings, similar or identical elements may be provided with the same reference numerals. DETAILED DESCRIPTION

[0524] Fig.18 A schematic representation of the underlying acoustic problem is shown. Ωi is the interior domain, and Ωe is the exterior domain. Γ is the interface between the two domains, and ΓS is the envelope in which the Sommerfeld condition is imposed (i.e. no wave can be reflected at infinity). The interface Γ is a convex shape CXS around the object OBJ, using the fast convex hull algorithm QHA. The control volume CVL of the unbounded domain UBD is limited by a dotted line to illustrate that the boundary is not reflective like a normal boundary. The object emits sound as a sound source SDS.

[0525] Fig.19 A simplified flow chart illustrating a computer-implemented method for acoustically modeling an object OBJ in an unbounded domain UBD as a numerical simulation model MDL to determine field variables FDV within a control volume CVL of the unbounded domain is shown. The method comprises the following steps:

[0526] (a) providing a coordinate system COS mapping the position POS in the control volume CVL;

[0527] (b) providing a geometrical specification of a sound source SDS in said control volume CVL based on said coordinate system COS;

[0528] (c) constructing a convex CXS around the sound source SDS using a fast convex hull algorithm QHA;

[0529] (d) Stretching the infinite element IFE starting from the convex CXS;

[0530] (e) By mapping the radial coordinates of the coordinate system COS to the auxiliary coordinates extending from the convex shape CXS to the infinite element IFE Transform the coordinate system COS into the field variable coordinate system FVC;

[0531] (f1) Use auxiliary coordinates To construct the field shape function FSF;

[0532] (f2) Applying the field shape function FSF to the infinite element IFE;

[0533] (f3) using the field shape function FSF in the numerical simulation model MDL to realize infinite elements;

[0534] (g) determining the field variables FDV by solving the numerical simulation model MDL using the solver SLV;

[0535] (h) displaying the field variable FDV and / or key figures calculated based on the field variable FDV or an image illustrating the field variable FDV to the user USR via the display DSP;

[0536] (i) providing the field variable FDV to an iterative object design process 10D of the object OBJ for designing the object OBJ with the improvement of acoustic emission as a design goal;

[0537] (j) generating said object OBJ according to the design as determined by said iterative object design process IOD. This generating GNT can be done by any manufacturing or production method.

[0538] The method can be implemented as a computer product CMP or stored thereon or as it can be downloaded, the computer product CMP being arranged and configured to perform the steps of the computer-implemented method according to the invention or one of its embodiments. The computer product CMP can be computer hardware comprising executable instructions stored in a data store. Another additional or alternative option is to provide the method according to the invention via a computer-readable medium CRM encoded with executable instructions, which, when executed, causes a computer or computer product to implement the method according to the invention or one of its embodiments. Another option is to provide the invention as a system SYS for determining acoustic parameters of an emission of an object, the system SYS comprising at least one processor CPU, which is prepared to execute the method according to the invention or one of its embodiments by uploading a computer executable code.

[0539] While the present invention has been described in detail with reference to preferred embodiments, it will be understood that the invention is not limited to the disclosed examples and that those skilled in the art can make numerous additional modifications and variations thereto without departing from the scope of the present invention.

[0540] It should be noted that "a" or "an" used throughout this application does not exclude a plurality, and "comprising" does not exclude other steps or elements. In addition, elements described in association with different embodiments can be combined. It should also be noted that the reference numerals of the present invention should not be interpreted as limiting the scope of the present invention.

[0541] Independent of grammatical usage of the term, individuals with masculine, feminine, or other gender identities are included in the term.

Claims

1. A computer-implemented method for modeling the acoustics of an object (OBJ) in an unbounded domain (UBD) as a numerical simulation model (MDL) to determine field variables (FDV) within a control volume (CVL) of the unbounded domain, the method comprising: (a) providing a coordinate system (COS) mapping the position (POS) in said control volume (CVL), (b) providing a geometrical specification of said object (OBJ), said object (OBJ) acting as a sound source (SDS) emitting sound in said control volume (CVL), wherein said geometrical specification is based on said coordinate system (COS), Characterized in that the method further comprises the additional steps of: (c) constructing a convex shape (CXS) around the sound source (SDS) using a fast convex hull algorithm (QHA), (d) stretching an infinite element (IFE) starting from the convex shape (CXS), (e) by mapping the radial coordinates of the coordinate system (COS) to auxiliary coordinates extending from the convex shape (CXS) to the infinite element (IFE) Transform the coordinate system (COS) into the field variable coordinate system (FVC), (f1) Using the auxiliary coordinates To construct the field shape function (FSF), (f2) applying the field shape function (FSF) to the infinite element (IFE); (f3) implementing the infinite element using the field shape function (FSF) in the numerical simulation model (MDL); (g) determining the field variables (FDV) by solving the numerical simulation model (MDL).

2. The method according to claim 1, wherein: Step (e) converts the auxiliary coordinates of the position outside the convex shape (CXS) Provided as: in: s, t, The three-dimensional coordinates (s, t are tangential coordinates, is the radial coordinate); Identify the point on the boundary Γ of the convex shape (CXS) that is closest to the external position, is the distance between the center of the sound source (SDS) and a specific point on Γ; defining a distance between the sound source (SDS) and the location; is the auxiliary coordinate.

3. The method according to any one of the preceding claims 1 or 2, comprising: (h) Displaying to a user (USR) the field variables (FDV) and / or key figures calculated from the field variables (FDV) or an image illustrating the field variables (FDV).

4. The method according to any one of the preceding claims 1, 2 or 3, comprising: (i) providing the field variables (FDV) to an iterative object design process (IOD) of the object (OBJ) for designing the object (OBJ) with the improvement of acoustic emission as a design goal.

5. A computer product (CMP) arranged and configured to perform the steps of the computer-implemented method according to any of the preceding claims 1 to 4.

6. A computer readable medium (CRM) encoded with executable instructions which, when executed, cause a computer product according to claim 5 to implement a method according to any one of claims 1 to 4.

7. A system (SYS) for determining acoustic parameters of an emission of an object, the system (SYS) comprising at least one processor (CPU) which is prepared to execute the method according to any of the preceding claims 1 to 4 by uploading a computer executable code.

Images