JP2022112624A - Sound field analysis apparatus, sound field analysis method, and program - Google Patents

Abstract

To efficiently handle an impedance boundary depending on a frequency in sound field analysis using a time domain finite element method.SOLUTION: A sound field analysis device 10 analyzes a sound field by using a time domain finite element method based on a first order normal differential equation relating to a sound pressure vector showing a distribution of a sound pressure in a space to be analyzed and a differential vector equivalent to a time first derivative of the sound pressure vector, the first order normal differential equation including an attenuation term for evaluating attenuation of sound waves at a boundary surrounding the space. A numerical value calculation section 120 calculates a value of the sound pressure vector at a plurality of time steps according to a time update equation of the sound pressure vector obtained by discretizing the time first derivative of the sound pressure vector by a difference approximation in the first order normal differential equation, and a time update equation of the differential vector obtained by approximating convolution of an inverse Fourier transform value of a specific acoustic admittance ratio depending on a frequency included in the attenuation term and the time first derivative of the sound pressure vector by an Auxiliary Differential Equations method.SELECTED DRAWING: Figure 1

Description

The present invention relates to a sound field analysis device, sound field analysis method and program.

A method of analyzing a sound field by numerical simulation is known. For example, Non-Patent Documents 1 and 2 disclose a method of analyzing a sound field using the Time-Domain Finite Element Method (TDFEM) based on the Ordinary Differential Equation (ODE). is doing. Specifically, Non-Patent Document 1 discloses a method of analyzing a sound field according to the following (Procedure 1) to (Procedure 4).

(Step 1)
It is obtained by introducing the diagonal mass matrix and the vector equivalent to the time-first derivative of the sound pressure into the semi-discretized equation for the sound field surrounded by the rigid wall boundary, the vibration boundary, or the impedance boundary. We are interested in a system of first-order ordinary differential equations Specifically, the simultaneous first-order ordinary differential equations are a vector p (hereinafter referred to as "sound pressure vector p") indicating the distribution of sound pressure in the space to be analyzed, and the time first-order differential of the sound pressure vector p A vector v equivalent to (hereinafter referred to as a "differential vector v") is expressed as the following equations (1) and (2) as unknowns.

Here, "M", "K", "C", and "D" in formulas (1) and (2) represent mass matrix, stiffness matrix, damping matrix, and diagonal mass matrix, respectively. The mass matrix M, stiffness matrix K, damping matrix C, and diagonal mass matrix D superimpose the element mass matrix M _e , the element stiffness matrix K _e , the element damping matrix C _e , and the element diagonal mass matrix D _e , respectively. It is the overall matrix constructed together. Further, in the equations (1) and (2), “f” represents an external force vector indicating the distribution of the external force applied to the space to be analyzed, “c ₀ ” represents the speed of sound, and “·” represents time-related represents the derivative. Also, y _n represents the frequency-independent specific acoustic admittance ratio. y _n is the reciprocal of the specific acoustic impedance ratio.

(Step 2)
A square or cubic linear linear element is set as the finite element of the air region. Then, the modified integral rule, which is a discretization error reduction method disclosed in Non-Patent Document 2, is applied to the set finite elements of the air region, and the mass matrix in equations (1) and (2) Compute M and the stiffness matrix K.

(Step 3)
A time progression scheme is constructed by discretizing the simultaneous first-order ordinary differential equations represented by equations (1) and (2) in the time direction, and the time progression of sound waves is calculated.

(Step 4)
In order to maintain numerical stability when calculating the time progression of sound waves in (procedure 3), the time first derivative of the sound pressure vector p in equation (1) is discretized by forward difference approximation with first-order accuracy, and (2) Discretize the first-order time differential of the differential vector v of the equation (2) by backward difference approximation with first-order accuracy. As a result, the following formulas (3) and (4) are obtained.

Here, "n" in equations (3) and (4) is a natural number representing the number of time steps, and "Δt" represents the time step width. In the method disclosed in Non-Patent Document 1, the progress of sound waves over time is calculated according to the time update formulas represented by formulas (3) and (4). Although equation (4) is an implicit equation with respect to ^vn , it can be transformed into an explicit equation with respect to ^vn by centralizing the attenuation matrix C.

The time-domain finite element method based on such a first-order ordinary differential equation has the following two advantages regarding computational efficiency.
(I) Since the time marching scheme is explicit, no simultaneous linear equation solution is required. Therefore, the computational load per time step is small.
(II) By application of the modified integral rule, fourth order accuracy in time and space can be achieved despite the fact that the elements used are first order elements.

Takumi Yoshida, Takeshi Okuzono, Kimihiro Sakagami, “Numerically stable explicit time-domain finite element method for room acoustics simulation using an equivalent impedance model”, Noise Control Engineering Journal, 66(3), 176-189 (1 May2018) B. Yue, MN. Guddati, “Dispersion-reducing finite elements for transient acoustics”, Journal of Acoustical Society of America, 118(4), 2132-2141(2005)

In acoustic simulations using computational mechanics methods including explicit time-domain finite element methods, how to handle sound absorption at the interface greatly influences the prediction results. In particular, proper handling of frequency dependence in boundary sound absorption is important for simulating actual acoustic phenomena. However, the acoustic analysis method based on the explicit time-domain finite element method considers the sound absorption at the boundary with a frequency-independent surface impedance model because of its ease of implementation and low computational load. is treated as an approximation. This is because the formulation of a highly accurate sound absorption model in time domain analysis involves convolution with a large computational load due to consideration of frequency dependent characteristics. Therefore, in sound field analysis using the time-domain finite element method, a method that can efficiently handle frequency-dependent impedance boundaries is required.

SUMMARY OF THE INVENTION The present invention is intended to solve the above problems. An object of the present invention is to provide a field analysis method and program.

In order to achieve the above object, the sound field analysis device according to the first aspect of the present invention includes:
A first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in a space to be analyzed and a differential vector equivalent to the time first order differential of said sound pressure vector, wherein the sound wave at the boundary surrounding said space A sound field analysis device for analyzing a sound field using a time domain finite element method based on the first-order ordinary differential equation including an attenuation term for evaluating the attenuation of
In the first-order ordinary differential equation, the time update formula of the sound pressure vector obtained by discretizing the time first-order differential of the sound pressure vector by difference approximation, and in the first-order ordinary differential equation, the attenuation term A time update formula for the differential vector obtained by approximating the convolution of the included frequency-dependent inverse Fourier transform value of the specific acoustic admittance ratio with the time first derivative of the sound pressure vector using the Auxiliary Differential Equations method, according to , a numerical calculator for calculating values of the sound pressure vector at a plurality of time steps.

