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

Abstract

To improve analysis accuracy 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. 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 applying a linear multistage method optimized so that target accuracy is achieved with respect to a time and space discretization error and an amplitude error is not included in the discretization error in the first order normal differential equation, and a time update equation of the differential vector obtained by discretizing a time first derivative of the differential vector by a difference approximation in the first order normal differential equation.SELECTED DRAWING: Figure 1

Description

Application for application of Article 30, Paragraph 2 of the Patent Law Issued by: Architectural Institute of Japan, Publication name: Summaries of technical papers of academic lectures at the 2020 Annual Meeting of the Architectural Institute of Japan (Kanto), Date of issue: July 20, 2020 Day

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 to 3 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)
Simultaneous first-order ordinary derivative obtained by introducing diagonal mass matrix and time-first derivative vector of sound pressure into semi-discretized equation for sound field surrounded by rigid wall boundary or vibrating boundary Targets 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", and "D" in formulas (1) and (2) represent the mass matrix, stiffness matrix, and diagonal mass matrix, respectively. The mass matrix M, stiffness matrix K, and diagonal mass matrix D are global matrices formed by overlapping the element mass matrix M _e , the element stiffness matrix K _e , and the element diagonal mass matrix D _e , respectively. 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.

(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 precision. 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 size. 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).

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.

Non-Patent Document 3 discloses a modified Adams method, which is a spatio-temporal 4th order accuracy method for handling arbitrarily shaped elements. In order to achieve high precision by the modified integral rule disclosed in Non-Patent Document 2, the shapes of usable elements are limited to squares or cubes, whereas the modified Adams The method can reduce discretization errors in time and space to fourth order accuracy when using elements with shapes other than squares or cubes. As a result, even when an indoor space or an acoustic material including a complicated shape is to be analyzed, it is possible to accurately reflect the characteristics due to the shape.

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) 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).

Improving computational efficiency is important in simulating acoustic wave propagation using computational mechanics methods including the time-domain finite element method. Conventional explicit time-domain finite element method acoustic analysis methods have excellent computational efficiency for air domain calculations by applying modified integral rules. In particular, the modified Adams method disclosed in Non-Patent Document 3 can handle arbitrarily shaped elements.

However, the conventional approach involves amplitude error, which is a type of discretization error. The occurrence of amplitude errors leads to deterioration of the accuracy of sound field analysis, especially when predicting high-frequency sound wave propagation over a long period of time. Therefore, a method capable of suppressing the occurrence of amplitude errors and improving analysis accuracy is desired.

The present invention is intended to solve the above problems, and provides a sound field analysis device, a sound field analysis method, and a program capable of improving analysis accuracy in sound field analysis using the time domain finite element method. intended to provide

In order to achieve the above object, the sound field analysis device according to the first aspect of the present invention includes:
A sound field is generated using a time domain finite element method based on a first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in the space to be analyzed and a differential vector equivalent to the time first order differential of the sound pressure vector. A sound field analysis device for analyzing
Obtained by applying a linear multi-stage method optimized so that the discretization error in time and space does not include an amplitude error and the target accuracy is achieved in the first-order ordinary differential equation and a time update formula for the differential vector obtained by discretizing the time first order differential of the differential vector in the first order ordinary differential equation using difference approximation, according to a plurality of a numerical calculator for calculating the value of the sound pressure vector at time steps of .

The time update formula for the sound pressure vector is such that the value of the sound pressure vector at the n-th time step is the value of the sound pressure vector at the (n−M) time step to the (n−1) time step. , and the value of the differential vector from the (n−N)th time step to the (n−1)th time step are multiplied by a weighting factor optimized to the value in which the imaginary part is not included in the discretization error It may be a formula calculated by adding above.

When the target accuracy is K-th order accuracy, the M and N may be set to a value equal to or greater than K−1.

Both the sum of the weighting factors for the sound pressure vector and the sum of the weighting factors for the differential vector may be equal to one.

The differential vector time update formula is such that the value of the differential vector at the n-th time step is combined with the value of the differential vector at the (n-1)th time step and the value of the sound pressure vector at the n-th time step. It may be an expression determined by a value.

