For Dirichlet or Neumann boundary conditions, that is, on a sound-soft boundary with p = p₀ (input as PRES on the D command) or on a sound-hard boundary with , the surface integration will be eliminated.

The Robin boundary condition on impedance boundary Γ_Z is given by:

(8–38)

where:

Substituting Equation 8–38 into Equation 8–10 yields:

(8–39)

For simplicity, the inviscid fluid will be investigated. The matrix forms in Equation 8–37 are rewritten as:

(8–40)

(8–41)

In acoustic design, the attenuation coefficient (that is, the absorption coefficient) is used to define the absorption properties of the material. The attenuation coefficient is the ratio of the absorbed sound power density to the incident sound power density. The ratio is expressed as follows:

(8–42)

where:

Because the module of reflection coefficient is given as , the equivalent impedance (assumed Im(z)=0) of the attenuating material is expressed as:

(8–43)

where:

z₀ = ρ₀ c ₀ = sound impedance

The impedance value of the Robin boundary condition can be defined by the sound impedance z₀ = ρ₀ c ₀ (input as INF on the SF command), the attenuation coefficient (input as ATTN on the SF command), or the general complex impedance (admittance for modal analysis) (input as IMPD on the SF command). The Robin boundary condition can be used for the harmonic response or modal analysis.

The Floquet principle is satisfied when the acoustic wave propagates in a periodic structure (input as Lab = FPBC on the BF command or Label = PLAN on the APORT command for plane wave incidence):

(8–44)

where:

= period of infinite periodic structure

= complex propagating constant (); = attenuation constant; = phase constant

A modified eigenvalue matrix equation solves the phase shift as eigenvalues for acoustic elements if the phase shift is not equal to 0 or :

(8–45)

where:

( and are defined in Equation 8–37)

= the internal nodes

= the independent nodes for coupled independent/dependent nodes on the boundary

= dependent nodes for coupled independent/dependent nodes on the boundary

as eigenvalue

The perfectly matched layers (PML) are located inside a cubic enclosure. The PML materials are specified in face, edge, and corner regions in a three-dimensional model.

In the frequency domain, the governing lossless and source-free momentum and conservation are given by:

(8–58)

(8–59)

The anisotropic PML material is assumed as follows:

(8–60)

(8–61)

where:

In order to satisfy the vector operation, the governing acoustic equations for PML are written as:

(8–62)

(8–63)

If the modified operator is expressed as:

(8–64)

Then the wave equation for a PML with a modified differential operator is expressed as:

(8–65)

The interface between the PML and the real acoustic medium can be rendered reflectionless for all frequencies and angles of incidence of a wave propagating from the acoustic medium to the PML if the elements of [Λ] are defined as follows:

(8–66)

where:

Multiplying s_xs_ys_z with Equation 8–65 results in the final wave equation needed to construct the PML:

(8–67)

where:

By using the Galerkin method and setting the pressure to zero (input as PRES on the D command) on the backed boundary, the "weak" formulation of an acoustic pressure wave in PML is expressed as:

(8–68)

Equation 8–68 can be written in matrix notation to produce the following discretized wave equation:

(8–69)

where:

SRe = real part of complex constant sxsysz

SIm = imaginary part of complex constant sxsysz

[ΨRe] = real part of complex tensor [Ψ]

[ΨIm] = imaginary part of complex tensor [Ψ]

The PML material is defined using FLUID30, FLUID220, and FLUID221 elements with KEYOPT(4) = 1. Because the PML material is constructed in Cartesian coordinates, the edges of the 3D PML region must be aligned with the axes of the global Cartesian coordinate system.

More than one 1D PML region may exist in a finite element model. The PML element coordinate system (PSYS) uniquely identifies each PML region. A parabolic attenuation distribution is used to minimize numerical reflections in the PML. Ansys, Inc. recommends using at least three PML element layers to obtain satisfactory accuracy. Some buffer elements between the PML region and objects should be utilized.

Because a PML acts as an infinite open domain, all boundary conditions and material properties need to be carried over into the PML region. The pressure must be set to zero (input as PRES on the D command) on the exterior surface of the PML, excluding Neumann boundaries (symmetric planes in the model). All types of excitation sources are prohibited in the PML region.

In transient analysis, the additional auxiliary variables are introduced in Equation 8–67, referred to by Assi [442].

Refer to Artificially Matched Layers in the Acoustic Analysis Guide for more information about using PML.

PML for axisymmetric elements are modeled in two-dimensional acoustic elements. The PML located in a cylindrical enclosure represents the rotation of two-dimensional elements in the rotation axis. Therefore, the equations employed for axisymmetric PML are derived in cylindrical coordinates.