The differential vector time update formula includes an auxiliary function obtained by rational function approximation of the specific acoustic admittance ratio,
The numerical calculation unit may calculate the time evolution of the auxiliary function by a time integration method.

The numerical calculation unit may calculate the time evolution of the auxiliary function using an explicit Runge-Kutta method with fourth-order precision as the time integration method.

The time update formula for the sound pressure vector is obtained by combining the value of the sound pressure vector at the nth time step with the value of the sound pressure vector at the (n−1)th time step and the (n−1)th time step. and the value of the differential vector at
The time update formula for the differential vector is such that the value of the differential vector at the nth time step is combined with the value of the sound pressure vector at the (n−1)th time step and the nth time step, and the (n−1)th ) A formula determined by the value of the differential vector at the time step and the value of the auxiliary function at the n-th time step.

The time update formula for the differential vector discretizes the time first derivative of the differential vector by backward difference approximation with first-order accuracy, and further discretizes the time first derivative of the sound pressure vector by central difference approximation with second-order accuracy. It may be obtained by

In order to achieve the above object, a sound field analysis method according to a second aspect of the present invention comprises:
A first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in a space to be analyzed and a differential vector equivalent to the time first order differential of said sound pressure vector, wherein the sound wave at the boundary surrounding said space A sound field analysis method for analyzing a sound field using a time domain finite element method based on the first-order ordinary differential equation including an attenuation term for evaluating the attenuation of
In the first-order ordinary differential equation, the time update formula of the sound pressure vector obtained by discretizing the time first-order differential of the sound pressure vector by difference approximation, and in the first-order ordinary differential equation, the attenuation term A time update formula for the differential vector obtained by approximating the convolution of the included frequency-dependent inverse Fourier transform value of the specific acoustic admittance ratio with the time first derivative of the sound pressure vector using the Auxiliary Differential Equations method, according to , to calculate the values of the sound pressure vector at multiple time steps.

In order to achieve the above object, a program according to the third aspect of the present invention,
the computer,
A first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in a space to be analyzed and a differential vector equivalent to the time first order differential of said sound pressure vector, wherein the sound wave at the boundary surrounding said space A program that functions as a sound field analysis device that analyzes a sound field using a time domain finite element method based on the first-order ordinary differential equation including an attenuation term for evaluating the attenuation of
said computer,
In the first-order ordinary differential equation, the time update formula of the sound pressure vector obtained by discretizing the time first-order differential of the sound pressure vector by difference approximation, and in the first-order ordinary differential equation, the attenuation term A time update formula for the differential vector obtained by approximating the convolution of the included frequency-dependent inverse Fourier transform value of the specific acoustic admittance ratio with the time first derivative of the sound pressure vector using the Auxiliary Differential Equations method, according to , a numerical calculation unit that calculates the values of the sound pressure vector at a plurality of time steps.

According to the present invention, frequency-dependent impedance boundaries can be efficiently handled in sound field analysis using the time domain finite element method.

1 is a block diagram showing the configuration of a sound field analysis device according to an embodiment of the present invention; FIG. It is a figure which shows the setting example of the mesh in embodiment. 4 is a block diagram showing the configuration of a numerical calculation unit in the sound field analysis device according to the embodiment; FIG. FIG. 10 is a diagram showing a display example of temporal changes in sound pressure obtained by the sound field analysis device according to the embodiment; FIG. 4 is a diagram showing a display example of the spatial distribution of sound pressure obtained by the sound field analysis device according to the embodiment; 4 is a flow chart showing the flow of processing executed by the sound field analysis device according to the embodiment; FIG. 4 is a diagram showing a cross section of an acoustic tube model used as a space to be analyzed in numerical experiments; FIG. 4 is a diagram showing a comparison between the sound absorption coefficient of a porous sound absorber calculated by a rational function model and a theoretical value in numerical experiments. FIG. 5 is a diagram showing a comparison between theoretical values and sound absorption coefficients of porous sound absorbers calculated by the method of the embodiment and the conventional method in numerical experiments.

BEST MODE FOR CARRYING OUT THE INVENTION Hereinafter, embodiments of the present invention will be described with reference to the drawings. The same reference numerals are given to the same or corresponding parts in the drawings.

(embodiment)
The sound field analysis device 10 according to the present embodiment analyzes the sound field by wave acoustic numerical simulation using the time domain finite element method, and can be used for the development of materials and measurement methods in the acoustic field, design support, etc. is a device capable of Here, the sound field means a space in which sound waves exist.

FIG. 1 shows the configuration of the sound field analysis device 10. As shown in FIG. As shown in FIG. 1 , the sound field analysis device 10 includes a control section 11 , a storage section 12 , an operation section 13 , a display section 14 and a communication section 15 .

The control unit 11 is connected to each unit of the sound field analysis device 10 via a control bus, and controls the overall operation of the sound field analysis device 10 . Specifically, the control unit 11 includes a CPU (Central Processing Unit), a ROM (Read Only Memory), and a RAM (Random Access Memory). A CPU is, for example, a microprocessor or the like, and is a central processing unit that performs various processes and calculations. In the control unit 11, the CPU reads the control program stored in the ROM and executes various processes while using the RAM as a work memory.

The storage unit 12 is a non-volatile memory such as flash memory or hard disk. The storage unit 12 stores programs and data used by the control unit 11 to perform various processes, such as an operating system and application programs. The storage unit 12 also stores data generated or acquired by the control unit 11 performing various processes.

The operation unit 13 includes input devices such as a keyboard, mouse, buttons, touch pad, and touch panel, and receives operations from the user. A user can input various instructions to the sound field analysis device 10 by operating the operation unit 13 . Upon receiving an operation instruction input by the user, the operation unit 13 transmits the received operation instruction to the control unit 11 .

The display unit 14 includes a display device such as a liquid crystal display or an organic EL (Electro Luminescence) display. The display unit 14 is driven by a display driving circuit (not shown) and displays various images under the control of the control unit 11 . For example, the display unit 14 displays an image representing the result of analyzing the sound field by numerical simulation.

The communication unit 15 has a communication interface for the sound field analysis device 10 to communicate with an external device. The communication unit 15 communicates with external devices through wired or wireless communication such as a wide area network such as the Internet, a LAN (Local Area Network), or a USB (Universal Serial Bus).

The control unit 11 functionally includes a setting reception unit 110, a numerical calculation unit 120, and an output unit 130, as shown in FIG. In the control unit 11, the CPU functions as each of these units by reading the program stored in the ROM into the RAM and executing it.