The first-order ordinary differential equation may include a mass matrix and a stiffness matrix calculated using numerical integration points independent of the shape of the space-discretizing elements.

In order to achieve the above object, a sound field analysis method according to a second aspect of the present invention comprises:
A sound field is generated using a time domain finite element method based on a first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in the space to be analyzed and a differential vector equivalent to the time first order differential of the sound pressure vector. A sound field analysis method for analyzing
Obtained by applying a linear multi-stage method optimized so that the discretization error in time and space does not include an amplitude error and the target accuracy is achieved in the first-order ordinary differential equation and a time update formula for the differential vector obtained by discretizing the time first order differential of the differential vector in the first order ordinary differential equation using difference approximation, according to a plurality of Calculate the value of the sound pressure vector at time steps of .

In order to achieve the above object, a program according to the third aspect of the present invention,
the computer,
A sound field is generated using a time domain finite element method based on a first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in the space to be analyzed and a differential vector equivalent to the time first order differential of the sound pressure vector. A program that functions as a sound field analysis device that analyzes
said computer,
Obtained by applying a linear multi-stage method optimized so that the discretization error in time and space does not include an amplitude error and the target accuracy is achieved in the first-order ordinary differential equation and a time update formula for the differential vector obtained by discretizing the time first order differential of the differential vector in the first order ordinary differential equation using difference approximation, according to a plurality of A numerical calculation unit that calculates the value of the sound pressure vector at time steps of .

According to the present invention, the accuracy of analysis is improved 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 model used as a space to be analyzed in numerical experiments; FIG. 5 is a diagram showing a comparison of amplitude errors included in 4.5 kHz sound waves calculated using the method of the embodiment and the modified Adams 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 sound pressure distribution 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 the setting of either a rigid wall boundary or a vibration 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.

Third, the setting reception unit 110 receives mesh settings. Here, a mesh is obtained by dividing (discretizing) a 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.

As will be described later, in this embodiment, the numerical calculation unit 120 calculates the value of the sound pressure vector p using numerical integration points that do not depend on the shapes of the elements. Therefore, the shape of the element is not limited to a square, and may be an arbitrary quadrilateral or an arbitrary hexahedron.

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 for the sound pressure vector p and the time update formula for the differential vector v obtained from the differential equation. A first-order ordinary differential equation is an equation including a mass matrix M and a stiffness matrix K, as shown in equations (1) and (2).

Here, a calculation method by the numerical calculation unit 120 will be described. In the present embodiment, the numerical calculation unit 120 uses Equation (4) described in Background Art as the time update equation for the differential vector v. Equation (4) is a time updating equation for the differential vector v obtained by discretizing the time first-order differential of the differential vector v in the equation (2) with the differential approximation of the first-order precision, as described in Background Art. be.

On the other hand, the numerical calculation unit 120 uses, as the time update formula for the sound pressure vector p, a time update formula different from the formula (3) described in the background art. Specifically, the numerical calculation unit 120 achieves the target accuracy for the discretization error in time and space, and applies a linear multistage method optimized so that the discretization error does not include the amplitude error. A time update formula for the sound pressure vector p obtained by is used.

Derivation of the time update formula for the sound pressure vector p in this embodiment will be described below. The numerical calculation unit 120 applies the modified integral rule disclosed in Non-Patent Document 2 to the calculation of the elements of the mass matrix M and the stiffness matrix K. Furthermore, the numerical calculator 120 applies a linear multi-step method to calculate the time evolution of the sound pressure vector p. Specifically, the numerical calculation unit 120 uses an explicit linear multi-stage time integration formula represented by the following formula (5) as a time update formula for the sound pressure vector p.

In the equation (5), M and N are natural numbers representing the number of time steps of the past sound pressure vector p and differential vector v used to calculate the time evolution. As shown in equation (5), the time update formula for the sound pressure vector p is to change the value of the sound pressure vector p ⁿ at the nth time step from the (n−M)th time step to the (n−1)th time step. The values of the sound pressure ^vectors p n ^{- 1} , p ^n-2 , . , v ⁿ⁻² , . . . , v ^n−N .