For axisymmetric PML, the modified operator, , is defined in cylindrical coordinates (r, θ, i) as follows:

(8–70)

where and are complex constants, and is given by:

(8–71)

where is the distance to the rotation axis, and is the distance from the interface between PML and the acoustic medium to the rotation axis. By using Galerkin method, the "weak" formulation of the acoustic pressure wave in axisymmetric PML is demonstrated as:

(8–72)

where 𝒶 = mode number.

Irregular perfectly matched layers (IPML) are located inside a convex enclosure. The IPML material is specified based on the local normal direction of the IPML interface. The pressure wave is assumed to propagate in the z'-direction of the local Cartesian coordinates (x', y', z'), which is set on the IPML interface by the program.

The wave equation for IPML with a modified differential operator is:

(8–73)

The modified operator is defined by:

(8–74)

The complex IPML parameter is constructed by:

(8–75)

Multiplying s_z by Equation 8–73 results in the final wave equation for IPML as follows:

(8–76)

where:

= diagonal tensor diag{sz, sz, 1/sz}

The "weak" formulation of an acoustic pressure wave in IPML can be obtained with the operator and complex tensor in the local Cartesian coordinate system:

(8–77)

The matrix notation for IPML is similar to Equation 8–69, except for using the local operator and the complex tensor .

The IPML material is defined using FLUID30, FLUID220, and FLUID221 elements with KEYOPT(4) = 2. The IPML region must be a convex shape. The IPML is not related to any user-defined global or local coordinate system defined.

More than one 1D IPML region may exist in a finite element model. A parabolic attenuation distribution is used to minimize numerical reflections in the IPML. Ansys, Inc. recommends using at least three IPML element layers to obtain satisfactory accuracy and the uniform thickness of the IPML region. Some buffer elements between the IPML region and objects should be utilized.

All boundary conditions and material properties need to be carried over into the IPML region. The pressure must be set to zero (input as PRES on the D command) on the exterior surface of the IPML, excluding Neumann boundaries (rigid walls in the model). All types of excitation sources are prohibited in the IPML region. The program can automatically apply zero pressure to the exterior surface of the IPML region, but the rigid walls (the Neumann boundaries) have to be flagged (with the SF,,RIGW command), if they exist in the model.

In transient analysis, the additional auxiliary variables are introduced in Equation 8–76 with the coordinate space transformation. Refer to Assi [442] for details.

Refer to Artificially Matched Layers in the Acoustic Analysis Guide for more information about using IPML.

Specified pressure (input as PRES on the D command) or normal velocity on the domain boundary offers a straightforward method to introduce the excitation source for harmonic response analysis (normal acceleration for transient analysis) (input as SHLD on the SF command) or the velocity for harmonic response analysis (acceleration for transient analysis) (input as EF on the BF command) (Equation 8–10) .

A mass source in the wave Equation 8–5 (input as MASS on the BF command) can be used to introduce excitations for a pressure wave.

A point mass source is given by:

(8–78)

where:

A line mass source is given by:

(8–79)

where:

A surface mass source is given by:

(8–80)

where:

The volumetric mass source Q in Equation 8–8 has the unit [kg/(m³s)].

Depending on how the load is applied on the finite element, the program will interpret the load as a point, line, surface, or volumetric mass source. For a transient analysis, the mass source rate () in Equation 8–8 should be defined.

The various incident wave sources are defined using the AWAVE command for harmonic response analysis.

A time-harmonic incident planar wave (input as PLAN on the AWAVE command) is given by:

(8–81)

where:

Pinc,0= magnitude of incident planar wave

φ = initial phase angle of incident planar wave (usually ignored)

The particle velocity is given by:

(8–82)

A time-harmonic radiated monopole pressure (input as MONO on the AWAVE command) is introduced by:

(8–83)

where:

k = wave number

If the normal velocity is along the direction, the velocity of particle is given by:

(8–84)

Assuming that the monopole is a pulsating sphere (r=a) with surface pressure p_a, the pressure can be expressed as:

(8–85)

If the normal velocity on the pulsating sphere surface is equal to v_a, the pressure is expressed as:

(8–86)

Other applicable incident sources (Whitaker [[405]]) are also defined by the AWAVE command.

To apply the incident wave to the model, the splitting pressure should be introduced on the exterior excitation surface when the total pressure is solved by FEM. The total pressure consists of the incident pressure p^inc and scattered pressure p^sc. This can be expressed as:

(8–87)