The setting reception unit 110 receives settings for analyzing the sound field in the sound field analysis device 10 . Here, the settings are conditions, parameters, and the like for the numerical calculation unit 120 to calculate the distribution of sound pressure by numerical simulation, and specifically, information regarding the space to be analyzed, boundary conditions, meshes, and the like. The setting reception unit 110 displays a setting screen on the display unit 14 for the user to input desired settings. The user operates the operation unit 13 to input settings while viewing the setting screen displayed on the display unit 14 . The setting reception unit 110 receives settings input by the user.

First, the setting reception unit 110 receives the setting of the space to be analyzed, which is the object of the sound field analysis. Specifically, the setting reception unit 110 receives settings for the dimension and shape of the space to be analyzed according to an operation input by the user via the operation unit 13 .

For example, when analyzing a sound field in a rectangular two-dimensional plane, the setting reception unit 110 receives the setting of a rectangular two-dimensional space as the space to be analyzed. Alternatively, when analyzing a sound field in a three-dimensional indoor space, the setting reception unit 110 receives a setting of a three-dimensional space having a shape corresponding to the indoor space as a space to be analyzed.

Second, the setting receiving unit 110 receives setting of boundary conditions. A boundary condition is a condition regarding a boundary surrounding a space to be analyzed. The setting reception unit 110 receives a setting of one of a rigid wall boundary, a vibration boundary, and an impedance boundary as a boundary condition according to an operation input by the user via the operation unit 13 .

A rigid wall boundary is a boundary that is acoustically rigid and the velocity of particles on the boundary is zero. A sound wave incident on a rigid wall boundary is completely reflected without being absorbed. In contrast, an oscillating boundary is a boundary on which particles can oscillate. Also, an impedance boundary is a boundary having an acoustic impedance. A sound wave incident on an impedance boundary is absorbed with an absorption rate corresponding to the acoustic impedance value _zn of the boundary. Therefore, sound waves are attenuated by reflection at impedance boundaries.

When an impedance boundary is included in at least a part of the boundaries surrounding the space to be analyzed, the first-order ordinary differential equation includes the following equation (2) in order to evaluate the attenuation of the sound wave reflected by the impedance boundary. A damping term containing the damping matrix C is required as follows. On the other hand, if the boundary surrounding the space to be analyzed does not include an impedance boundary, the attenuation term including the attenuation matrix C in equation (2) becomes unnecessary. In the following, an example will be described in which an impedance boundary is included in at least a part of the boundary surrounding the space to be analyzed. When setting an impedance boundary as a boundary condition, the setting reception unit 110 receives a setting of the acoustic impedance value _zn at the impedance boundary from the user.

Third, the setting reception unit 110 receives mesh settings. Here, the mesh is obtained by dividing (discretizing) the space to be analyzed into a plurality of finite elements having simple shapes such as quadrilaterals in order to numerically analyze the space to be analyzed. Note that finite elements are also simply referred to as “elements”. The setting accepting unit 110 accepts settings for the shape and size of a plurality of elements forming the mesh according to an operation input by the user via the operating unit 13 .

FIG. 2 shows an example in which a two-dimensional analysis target space 20 is divided by rectangular elements. When the space 20 to be analyzed is two-dimensional, for example, as shown in FIG. divided into pieces. That is, the space 20 to be analyzed is discretized with a total of (X×Y) elements. Alternatively, when the space to be analyzed is three-dimensional, although illustration is omitted, the space to be analyzed is divided into the horizontal direction, the vertical direction, and the height direction by three-dimensional elements such as cubes and rectangular parallelepipeds. be done.

The sound field analysis apparatus 10 analyzes the distribution of sound pressure in the space to be analyzed by calculating the value of the sound pressure in units of the discretized elements. As the size of the elements is set smaller, the error due to the discretization of the space becomes smaller, so that the analysis accuracy improves, but the amount of calculation increases.

Fourth, the setting reception unit 110 receives settings of the time step size Δt, the total number of steps Ns, and the external force vector f. The time step size Δt is the unit of time step for calculating the time progression of sound waves in the space to be analyzed. The smaller the time step width Δt is set, the smaller the error associated with the discretization of time, which improves the accuracy of analysis, but increases the amount of calculation. The total number of steps Ns is a value indicating how many steps forward the time progression of sound waves in the space to be analyzed is to be calculated from the initial state. The numerical calculation unit 120 calculates the sound pressure over a time interval of Δt×Ns. The external force vector f is a vector representing the distribution of the external force applied to the space to be analyzed in the time interval from the start of the calculation to the end of the calculation.

In this manner, the setting reception unit 110 receives settings for analyzing the sound field by the time domain finite element method from the user via the operation unit 13 . The setting reception unit 110 is implemented by the control unit 11 cooperating with the operation unit 13 and the like.

Returning to FIG. 1, numerical calculation section 120 performs numerical calculation under the settings accepted by setting acceptance section 110 . More specifically, the numerical calculation unit 120 calculates a first-order constant for a sound pressure vector p indicating the distribution of sound pressure in the space to be analyzed and a differential vector v equivalent to the first-order time differential of the sound pressure vector p. The values of the sound pressure vector p at a plurality of time steps are calculated according to the time update formula of the sound pressure vector p and the time update formula of the differential vector v obtained from the differential equation. The first-order ordinary differential equations are shown in equations (1) and (2) and include attenuation terms to estimate the attenuation of sound waves at boundaries.

Here, a calculation method by the numerical calculation unit 120 will be described. In the present embodiment, the numerical calculation unit 120 uses Equation (3) described in Background Art as the time update equation for the sound pressure vector p. Equation (3), as described in the background art, is a time update of the sound pressure vector p obtained by discretizing the time first order differential of the sound pressure vector p in the equation (1) by the difference approximation of the first order accuracy. is the formula. On the other hand, the numerical calculation unit 120 uses a time update formula different from the formula (4) described in the background art as a time update formula for the differential vector v in order to handle the frequency dependence of sound absorption at the boundary.

Derivation of the time update formula for the differential vector v in this embodiment will be described below. When dealing with frequency-dependent impedance boundaries, equation (2) described in Background Art is rewritten as equation (5) below. Here, "*" in the equation (5) represents convolution defined by the following equation (6). Also, in the equation (6), y with "^" represents the inverse Fourier transform value of the frequency-dependent specific acoustic admittance ratio y(ω).