Also, in equation (5), the coefficients a _i and b _j are weighting coefficients for the values of the sound pressure vector p and the differential vector v at each time step. In order to ensure stability of calculation, the coefficients a _i and b _j must satisfy the following equation (6). That is, _both the _{sum a 1 +} a ₂ + .

The values of M and N in equations (4) and (5), i.e. the number of past time steps for calculating the sound pressure vector p ⁿ and differential vector v ⁿ at the nth time step, are the discretization error to be achieved is set according to the accuracy of Specifically, in order to achieve the K-th order accuracy (K is a natural number) as the target accuracy of the discretization error in time and space, the values of M and N must be K−1 or more.

As an example, if we want to achieve fourth order accuracy for the discretization error, we need M≧3 and N≧3. A case where M=N=3 will be described below as an example. When M=N=3, the following equation (7) is derived by eliminating the differential vector v of each time step from equation (5) using equation (4). At this time, the term of the external force vector ^fn in the equation (4) is removed.

As a condition for the discretization error to contain no amplitude error, the weighting coefficients a _i and b _j must be optimized to values that do not contain imaginary parts in the discretization error. Specifically, when the coefficients of the sound pressure vector are symmetrical about the sound pressure vector at time n, the imaginary part is not included in the discretization error. That is, when M=N=3, the condition for the discretization error not to include the amplitude error is that the coefficient 1 of p ⁿ⁺² and the coefficient a ₃ of p ⁿ⁻² in the first half of equation (7) are equal, The coefficient (a ₁ +1) of p ⁿ⁺¹ and the coefficient (a ₃ −a ₂ ) of p ⁿ⁻¹ are equal, and the coefficient b ₁ of p ⁿ⁺¹ and the coefficient b of p ⁿ⁻¹ in the second half of equation (7) ₃ must be equal. Therefore, it is necessary to satisfy the following expression (8).

The following equation (9) is determined from equations (6) and (8).

From the expression (9), the expression (7) is expressed as the following expression (10).

Here, assuming plane wave propagation in free space discretized by cubic elements of element size h, the dispersion characteristics of equation (10) are evaluated. A plane wave in the time domain in a three-dimensional sound field is represented by the following equation (11). For details of evaluation of dispersion characteristics, see, for example, “Takeshi Okuzono, Takumi Yoshida, Kimihiro Sakagami, Toru Otsuru, An explicit time-domain finite element method for room acoustics simulations: Comparison of the performance with implicit methods, Applied Acoustics, 104, 76-84 (2016)”.

(11), k represents the wavenumber and ω represents the angular frequency. θ and φ represent elevation and azimuth angles in a spherical coordinate system. Also, i is an imaginary unit. From the expression (11), the expression (10) is expressed as the following expression (12).

From the equation (12), a dispersion relational expression expressed by the following equation (13) is obtained. Note that M _c and K _c in the formula (13) are expressed as the following formulas (14), (15), and (16).

In equation (13), ^{ch and k h represent} the numerical sound velocity and wave number, respectively. In equations (14) and (15), α _m and α _k are integration points used for numerical integration of the element matrices of the mass matrix M and stiffness matrix K, respectively. To achieve spatial fourth-order accuracy, the numerical integration points shown in equation (17) below are used. The numerical integration point shown in Equation (17) is a numerical integration point that does not depend on the shape of the elements that discretize the space to be analyzed, as disclosed in Non-Patent Document 3.

By Taylor-expanding the equation (13) with respect to k ^h using the integration points of the equation (17), the following equation (18) is obtained as an equation representing the dispersion error (discretization error).

Therefore, in order to achieve the fourth-order temporal accuracy as the target accuracy, it is sufficient that the third-order or lower terms related to the wave number k become zero in the equation (18). Therefore, it suffices if b1 becomes _a value indicated by the following equation (19).

The variance error (discretization error) up to the fifth-order term of the explicit time domain finite element method using the equation (19) is expressed as the following equation (20). In addition, in the equation (20), the value of A is expressed as in the following equation (21).