For a planar wave, the incident wave and scattered wave can be written respectively as:

(8–88)

(8–89)

where:

Assuming the normal scattered wave vector is in the opposite direction of the normal incident wave vector, this can be expressed as:

(8–90)

After substituting Equation 8–82 into the boundary integration of Equation 8–8, this can be expressed as:

(8–91)

where:

The damping matrix and source vector in Equation 8–37 for lossless and source-free fluid with the planar incident wave on the exterior boundary is given by:

(8–92)

(8–93)

A similar procedure can be used for a spherical incident wave, but the scattered wave is assumed to be a planar wave in the direction. The damping matrix and source vector in Equation 8–37 for lossless and source-free fluid with the spherical incident wave on the exterior boundary is given by:

(8–94)

(8–95)

If we assume that the incident wave perpendicularly projects on the transparent exterior surface of the model, the analogous derivation can be achieved with the Robin boundary condition (Equation 8–38), using the wave admittance Y₀ as the impedance boundary. Since the incoming incident wave is propagating into the model, the Robin boundary condition is expressed as:

(8–96)

where:

Assuming that the scattered wave has the same wave admittance as the incident wave, the Robin boundary condition for the outgoing scattered wave is expressed as:

(8–97)

where:

Substituting Equation 8–96 into the boundary integration of Equation 8–8 and utilizing Equation 8–97 leads to:

(8–98)

The damping matrix and source vector in Equation 8–37 for lossless and source-free fluid with the incident excitation on the impedance boundary is expressed as:

(8–99)

(8–100)

The above equations are useful for the analysis of the wave scattering and the pipe structure, respectively.

The planar scattered pressure wave propagating along the opposite direction of the incident wave is assumed to be on the exterior surface of the model, while the incident wave projects into the simulation domain from outside of the model. Obviously, this assumption leads to a numerical error once the non-planar scattered wave occurs. Hence, the PML is used to absorb the arbitrary scattered outgoing waves.

Assuming that the PML region is a scattered pressure region and that other regions are total pressure regions, the "weak" forms for the total pressure and the scattered pressure are respectively written as:

(8–101)

In PML scattered pressure region:

(8–102)

Total pressure can be split into scattered pressure and incident pressure on the interface between the general acoustic region and the PML region. This can be expressed as:

(8–103)

The degrees of freedom on nodes with PML attributes are scattered pressure. The surface integration of scattered pressure on the interface will be cancelled out in Equation 8–101 and Equation 8–102. Assuming that the inside model is source free, the general matrix notation in the normal acoustic region, including the pressure on interface nodes, can be written as:

(8–104)

Substituting Equation 8–103 into Equation 8–104 leads to the matrix equation for the interior total pressure and the interface scattered pressure, expressed as:

(8–105)

The matrix equation for whole domain is expressed as:

(8–106)

where:

= complex stiffness matrix for interior total pressure from Equation 8–105

= complex stiffness matrix on interface from Equation 8–105

= complex stiffness matrix for interface scattered pressure from Equation 8–105

= complex stiffness matrix for PML interior scattered pressure from Equation 8–69

finc = complex source vector from Equation 8–105

The various analytic modes in a guided duct are defined using the APORT command for harmonic response analysis.

A time-harmonic incident planar wave (Label = PLAN on APORT) is given by:

(8–107)

where:

= complex initial pressure

= vector wave number

= position vector

A mode existing in a rectangular duct (Label = RECT on APORT) is given by:

(8–108)

where:

= width of duct

= height of duct

= propagating constant in z-direction

A mode existing in a circular duct (Label = CIRC on APORT) is given by:

(8–109)

where:

= mth order Bessel function of the first kind

= nth zero of the derivative of the mth order Bessel function of the first kind

= radius of the circular duct

A mode existing in a coaxial duct (Label = COAX on APORT) is given by:

(8–110)

where:

= coefficients determined by pressure and rigid wall boundary condition at the inner radius r = a

= mth order Bessel function of the first kind

= mth order Bessel function of the second kind

= nth zero of the equation

= inner radius of the coaxial duct

= outer radius of the coaxial duct

A power source in the wave equation (Equation 8–239) (input as MASS on the BF command) can be used to introduce excitations for a diffusion equation solution.

A point power source is given by:

(8–111)

where:

A line power source is given by:

(8–112)

where:

A surface power source is given by:

(8–113)

where:

The volumetric power source P in Equation 8–239 has units of [W/m³].

In non-uniform acoustic media the mass density and sound speed vary with the spatial position. Wave Equation 8–5 in lossless media is rewritten as:

(8–114)

According to the ideal gas law, the equation of state and the speed of sound in an ideal gas are respectively expressed as:

(8–115)

(8–116)

where:

CP = heat coefficient at a constant pressure per unit mass (input as C on the MP command)

CV = heat coefficient at a constant volume per unit mass (input as CVH on the MP command)

R = gas constant

T = temperature

Pstate = the absolute pressure of the gas measured in atmospheres

Assuming density ρ₀ (input as DENS on the MP command) and sound speed c₀ (input as SONC on the MP command) at the reference temperature T₀ (input on TREF command), and the reference static pressure (input as P_sref on R command) the density and sound speed in media are cast as follows:

(8–117)

(8–118)

The spatial temperature is input as TEMP on BF command and the nodal static pressure is input as SPRE on the BF command.

Analogous to Equation 8–38, the boundary condition across an impedance sheet in the acoustic domain is given by:

(8–139)

where:

Assuming that the acoustic domain is divided into Ω⁺ and Ω^- with the impedance sheet on the surface, substituting Equation 8–139 into Equation 8–10 casts the damping matrix as Equation 8–40. For a thin perforated layer, the sheet impedance (input as IMPD on the BF command) is given by:

(8–140)

where:

ZL = load impedance

d = thickness of perforated layer

= sound impedance of the perforated material defined by Equation 8–123

= complex wave number of the perforated material (ceff is defined by Equation 8–122)

For a thin perforated material layer backed with a rigid wall, the impedance is given by:

(8–141)

An acoustic waves propagating in viscous-thermal media have a complex propagating constant in the frequency domain. The attenuation of the acoustic wave is proportional to the shear and bulk viscosity and the thermal conduction coefficient of the media.

 

8.3.5.4.1. Boundary Layer Impedance (BLI) Model

In a rigid-walled, bounded domain, the interaction between he acoustic movement and both the diffusion of heat and the diffusion of the shear wave leads to a reactive and absorbing process on the walls. The entropic and vortical perturbations diffuse into the rigid wall in a direction normal to the boundary and decay to zero depending on the thickness of the boundary layers. A hybrid numerical and analytic solution using the boundary layer impedance (BLI) model is proposed with the following assumptions ([411]):

• The thicknesses of the boundary layers are much smaller than the dimensions of the numerical domain, while the pressure variation is constant over the boundary layers’ thicknesses.

• The acoustic pressure is approximated by its zero expansion in the boundary layers and its flow is tangential to the wall.

• Spatial variations of both velocity and temperature in the normal direction are much greater than spatial variations in the tangential directions.

The fundamental Helmholtz governing equation and boundary condition in the BLI model is given by:

(8–142)

where:

k0 = ω / c0

lvh = lv+(γ-1)lh

kw = the transverse wavenumber

l’v = μ / (ρ0c0)

lh = λ / (ρ0c0Cp)

γ = Cp / Cv = the specific heat ratio

ω = the angular frequency

ρ0 = the mass density of the fluid

c0 = the sound speed of the fluid

μ = the dynamic viscosity

η = the bulk viscosity

λ = the thermal conductivity

Cp = the heat coefficient at constant pressure per unit mass

Cv = the heat coefficient at constant volume per unit mass

 

8.3.5.4.2. Low Reduced Frequency (LRF) Model

The low reduced frequency (LRF) model is derived from simplified Navier-Stokes equations with the following assumptions ([412]):

• The acoustic wavelength is much greater than the length scale of the geometry

• The acoustic wavelength is much greater than the boundary layer thickness

The Helmholtz wave equation with modified density and bulk module for the LRF model is as follows:

(8–143)

where:

ρ0 = the mass density

K0 = the modified bulk module

l = the length scale

σ = Cpη / κ = Prandtl number

γ = specific heat ratio

The viscous and thermal effects are taken into account by the function B(s/l) and n(sσ/l), respectively.

The analytic solution of both the functions B(s/l) and n(sσ/l) can be derived for three particular cases: thin layer, tube with rectangular cross-section, and tube with circular cross section.

A thin layer is specified by thickness h. The functions B(s/l) and n(sσ/l) are derived by:

(8–144)

(8–145)

where:

A tube with a rectangular cross-section is specified by width w and height h. The functions B(s/l) and n(sσ/l) are derived by:

(8–146)

(8–147)

where:

A tube with a circular cross section is specified by radius R. The functions B(s/l) and n(sσ/l) are derived by:

(8–148)

(8–149)

where:

J0 = the zero-order Bessel function

J2 = the second-order Bessel function