Here, in general, acoustic analysis by the time domain finite element method requires a large number of calculation time steps. Therefore, since the computational cost of convolution in equation (5) is enormous, it is difficult to practically use a method involving convolution. Therefore, efficient convolution techniques are needed to deal with frequency-dependent impedance boundaries in the time-domain finite element method.

Therefore, in this embodiment, the ADE (Auxiliary Differential Equations) method, which is a recursive convolution method with excellent computational efficiency, is used to efficiently incorporate frequency dependence into the time domain finite element method. Details of the ADE method are disclosed, for example, in "D. Dragna, A generalized recursive convolution method for time-domain propagation in porous media, Journal of Acoustical Society of America, 1030-1042 (2015)".

When using the ADE method, the frequency-dependent specific acoustic admittance ratio y(ω) is approximated by a rational function as shown in Equation (7) below. (7), ω and j are the angular frequency and the imaginary unit j ² =−1, respectively. Also, n _rp , n _cp are the number of real poles λ _i and the number of complex conjugate pole pairs α _i ±jβ _i used for rational function approximation, respectively. y _∞ , A _i , B _i , C _i are coefficients.

The poles λ _i , α _i , β _i and the coefficients y _∞ , A _i , B _i , C _i in the equation (7) are all real numbers, and their values are obtained by using a data fitting method. For example, a calculated value obtained by calculating y(ω) at a plurality of frequencies using equation (7) in which appropriate initial values are given to the poles and coefficients described above is compared with the theoretical or experimental value. Then, using a data fitting technique, poles λ _i , α _i , β _i and coefficients y _∞ , A Determine _i , B _i , C _i . Note that λ _i and α _i are positive values in order to satisfy causality.

By using the rational function-approximated specific acoustic admittance ratio y(ω), the convolution of the equation (5) becomes the auxiliary functions φ _i , ψ _i ⁽¹⁾ , ψ _i ^{( 2)} can be used for approximation. The auxiliary functions φ _i , φ _i ⁽¹⁾ , φ _i ⁽²⁾ are called accumulators.

Here, the time evolution of the auxiliary functions φ _i , ψ _i ⁽¹⁾ and ψ _i ⁽²⁾ in the equation (8) are represented by the following equations (9), (10) and (11) respectively. It is calculated by an ordinary differential equation.

The numerical calculation unit 120 calculates the time evolution of the auxiliary functions φ _i , φ _i ⁽¹⁾ and φ _i ⁽²⁾ by the time integration method. Well-known methods such as the Runge-Kutta method, the Euler method, the Crank-Nicholson method, and the Huyn method can be used as the time integration method.

The accuracy and stability of the time integration method affect the accuracy and stability of the entire calculation. Therefore, it is necessary to choose an appropriate time integration method for efficient simulation. As shown in the numerical experiments described later, in the explicit time-domain finite element method, the explicit Runge-Kutta method with fourth-order accuracy can be used to calculate the equation (8) with high accuracy and stability. Therefore, the case where the explicit Runge-Kutta method with fourth-order precision is used as the time integration method will be described below as an example. For details of the fourth-order precision explicit Runge-Kutta method in the time-domain finite element method, see, for example, "P.O.J. Scherer, Computational Physics: Simulation of Classical and Quantum Systems, Third ed. (Springer Nature, Berlin/Heidelberg, 2017)." disclosed in

We then construct an explicit scheme for dealing with frequency-dependent impedance boundaries. Specifically, the following equation (12) is derived by substituting equation (8) into equation (5).

In order to stably solve the equation (12), the time first derivative of the differential vector v is discretized by the backward difference approximation of the first order accuracy. Furthermore, in order to calculate the frequency-dependent impedance boundary with high accuracy, the first-order time derivative of the sound pressure vector p is discretized by a second-order accuracy central difference approximation. At that time, the first-order time differential of the sound pressure vector p is discretized as shown in the following equation (13) based on the equation (3).

Furthermore, for a fully explicit scheme, the mass matrix M in Eq. (13) is centralized and transformed into a diagonal mass matrix D. As a result, the time update formula for the differential vector v is expressed as the following formula (14). Equation (14) is a differential vector v is the time update formula of

Numerical calculation unit 120 calculates a plurality of time steps according to equation (14), which is the time update equation for differential vector v and equation (3), which is the time update equation for sound pressure vector p. Calculate the value of the sound pressure vector p at

FIG. 3 shows a more detailed configuration of the numerical calculation unit 120. As shown in FIG. As shown in FIG. 3, the numerical calculation unit 120 has the functions of a matrix setting unit 121, a sound pressure vector calculation unit 123, an auxiliary function calculation unit 124, a differential vector calculation unit 125, and a repetition unit 127. .

The matrix setting unit 121 sets the mass matrix M, the stiffness matrix K, the damping matrix C, and the diagonal mass matrix D included in the equations (3) and (14). The mass matrix M, stiffness matrix K, damping matrix C, and diagonal mass matrix D are matrices of sizes corresponding to the number of elements that discretize the space to be analyzed, and the values of the elements included in each matrix are , is determined depending on the dimension and shape of the space to be analyzed and the shape of the elements to be discretized. Therefore, the matrix setting unit 121 sets the mass matrix M, the stiffness matrix K, the damping matrix C, and the diagonal mass matrix D according to the settings received from the user by the setting receiving unit 110 . Thereby, the matrix setting unit 121 establishes a time update formula for the sound pressure vector p and the differential vector v.

More specifically, the matrix setting unit 121 applies the modified integral rule disclosed in Non-Patent Document 2 to calculate the mass matrix M and the stiffness matrix K. As a calculation method, for example, following Non-Patent Document 2, the matrix setting unit 121 calculates the Courant number τ from the time step width Δt and the element length h set by the setting reception unit 110, and the numerical integration point α _m and Set α _k . Then, the matrix setting unit 121 calculates the element mass matrix M _e and the element stiffness matrix K _e by the Gauss-Legendre quadrature method using the set numerical integration points α _m and α _k , and calculates the mass matrix M and the stiffness matrix K set.

Here, the method of calculating the element mass matrix M _e and the element stiffness matrix K _e from the numerical integration points α _m and α _k and further obtaining the mass matrix M and the stiffness matrix K are disclosed in Non-Patent Documents 1 and 2, etc. well-known techniques can be used. In addition, the matrix setting unit ₁₂₁ calculates the element diagonal mass matrix De, which is a diagonal matrix, using a well-known method disclosed in Non-Patent Documents 1, 2, etc., and sets the diagonal mass matrix D do.