Equation (20) shows that the explicit time-domain finite element method using the parameters of equation (19) achieves fourth order accuracy in time and space. Also, since the equation (20) does not have an imaginary part, no amplitude error occurs. Note that the equation (13) becomes a real function that does not include the imaginary part due to the condition of the equation (8), so the method of the present embodiment does not show the imaginary part even in the dispersion error of the sixth order or higher. Note that the dispersion error (discretization error) in the three-dimensional analysis of the modified Adams method including the amplitude error is calculated by the following equation (22).

In the variance error in the equation (22), an imaginary part, which is an imaginary term, appears in the quintic term. Therefore, it can be seen that an amplitude error occurs. Also, from equations (20) and (22), both the method of the present embodiment and the modified Adams method achieve fourth-order accuracy in time and space. It can be seen that the method of the present embodiment has performance equivalent to the modified Adams method with respect to numerical variance.

The numerical calculation unit 120 substitutes the values of the formulas (9) and (19) into the weighting coefficients a _i and b _j of the formula (5), which is the time update formula for the sound pressure vector p derived in this way. The values of the sound pressure vector p at a plurality of time steps are calculated according to the equation and the equation (4), which is the time update equation for the differential vector v.

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 functions of a matrix setting unit 121 , a sound pressure vector calculation unit 123 , 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, and the diagonal mass matrix D included in the formulas (5) and (4). The mass matrix M, stiffness matrix K, 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 the values of the elements to be analyzed. It is determined depending on the dimension and shape of the space and the shape of the discretized element. Therefore, the matrix setting unit 121 sets the mass matrix M, the stiffness matrix K, 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 integration rule disclosed in Non-Patent Document 2. That is, the matrix setting unit 121 calculates the mass matrix M and the stiffness matrix K using the numerical integration points α _m and α _k that do not depend on the shape of the elements that discretize the space to be analyzed shown in Equation (17). do. Specifically, the matrix setting unit 121 uses the numerical integration points α _m and α _k to calculate the element mass matrix M _e and the element stiffness matrix K _e by the Gauss-Legendre quadrature method, Set the stiffness matrix K.

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.

The sound pressure vector calculation unit ₁₂₃ calculates the _n -th Compute the sound pressure vector p ⁿ at the time step. 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, equation (5) expresses the value of the sound pressure vector p ⁿ at the n-th time step as the sound pressure vector p ^{n −M} ~ p ^n-1 and the values of differential vectors v ^n-N ~ v ^n-1 from the (nN)th time step to the (n-1)th time step are imaginary to the discretization error. It is an expression for calculating by multiplying optimized weighting coefficients a _i and b _j to values that do not include parts and then adding them. Therefore, the sound pressure vector calculator 123 calculates the value of the sound pressure vector p ⁿ at the n-th time step as the sound pressure vector p ^{n- 1} , ^{p n−2 , .} calculated from the value of ^nN . Specifically, when N=M=3, the sound pressure vector calculator 123 calculates the values of the sound pressure vectors p ⁿ at the n-th time step as sound pressure vectors p ⁿ⁻¹ , p ⁿ⁻² , p It is calculated from the value of ^n-3 and the values of the differential vectors v ^n-1 , v ^n-2 and v ^n-3 .

The differential vector calculator 125 calculates the differential vector v ⁿ at the n-th time step according to Equation (4), 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 (4) is obtained by ^dividing the value of the differential vector vn at the nth time step into the value of the differential vector vn−1 at the (n−1)th time step and the value of the differential vector vn ⁻¹ at the nth time step and the value of the sound pressure vector ^pn . Therefore, the differential vector calculator 125 calculates the value of the differential vector vn at the ^{n -} th time step, the value of the differential vector vn-1 at the (n-1)th time step, and the sound at the n-th time step. and the value of the pressure vector ^pn .

The repetition unit 127 performs a process of 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, and the differential vector calculation unit 125 and the process of calculating the value of the differential vector v ⁿ at the step are repeated. Specifically, the repeating unit 127 changes the value of the sound pressure vector p ^k and the differential vector v ^k when n=k while changing the value of n from 1 to the total number of steps Ns. Values of sound pressure vectors p ^k-1 to p ^k-M calculated when k-1 to kM, and differential vector v ^k-1 calculated when n=k-1 to kN ˜v ^k−N values. 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.

Sound pressure vectors p ^n-1 to p ^n-M and differential vectors v ^n-1 to v ^n-N when n=1, and sound pressure vectors p ^n-2 to p ^n- when n=2 ^M and derivative vectors v ⁿ⁻² to v ^n−N , . . . etc. 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 derivative vector ^vn according to the equations (5) and (4).

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 sets the dimension and shape of the space to be analyzed, the boundary conditions, the shape and size of the elements, the time step width Δt, the total number of steps Ns, the external force vector f, and the like. from the user via

Upon receiving the settings, the control unit 11 functions as the matrix setting unit 121 and sets the mass matrix M, stiffness matrix K, and diagonal mass matrix D (step S12). As in Non-Patent Document 3, the control unit 11 sets “α _m =±√(4/3 ₎ ” and “α _{k =} ± √(2/3)" is used. 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.

After setting the mass matrix M, the stiffness matrix K, and the diagonal mass matrix D, 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 controller 11 functions as the sound pressure vector calculator 123 and calculates the sound pressure vector ^pn (step S14). Specifically, the control unit 11 controls sound pressure vectors p ^n-1 , p ^{n -2} , . and the values of differential vectors v _n ⁻¹ , v ^{n −2} , . By substituting _j into the optimized equation (5), the sound pressure vector p ⁿ at the n-th time step is calculated.

After calculating the sound pressure vector ^pn , the control unit 11 secondly functions as the differential vector calculator 125 to calculate the differential vector ^vn (step S15). 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 into equation (4), the differential vector v ⁿ at the n-th time step is calculated.

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 S16). If the value of n is smaller than the total number of steps Ns (step S16; YES), the controller 11 increments the value of n (step S17). 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 S16 again.

Specifically, the control unit 11 uses the sound pressure vector p n and the differential vector v n already calculated for the value of n before incrementing, and uses the sound pressure vector p ⁿ and the differential vector v ⁿ for the value of ⁿ after incrementing. Compute the vector ^vn . 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 S16 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 S16; NO), the control section 11 functions as the output section 130 and outputs the calculation result (step S18). For example, the control unit 11 repeats the processing of steps S14 to S16 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 apparatus 10 according to the present embodiment achieves the target accuracy with respect to the discretization error, and uses the linear multistage method optimized so that the discretization error does not include the amplitude error. A plurality of time steps according to the time update formula for the sound pressure vector p obtained by applying the Calculate the value of the sound pressure vector p at As a result, the sound field analysis apparatus 10 according to the present embodiment can suppress the occurrence of amplitude errors and improve analysis accuracy.

In particular, the time domain finite element method produces a discretization error, which is an error due to the discretization of time and space. Therefore, in order to suppress the discretization error and make the calculation results reasonable, a sufficiently finer mesh than the wavelength of the sound wave to be analyzed is required. However, when analyzing sound waves in the audible range, a very fine computational mesh is required, increasing the computational load. Therefore, in order to utilize computational mechanics methods in the field of acoustics, it is necessary to improve computational efficiency by developing a computational algorithm that can suppress errors even with relatively coarse discretization and perform appropriate computations. is important.

The modified Adams method disclosed in Non-Patent Document 3 can handle arbitrarily shaped elements and can reduce time and space discretization errors to fourth-order accuracy, but it contains amplitude errors. . However, the modified Adams method contains an amplitude error. Therefore, when predicting the propagation of high-frequency sound waves over a long period of time, such as when predicting the reverberation time, which is a typical parameter of architectural acoustics, the prediction accuracy deteriorates. In contrast, the method of the present embodiment achieves fourth order accuracy in time and space, comparable to the modified Adams method. Moreover, the method of the present embodiment does not include the amplitude error, which is a problem in the modified Adams method. Therefore, the sound field analysis device 10 according to the present embodiment can analyze the sound field with high accuracy using a more versatile technique.