Furthermore, the matrix setting unit 121 sets the attenuation matrix C by calculating the element attenuation matrix C _e using the acoustic impedance values _zn set at the impedance boundaries. As a method of calculating the element attenuation matrix C _e from the acoustic impedance values z _n and obtaining the attenuation matrix C, a well-known method disclosed in Non-Patent Document 1 or the like can be used.

The sound pressure vector calculator 123 calculates the sound pressure vector ^pn at the n-th time step according to formula (3), which is a time update formula for the sound pressure vector p. Here, the sound pressure vector p ⁿ at the n-th time step has as elements the value of the sound pressure at the n-th time step from the start of calculation in each of the plurality of elements set in the space to be analyzed. is a vector.

As described above, the equation (3) ^expresses the value of the sound pressure vector p_n at the n-th time step as the value of the sound pressure vector p_n-1 at the (n-1)th time step and the value of the sound pressure vector p_n- ¹ at the (n-1)th time step. -1) the value of the differential vector v ^n-1 of the time step; Therefore, the sound pressure vector calculator 123 calculates the value of the sound pressure vector p n at the n-th time step, the value of the sound pressure vector p ^{n -1} at the (n−1) time step, and the (n−1)-th time step. 1) Calculate from the value of the differential vector v ⁿ⁻¹ of the time step.

The auxiliary function calculation unit 124 calculates the auxiliary functions φ _i , ψ _i ⁽¹⁾ , φ _i ⁽ 1) , φ i (1) , according to the time update formulas (9) to (11) of the auxiliary functions φ _i , φ _i ⁽¹⁾ , φ i (2). Compute the time evolution of ψ _i ⁽²⁾ . As described above, Equation (14), which is a time update equation for the differential vector v, is an auxiliary function φ _i , ψ _i ⁽¹⁾ , ψ _i obtained by rational function approximation of the specific acoustic admittance ratio y(ω). ⁽²⁾ is included. The auxiliary function calculator 124 calculates the time evolution of the auxiliary functions φ _i , ψ _i ⁽¹⁾ , ψ _i ⁽²⁾ using the fourth-order precision explicit Runge-Kutta method.

Specifically, the auxiliary function calculator 124 calculates the values of the auxiliary functions φ _i , ψ _i ⁽¹⁾ , ψ _i ⁽²⁾ at the (n−1)th time step, By solving the ordinary differential equations (9) to (11) using the first time derivative of the sound pressure vector p ⁿ⁻¹ at the step and the explicit Runge-Kutta method with fourth-order accuracy, the n-th time step Compute the values of the auxiliary functions φ _i , φ _i ⁽¹⁾ , φ _i ⁽²⁾ at . The differential vector v ^n-1 at the (n-1)th time step is used as the time first-order differential of the sound pressure vector p ^n-1 at the (n-1)th time step.

The differential vector calculator 125 calculates the differential vector v ⁿ at the n-th time step according to Equation (14), which is the time update formula for the differential vector v. Here, the differential vector v ⁿ at the n-th time step is the time first-order differential of the sound pressure at the n-th time step from the calculation start time step in each of the plurality of elements set in the space to be analyzed. A vector whose elements are the corresponding values.

As described above, the equation (14) expresses the value of the differential vector v ⁿ at the n-th time step as the sound pressure vectors p ^n-1 and p ⁿ at the (n-1)-th time step and the n-th time step. , the value of the differential vector v ^n-1 at the (n-1)th time step, the values of the auxiliary functions φ _i , ψ _i ⁽¹⁾ , ψ _i ⁽²⁾ at the n-th time step, It is a formula determined by Therefore, the differential vector calculator 125 calculates the value of the differential vector v ⁿ at the n-th time step, the values of the sound pressure vectors p ⁿ⁻¹ and p ⁿ calculated by the sound pressure vector calculator 123, and the auxiliary function calculation. It is calculated from the values of the auxiliary functions φ _i , φ _i ⁽¹⁾ and φ _i ⁽²⁾ calculated by the unit 124 and the value of the differential vector v ⁿ⁻¹ at the (n−1)th time step.

The repetition unit 127 performs processing for calculating the value of the sound pressure vector p ⁿ at the n-th time step by the sound pressure vector calculation unit 123 while incrementing the value of n by 1, A process of calculating the values of the auxiliary functions φ _i , ψ _i ⁽¹⁾ and ψ _i ⁽²⁾ at each step, and a process of calculating the value of the differential vector v ⁿ at the n-th time step by the differential vector calculator 125 . ,repeat. More specifically, the repeating unit 127 changes the value of n from 1 to the total number of steps Ns, and calculates the sound pressure vector p ^k and the auxiliary functions φ _i , ψ _i ⁽¹⁾ when n=k. , ψ _i ⁽²⁾ and the differential vector v ^k are the sound pressure vector p ^k−1 calculated when n=k−1, the auxiliary functions φ _i , ψ _i ⁽¹⁾ , ψ _i ⁽²⁾ and the value of the differential vector v ^k−1 . Thereby, the repeater 127 calculates the values of the sound pressure vectors p ¹ to p ^Ns and the differential vectors v ¹ to v ^Ns at a plurality of time steps from the first time step to the Ns-th time step.

Note that the values of the sound pressure vector p ⁿ and differential vector v ⁿ when n=0 and the auxiliary functions φ _i , φ _i ⁽¹⁾ and φ _i ⁽²⁾ are not defined. Therefore, the values of these vectors are preset to appropriate values such as 0 as initial conditions. Then, the repeater 127 calculates the sound pressure vector ^pn and the differential vector ^vn according to the equations (3) and (14).

In this way, the numerical calculation unit 120 calculates, step by step, the sound pressure vector p representing the distribution of the sound pressure in the space to be analyzed and the differential vector v equivalent to the time first derivative of the sound pressure vector p. Analyze time progression. Numerical calculation unit 120 is realized by cooperation of control unit 11 with storage unit 12 and the like.

Returning to FIG. 1 , the output unit 130 outputs output information based on the sound pressure vectors p at multiple time steps calculated by the numerical calculation unit 120 . The output information based on the sound pressure vector p at a plurality of time steps is based on the sound pressure vectors p ¹ to p ^Ns from the first time step to the Ns-th time step obtained by numerical calculation by the numerical calculation unit 120. It is information that can be expressed. The output unit 130 outputs, as output information, information representing at least one of a temporal change in sound pressure and a spatial distribution of sound pressure.

For example, as shown in FIG. 4, the output unit 130 displays, on the display unit 14, information indicating temporal changes in sound pressure at a specific position in the space to be analyzed as output information. The output unit 130 outputs the position specified by the user among the sound pressure values included in the sound pressure vectors p ¹ to p ^Ns from the first time step to the Ns-th time step calculated by the numerical calculation unit 120. Arrange sound pressure values in chronological order. As a result, the output unit 130 generates an image representing the temporal change in sound pressure at a specific position, and outputs the image to the display unit 14 as shown in FIG.

As another example, as shown in FIG. 5, the output unit 130 displays, as output information, information indicating the spatial distribution of sound pressure in the space to be analyzed at a specific time step on the display unit 14 . The output unit 130 outputs the time step specified by the user among the sound pressure values included in the sound pressure vectors p ¹ to p ^Ns from the first time step to the Ns-th time step calculated by the numerical calculation unit 120. Arrange the values of the sound pressure vector in to express the spatial distribution. As a result, the output unit 130 generates an image representing the spatial distribution of sound pressure at a specific time step, and outputs the image to the display unit 14 as shown in FIG. In FIG. 5, an image representing the distribution in one-dimensional directions (vertical direction, horizontal direction, etc.) is displayed as an example. may be displayed.

Alternatively, the output unit 130 may output both the temporal change and the spatial distribution of the sound pressure. For example, the output unit 130 generates a moving image representing temporal changes in the spatial distribution of sound pressure at a plurality of time steps in order to visualize how sound waves propagate in the space to be analyzed, and the display unit 14 You can display output to . By confirming the temporal change and spatial distribution of the sound pressure displayed on the display unit 14 in this way, the user can obtain information on how sound waves propagate or are distributed in the space to be analyzed. .

Furthermore, the output unit 130 may output, as output information, an acoustic parameter (for example, reverberation time) representing the acoustic characteristics of the space to be analyzed, in addition to the temporal change and the spatial distribution. The user can appropriately set through the operation unit 13 what kind of information is to be output from the output unit 130 as the output information. The output unit 130 is implemented by the control unit 11 cooperating with the display unit 14 and the like.

The flow of sound field analysis processing executed by the sound field analysis apparatus 10 configured as described above will be described with reference to the flowchart shown in FIG. The sound field analysis process shown in FIG. 6 is appropriately started according to an instruction input from the user in a state where the sound field analysis apparatus 10 can execute the sound field analysis process.

When the sound field analysis process is started, the control unit 11 functions as the setting receiving unit 110 and receives setting of conditions and parameters for the sound field analysis (step S11). Specifically, the control unit 11 controls the dimension and shape of the space to be analyzed, boundary conditions, element shape and size, time step width Δt, total number of steps Ns, external force vector f, acoustic impedance value z at the impedance boundary Settings such as _n are received from the user via the operation unit 13 .

Upon receiving the settings, the control unit 11 functions as the matrix setting unit 121 and sets the mass matrix M, the stiffness matrix K, the damping matrix C, and the diagonal mass matrix D (step S12). Specifically, the control unit 11 sets numerical integration points α _m and α _k from the time step width Δt set in step S11 and the element length h of the mesh. Then, the control unit 11 calculates the element mass matrix M _e and the element stiffness matrix K _e from the numerical integration points α _m and α _k , respectively, and sets the mass matrix M and the stiffness matrix K. Furthermore, the control unit 11 calculates the element attenuation matrix C _e from the acoustic impedance value z _n set in step S11, and sets the attenuation matrix C.

After setting the matrix, the control unit 11 initializes the value of n to 1 (step S13). Then, the control unit 11 shifts to the process of calculating the sound pressure vector p and the differential vector v step by step while increasing the value of n.

First, the control section 11 functions as the sound pressure vector calculation section 123 and calculates the sound pressure vector ^pn according to the equation (3) (step S14). Specifically, the control unit 11 converts the sound pressure vector p ^n-1 at the (n-1)th time step and the differential vector v ^n-1 at the (n-1)th time step into ( 3) Calculate the sound pressure vector p ⁿ at the n-th time step by substituting into the equation.

After calculating the sound pressure vector p ⁿ , the control unit 11 secondly functions as an auxiliary function calculation unit 124 to calculate the auxiliary functions φ _i , ψ _i ⁽¹⁾ , ψ _i ⁽²⁾ at the n-th time step as Calculate (step S15). Specifically, the control unit 11 calculates the values of the auxiliary functions φ _i , ψ _i ⁽¹⁾ , ψ _i ⁽²⁾ at the (n−1)th time step and Assistance at the n ^- th time step by solving the ordinary differential equations of Eqs. Compute the values of the functions φ _i , φ _i ⁽¹⁾ , φ _i ⁽²⁾ . The differential vector v ^n-1 at the (n-1)th time step is used as the time first-order differential of the sound pressure vector p ^n-1 at the (n-1)th time step.

After calculating the auxiliary functions φ _i , φ _i ⁽¹⁾ and φ _i ⁽²⁾ , the control unit 11 thirdly functions as the differential vector calculator 125 and calculates the differential vector v ⁿ according to equation (14). (Step S16). Specifically, the control unit 11 controls the differential vector v ⁿ⁻¹ at the ( ⁿ −1)th time step, the external force vector fn at the nth time step, and the nth time calculated at step S14. By substituting the sound pressure vector p ⁿ at the step and the auxiliary functions φ _i , ψ _i ⁽¹⁾ , ψ _i ⁽²⁾ at the n-th time step calculated in step S15 into equation (14), Calculate the differential vector v ⁿ at the nth time step.

After calculating the differential vector vn, the control unit 11 determines whether or not the value of ⁿ is smaller than the total number of steps Ns (step S17). If the value of n is smaller than the total number of steps Ns (step S17; YES), the controller 11 increments the value of n (step S18). Then, the control unit 11 functions as the repeating unit 127, returns the process to step S14, and executes the processes of steps S14 to S17 again.

Specifically, the control unit 11 converts the sound pressure vector ^pn , the auxiliary functions φ _i , ψ _i ⁽¹⁾ , ψ _i ⁽²⁾ and the differential vector v ⁿ already calculated with the value of n before incrementing to are used to calculate the sound pressure vector p ⁿ , the auxiliary functions φ _i , φ _i ⁽¹⁾ , φ _i ⁽²⁾ and the derivative vector v ⁿ at the incremented value of n. Then, the control unit 11 determines whether or not the value of n has reached the total number of steps Ns. If the value of n has not reached the total number of steps Ns, the value of n is incremented and again step S14 to step S14. The process of S17 is executed. In this way, the control unit 11 calculates the sound pressure vector p ⁿ and the differential vector v ⁿ while increasing the value of n by 1 until the value of n reaches the total number of steps Ns. Calculate the time progress of a sound wave in space.