(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 the model used as the space to be analyzed in the numerical experiments. In the numerical experiment, as shown in FIG. 7, a geometrically unattenuated sound field in which a plane wave propagates in a long and narrow pipe with a length of 100 m and a width of 0.02 m was analyzed. The plane of the tube at x=0 m was vibrated by applying a Gaussian pulse with an upper limit frequency of 4.5 kHz. In addition, the plane at x=0 m of the tube was taken as an absorbing boundary by giving it the characteristic impedance of air.

The time response for 0.2 s in the range of x=1 to 60 m in such a tubular model was calculated by the method of the embodiment and the modified Adams method, respectively. Here, the element length was set to 0.01 m. Also, the time step size was set to 1/71000 s and 1/94000 s for the method of the embodiment and the modified Adams method, respectively. The sound pressure level was calculated by transforming the analysis result into a frequency response by discrete Fourier transform.

In this numerical experiment, the numerical amplitude error e _n (x) at a point xm away from the plane of x=0m, which is the sound source, is defined as in the following equation (23). (23), L _n (x) is the numerical solution for the sound pressure level at a point xm away from the sound source in the pipe. Specifically, the amplitude error e _n (x) is defined as the difference between the sound pressure level L _n (1) at a point 1 m from the sound source in the pipe and the sound pressure level L _n (x) at a point x m from the sound source.

FIG. 8 shows a comparison of the amplitude error e _n (x) included in the 4.5 kHz sound wave calculated using the method of the embodiment and the modified Adams method. As indicated by the solid line in FIG. 8, in the modified Adams method, an amplitude error of 4 dB or more occurs despite sound wave propagation for about 0.2 s. In order to reduce this amplitude error to 0.5 dB or less, it is necessary to reduce the time step width to 1/142000 s or less.

On the other hand, according to the method of the embodiment, no amplitude error occurs, as indicated by the dashed line in FIG. Therefore, it was verified by numerical simulations of plane wave propagation in an elongated tube that the 4.5 kHz sound waves calculated by the method of the embodiment do not contain amplitude errors.

From the results of such numerical experiments, it was found that, in the case of analysis with the same accuracy, the method of the embodiment can perform calculations in twice as fast calculation time using approximately twice as many time steps as the modified Adams method. This advantage increases as the analysis frequency increases and as the analysis time length increases.

(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 embodiment, the case where the values of M and N in equations (4) and (5) are 3, that is, the case where fourth-order accuracy is achieved with respect to time and space discretization errors has been described as an example. However, the values of M and N can be set to any value depending on the accuracy to be achieved.

When the values of M and N are arbitrary, the weighting coefficients a _i and b _j in equation (5), which is the time update equation for the sound pressure vector p, are adjusted according to the values of M and N to achieve the target accuracy and the discretization error is optimized to a value that does not include the amplitude error. In order not to include the amplitude error in the discretization error, the differential vectors v ⁿ⁻¹ , v ^{n −2} , . (7)), the coefficients of the sound pressure vector need to be symmetrical with the time step n as the origin. Specifically, the weight _coefficients a _i ( _{i = 1} , ₂ , . and the _weighting coefficients b _j (j _{= 1} , ₂ , . In addition, when the target accuracy is K-th order accuracy, in the discretization error formula (Equation (18) in the case of M = N = 3), the terms of the K-1th order or lower related to the wave number k are equal to 0 must meet the conditions. The weighting coefficients a _i and b _j are optimized to values that satisfy these conditions as well as the stability condition of equation (6).

Alternatively, the setting reception unit 110 may receive the setting of the target accuracy of the discretization error from the user, and set the values of M and N according to the received target accuracy. For example, if the received target accuracy is K-th order accuracy, the setting reception unit 110 sets the values of M and N to values of K−1 or more. In this case, the numerical calculation unit 120 sets the time update formula for the sound pressure vector p according to the values of M and N set by the setting reception unit 110 . Then, the numerical calculation unit 120 calculates the time evolution of the sound pressure vector p according to the set time update formula for the sound pressure vector p and the time update formula for the differential vector v shown in formula (4).

In addition, even if M and N are not set to large values, it is possible to reduce the variance error (discretization error) by optimizing the numerical integration points α _m and α _k and the weighting coefficients a _i and b _j . is. Specifically, when M=N=3, for example, numerical integration points α _m and α _k at which the dispersion error is minimized at a specific frequency (wavenumber) are used, and then dispersion The weighting coefficients a _i (i=1, 2, . . . , M) and the weighting coefficients b _j (j=1, 2, , N) may be optimized. This can significantly reduce the dispersion error at a particular frequency, even with low order values of M and N.

In the above embodiment, the time update formula for the differential vector v shown in Equation (4) is obtained by discretizing the time first-order differential of the differential vector v by backward difference approximation with first-order accuracy. However, the time first-order differential of the differential vector v may be discretized by other methods.

In the above embodiment, the setting receiving unit 110 receives setting of either the rigid wall boundary or the vibration boundary as the boundary condition. However, the setting reception unit 110 may receive the setting of the impedance boundary as the boundary condition. Here, 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 at least a part of the boundaries surrounding the space to be analyzed includes an impedance boundary, in order to evaluate the attenuation of the sound wave reflected by the impedance boundary, the first-order ordinary differential equation shown in equation (2) is: A damping term containing the damping matrix C is required. Therefore, the sound pressure vector calculation unit 123 calculates the differential vector ^vn according to a modified expression including an attenuation term instead of the expression (4).

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 125 Differential vector calculation unit 127 Repeat unit 130 Output Department

Claims (8)

A sound field is generated using a time domain finite element method based on a first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in the space to be analyzed and a differential vector equivalent to the time first order differential of the sound pressure vector. A sound field analysis device for analyzing
Obtained by applying a linear multi-stage method optimized so that the discretization error in time and space does not include an amplitude error and the target accuracy is achieved in the first-order ordinary differential equation and a time update formula for the differential vector obtained by discretizing the time first order differential of the differential vector in the first order ordinary differential equation using difference approximation, according to a plurality of a numerical calculator that calculates the value of the sound pressure vector at time steps of
Sound field analyzer. The time update formula for the sound pressure vector is such that the value of the sound pressure vector at the n-th time step is the value of the sound pressure vector at the (n−M) time step to the (n−1) time step. , and the value of the differential vector from the (n−N)th time step to the (n−1)th time step are multiplied by a weighting factor optimized to the value in which the imaginary part is not included in the discretization error is an expression that computes by adding over
The sound field analysis device according to claim 1. If the target accuracy is K-th order accuracy, the M and N are set to a value equal to or greater than K−1;
The sound field analysis device according to claim 2. the sum of the weighting factors for the sound pressure vector and the sum of the weighting factors for the derivative vector are both equal to 1;
The sound field analysis device according to claim 2 or 3. The differential vector time update formula is such that the value of the differential vector at the n-th time step is combined with the value of the differential vector at the (n-1)th time step and the value of the sound pressure vector at the n-th time step. is an expression determined by a value and
The sound field analysis device according to any one of claims 1 to 4. The first-order ordinary differential equation includes a mass matrix and a stiffness matrix calculated using numerical integration points that do not depend on the shape of the elements that discretize the space.
The sound field analysis device according to any one of claims 1 to 5. A sound field is generated using a time domain finite element method based on a first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in the space to be analyzed and a differential vector equivalent to the time first order differential of the sound pressure vector. A sound field analysis method for analyzing
Obtained by applying a linear multi-stage method optimized so that the discretization error in time and space does not include an amplitude error and the target accuracy is achieved in the first-order ordinary differential equation and a time update formula for the differential vector obtained by discretizing the time first order differential of the differential vector in the first order ordinary differential equation using difference approximation, according to a plurality of calculating the value of the sound pressure vector at time steps of
Sound field analysis method. the computer,
A sound field is generated using a time domain finite element method based on a first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in the space to be analyzed and a differential vector equivalent to the time first order differential of the sound pressure vector. A program that functions as a sound field analysis device that analyzes
said computer,
Obtained by applying a linear multi-stage method optimized so that the discretization error in time and space does not include an amplitude error and the target accuracy is achieved in the first-order ordinary differential equation and a time update formula for the differential vector obtained by discretizing the time first order differential of the differential vector in the first order ordinary differential equation using difference approximation, according to a plurality of functioning as a numerical calculator that calculates the value of the sound pressure vector at time steps of
program.

Images