Finally, when the value of n reaches the total number of steps Ns (step S17; NO), the control section 11 functions as the output section 130 and outputs the calculation result (step S19). For example, the control unit 11 repeats the processing of steps S14 to S17 to generate an image showing the time waveforms and sound pressure distributions of sound waves represented by sound pressure vectors p ¹ to p ^Ns calculated. Then, as shown in FIGS. 4 and 5, the control unit 11 displays an image showing the temporal change and spatial distribution of the sound pressure on the display unit 14. FIG.

As described above, the sound field analysis device 10 according to the present embodiment is a device for analyzing a sound field using the time domain finite element method based on the first-order ordinary differential equation. A time update formula for the sound pressure vector p obtained by discretizing the first derivative by difference approximation, the inverse Fourier transform value of the frequency-dependent specific acoustic admittance ratio y(ω), and the time first derivative of the sound pressure vector p The value of the sound pressure vector p at a plurality of time steps is calculated according to the time update formula of the differential vector v obtained by approximating the convolution with the ADE method. Since the ADE method approximates the convolution, which requires a large computational load, the sound field analysis apparatus 10 according to the present embodiment can efficiently handle frequency-dependent impedance boundaries.

In particular, in the conventional method using the frequency-independent surface impedance model, it was necessary to divide the calculation into multiple bands in order to reflect the frequency dependence of sound absorption. On the other hand, the sound field analysis apparatus 10 according to the present embodiment can reflect the frequency dependence of sound absorption with higher accuracy in one calculation.

In addition, the sound field analysis apparatus 10 according to the present embodiment solves the ordinary differential equations (9), (10), and (11) associated with the ADE method by the fourth-order accuracy explicit Runge-Kutta method, It is possible to accurately and stably reflect the frequency dependence of sound absorption.

(Numerical experiment)
Next, numerical experiments were conducted to evaluate how accurately the sound pressure can be estimated by the sound field analysis apparatus 10 according to the above embodiment. Thus, the accuracy improvement effect of the sound field analysis device 10 according to the above embodiment was verified.

FIG. 7 shows a cross section of the acoustic tube model used as the space to be analyzed in the numerical experiments. In the numerical experiment, as shown in FIG. 7, the sound field inside an acoustic tube with a width of 0.01 m×0.01 m and a length of 0.5 m having a frequency-dependent impedance boundary at the end of the tube was analyzed.

The frequency-dependent impedance boundary was given the surface impedance of the back rigid-wall porous sound absorber modeled by Miki's equation. For details of Miki's formula, see "Y. Miki, "Acoustical properties of porous materials - Modification of Delany-Bazley models -," J. Acoust. Soc. Jpn. (E), 11, 19-24 (1990). ."It is described in. The thickness and unit thickness flow resistance of the porous material were set to 0.025 m, respectively, and R=13,000 Pa s/m ² . The specific acoustic admittance ratio y(ω) of this porous sound absorber was approximated by a rational function with n _rp =n _cp =4.

FIG. 8 shows a comparison between the sound absorption coefficient of the porous sound absorber calculated by the rational function model and the theoretical value. In FIG. 8, the vertical axis represents the sound absorption coefficient, and the horizontal axis represents the frequency. As shown in FIG. 8, the calculated values of the rational function model indicated by the dashed line and the theoretical values indicated by the solid line are in good agreement.

The acoustic tube model was discretized with 6-node cubic elements with an element length of 0.005 m. As an initial condition of the analysis, we assumed a plane wave incident on the entrance of the acoustic tube. The analysis upper limit frequency was set to 10 kHz. Also, the time step width Δt was set to Δt=l·Δt _LIMIT . Here, Δt _LIMIT is the stability limit value of the time step size Δt of the explicit time domain finite element method in the calculation of the air domain when using cubic elements (Non-Patent Document 1). Also, l is a coefficient obtained numerically. For l=1, the stability was comparable to the air region, and for l<1, the stability was worse than the air region.

FIG. 9 shows a comparison between the sound absorption coefficient of the porous sound absorber calculated by the method of the embodiment and the theoretical value. In FIG. 9, as in FIG. 8, the vertical axis represents the sound absorption coefficient, and the horizontal axis represents the frequency. As shown in FIG. 9, there is good agreement between the calculated values of the method of the embodiment, indicated by the dashed line, and its theoretical values, indicated by the solid line.

In addition, in FIG. 9, in order to show the effectiveness of the method of the embodiment, an example of the calculation result by the conventional method is also shown by a dashed line. The conventional method shows a case where boundary sound absorption is analyzed by dividing it into six octave bands from 250 Hz to 8 kHz using a frequency-independent surface impedance model. It can be seen that the conventional method requires six calculations and does not properly reflect the frequency characteristics of sound absorption. On the other hand, by using the method of the embodiment, it can be confirmed that the broadband sound absorption characteristics can be reflected with high accuracy in one calculation. Moreover, the method of the embodiment can be stably calculated with l=0.85, and is excellent in terms of stability.

(Modification)
Although the embodiment of the present invention has been described above, the above embodiment is an example, and the scope of application of the present invention is not limited to this. That is, the embodiments of the present invention can be applied in various ways, and all embodiments are included in the scope of the present invention.

For example, in the above-described embodiment, the time update formula for the sound pressure vector p shown in formula (3) was obtained by discretizing the time first-order differential of the sound pressure vector p by forward difference approximation with first-order precision. In addition, the time update formula for the differential vector v shown in equation (14) discretizes the time first derivative of the differential vector v by backward difference approximation with first-order precision, and furthermore, the time first derivative of the sound pressure vector p is quadratically The accuracy was obtained by discretizing with the median difference approximation. However, the time-first derivative of the sound pressure vector p and the differential vector v may be discretized by other methods.

In the above embodiment, the matrix setting unit 121 calculates the Courant number τ from the time step size Δt and the element length h, following the modified integration rule disclosed in Non-Patent Document 2, and calculates the numerical integration points α _m , α _k It was set. However, the matrix setting unit 121 may set the numerical integration points α _m and α _k using another method. For example, the matrix setting unit 121 is disclosed in “Takumi Yoshida, Takeshi Okuzono, Kimihiro Sakagami, “Time domain room acoustic solver with fourth-order explicit FEM using modified time integration,” Applied Sciences, 10, 3750 (2020).” Alternatively, numerical integration points α _m and α _k that do not depend on the shape of the elements that discretize the space to be analyzed may be used. In this case, the numerical integration points α _m and α _k are set as “α _m =±√(4/3)” and “α _k =±√(2/3)” without using the Courant number τ. .

In the above embodiment, the sound field analysis apparatus 10 includes the display unit 14, and the output unit 130 displays output information based on the sound pressure vector p at a plurality of time steps calculated by the numerical calculation unit 120. 14 was displayed. However, the output section 130 may output the output information to a device external to the sound field analysis device 10 via communication by the communication section 15 . Then, an external device may display the output information output from the output unit 130, as shown in FIG. 4 or 5, for example. When the output information is output to an external device in this way, the sound field analysis device 10 does not have to include the display section 14 .

In the above-described embodiment, in the control unit 11 of the sound field analysis apparatus 10, the CPU executes programs stored in the ROM or the storage unit 12, so that each of the setting reception unit 110, the numerical calculation unit 120, and the output unit 130 worked. However, the controller 11 may be dedicated hardware. Dedicated hardware is, for example, a single circuit, a composite circuit, a programmed processor, an ASIC (Application Specific Integrated Circuit), an FPGA (Field-Programmable Gate Array), a combination thereof, or the like. When the control unit 11 is dedicated hardware, each function of each unit may be realized by separate hardware, or the functions of each unit may be collectively realized by single hardware. Also, some of the functions of each unit may be realized by dedicated hardware, and other parts may be realized by software or firmware. In this way, the control unit 11 can implement each of the functions described above using hardware, software, firmware, or a combination thereof.

By applying an operation program that defines the operation of the sound field analysis device according to the present invention to an existing personal computer, information terminal device, or the like, the personal computer, information terminal device, or the like can be used as the sound field analysis device according to the present invention. It is also possible to make it work.

Any method of distributing such programs may be used. For example, computer-readable recording media such as CD-ROMs (Compact Disk ROMs), DVDs (Digital Versatile Disks), MOs (Magneto Optical Disks), and memory cards may be used. It may be stored and distributed, or may be distributed via a communication network such as the Internet.

The present invention is capable of various embodiments and modifications without departing from the broader spirit and scope of the invention. Moreover, the above-described embodiments are for explaining the present invention, and do not limit the scope of the present invention. That is, the scope of the present invention is indicated by the claims rather than the embodiments. Various modifications made within the scope of the claims and within the meaning of the invention equivalent thereto are considered to be within the scope of the present invention.

10 Sound field analysis device 11 Control unit 12 Storage unit 13 Operation unit 14 Display unit 15 Communication unit 20 Space 110 Setting reception unit 120 Numeric calculation unit 121 Matrix setting unit 123 Sound pressure vector calculation unit 124 Auxiliary function calculation unit 125 Differential vector calculation unit 127 repeating unit 130 output unit

Claims (7)

A first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in a space to be analyzed and a differential vector equivalent to the time first order differential of said sound pressure vector, wherein the sound wave at the boundary surrounding said space A sound field analysis device for analyzing a sound field using a time domain finite element method based on the first-order ordinary differential equation including an attenuation term for evaluating the attenuation of
In the first-order ordinary differential equation, the time update formula of the sound pressure vector obtained by discretizing the time first-order differential of the sound pressure vector by difference approximation, and in the first-order ordinary differential equation, the attenuation term A time update formula for the differential vector obtained by approximating the convolution of the included frequency-dependent inverse Fourier transform value of the specific acoustic admittance ratio with the time first derivative of the sound pressure vector using the Auxiliary Differential Equations method, according to , a numerical calculator that calculates the values of the sound pressure vector at a plurality of time steps;
Sound field analyzer. The differential vector time update formula includes an auxiliary function obtained by rational function approximation of the specific acoustic admittance ratio,
The numerical calculation unit calculates the time evolution of the auxiliary function by a time integration method,
The sound field analysis device according to claim 1. The numerical calculation unit calculates the time evolution of the auxiliary function using the Runge-Kutta method, the Euler method, the Crank-Nicholson method, or the Huyn method as the time integration method.
The sound field analysis device according to claim 2. The time update formula for the sound pressure vector is obtained by combining the value of the sound pressure vector at the nth time step with the value of the sound pressure vector at the (n−1)th time step and the (n−1)th time step. and the value of the differential vector at
The time update formula for the differential vector is such that the value of the differential vector at the nth time step is combined with the value of the sound pressure vector at the (n−1)th time step and the nth time step, and the (n−1)th ) is defined by the value of the differential vector at the time step and the value of the auxiliary function at the nth time step,
The sound field analysis device according to claim 2 or 3. The time update formula for the differential vector discretizes the time first derivative of the differential vector by backward difference approximation with first-order accuracy, and further discretizes the time first derivative of the sound pressure vector by central difference approximation with second-order accuracy. obtained by
The sound field analysis device according to any one of claims 1 to 4. A first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in a space to be analyzed and a differential vector equivalent to the time first order differential of said sound pressure vector, wherein the sound wave at the boundary surrounding said space A sound field analysis method for analyzing a sound field using a time domain finite element method based on the first-order ordinary differential equation including an attenuation term for evaluating the attenuation of
In the first-order ordinary differential equation, the time update formula of the sound pressure vector obtained by discretizing the time first-order differential of the sound pressure vector by difference approximation, and in the first-order ordinary differential equation, the attenuation term A time update formula for the differential vector obtained by approximating the convolution of the included frequency-dependent inverse Fourier transform value of the specific acoustic admittance ratio with the time first derivative of the sound pressure vector using the Auxiliary Differential Equations method, according to , calculating values of said sound pressure vector at multiple time steps;
Sound field analysis method. the computer,
A first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in a space to be analyzed and a differential vector equivalent to the time first order differential of said sound pressure vector, wherein the sound wave at the boundary surrounding said space A program that functions as a sound field analysis device that analyzes a sound field using a time domain finite element method based on the first-order ordinary differential equation including an attenuation term for evaluating the attenuation of
said computer,
In the first-order ordinary differential equation, the time update formula of the sound pressure vector obtained by discretizing the time first-order differential of the sound pressure vector by difference approximation, and in the first-order ordinary differential equation, the attenuation term A time update formula for the differential vector obtained by approximating the convolution of the included frequency-dependent inverse Fourier transform value of the specific acoustic admittance ratio with the time first derivative of the sound pressure vector using the Auxiliary Differential Equations method, according to , a numerical calculator that calculates the values of the sound pressure vector at multiple time steps;
program.

Images