COMSOL 6.4 - Flow Duct — With Boundary Mode Analysis

The acoustic field in a model of an axially symmetric lined aero-engine duct, based on modal sound transmission, is analyzed. The source is generated by a single mode excitation at a boundary. Sources and nonreflecting conditions are applied using port boundary conditions. The model analysis is performed in three steps: first computing the background mean flow (compressible irrotational potential flow), then analyzing the propagating modes with a boundary mode analysis, and finally solving the acoustic field in the lined flow-duct with the linearized potential flow equations.

This model represents an extension of the Flow Duct model where the modes used at the ports are computed numerically using the Linearized Potential Flow, Boundary Mode interface. Here assuming the hard wall conditions for the modes. The model illustrates how to reference the computed modes that result from a Mode Analysis study, and use them in the Port boundary conditions. Yet another variant of this tutorial exists, the Flow Duct — Modes with Impedance Condition, where the modes used at the ports are computed including the wall lining (impedance condition).

Model Definition

The 2D axisymmetric duct geometry representation used in this model, shown in Figure 1, is taken from Ref. 1. It is an approximate model of the inlet section of a turbofan engine in the very common CFM56 series.

Figure 1: The duct geometry including reference planes used in the model.

The spinner and duct-wall profiles are given, respectively, by the equations

where 0 ≤ z′ = z/L ≤ 1, and L = 1.86393 is the duct length. A noise source is imposed at z′ = 0, henceforth referred to as the source plane. This is where the fan would be located in the actual engine geometry. The plane z = L corresponds to the fore end of the engine and is referred to as the inlet plane. The attenuation of the liner for specific flow conditions is computed from the source plane to the inlet plane. A cylindrical domain, adjoined at the inlet plane and extending to the terminal plane, extends the modeling domain into a region where you can consider the mean flow as being uniform. This allows you to impose the simple boundary condition of a constant velocity potential and a vanishing tangential velocity for the background flow. For the acoustic problem, port boundary conditions are used at the source plane and the terminal plane to set up ideal nonreflecting conditions as well as imposing the source.

The model will analyze so-called modal sound transmission, where a single propagating mode is used as source. In this particular example the first radial mode is used as the source, see Ref. 1 for details. All propagating modes are used when setting up the ports to ensure good nonreflecting performance. The sound transmission loss is computed from the source plane to the inlet plane. The power of the incident mode is defined through a predefined variable and the power of the transmitted sound at the inlet plane is computed as the integral of the axial intensity.

Model Conditions

Assume that the flow in the axisymmetric duct is compressible, inviscid, perfectly isentropic, and irrotational. This is an assumption often used for the study of duct or engine acoustics. In this case the background mean flow is well described by the Compressible Potential Flow interface and the acoustic field is well described by the Linearized Potential Flow, Frequency Domain interface.

For more theory information on the governing equations, see the aeroacoustics theory chapter in the Acoustics Module User’s Guide.

Compressible Potential Flow

This study examines two cases for the mean-flow normal velocity component at the source plane Vz, which (owing to the choice of reference speed) alternatively can be referred to as the source-plane axial Mach number M = −0.5, approximately representative of a passenger aircraft at cruising speed, and M = 0.

The governing equations are nondimensionalized in the present study. For the reference quantities in this model, choose the duct radius, the mean-flow speed of sound, and the mean-flow density at the source plane. Hence, all three of these quantities take the value 1.

The remaining boundary conditions for the mean flow consist of a natural boundary condition specifying the mass-flow rate through the source plane via the normal velocity and the density; slip conditions (vanishing tangential velocity) at the duct wall and at the spinner; and axial symmetry at r = 0.

Linearized Potential Flow

For the aeroacoustic field, the model considers two different boundary conditions at the duct wall:

•

Sound hard — the normal component of the acoustic particle velocity vanishes at the boundary.

•

Impedance — the normal component of the acoustic particle velocity is related to the particle displacement through the equation

where Z is the impedance, u0 is the mean background flow, p is the acoustic pressure, and u is the acoustic velocity. This condition is often referred to as the Ingard–Myers impedance condition. This boundary condition, first derived by Myers (Ref. 2), was later recast in a weak form by Eversman (Ref. 3); it is this weak version, which is directly suitable for finite element modeling, that is implemented in the Acoustics Module’s Linearized Potential Flow, Frequency Domain interface. The impedance boundary condition represents a lined duct wall. In this model, following Ref. 1, the impedance is taken to be Z = 2 − i.

The spinner, in contrast, is always assumed to be acoustically hard.

One of the configurations from Ref. 1 is studied in this model. This is the case where the dimensionless angular frequency (nondimensionalized through division by R∞/c∞) is ω = 16, and the azimuthal mode number is m = 10. If you want to obtain a deeper understanding of the duct’s aeroacoustic characteristics, you can, of course, perform a systematic exploration of parameter space by varying these quantities independently. Several more cases are examined in the reference paper.

Results and Discussion

The Mean-Flow Field

For the nontrivial case of a source-plane axial Mach number of M = −0.5, the resulting mean-flow field appears in Figure 2. Note that the velocity potential is uniform well beyond the terminal plane, thus justifying the boundary condition imposed there. Furthermore, as could be expected, deviations from the mean density value appear primarily near the nonuniformities of the duct geometry, such as at the tip of the spinner.

As a complement, a more quantitative picture of the variations of the mean-flow velocity and density profiles along the axial direction (for r = 0.8) appear in the cross-section plots in Figure 3.

Figure 2: Mean-flow velocity potential and density for source-plane Mach number M = −0.5.

Figure 3: Mean-flow cross section plot at a sample radius of 0.8.

Source-Plane Modes

With the solution for the mean-flow field at hand, it is possible to calculate the corresponding eigenmodes for the acoustic field at the ports. To use the modes at the ports it is necessary to know if they are incident or outgoing. In the no-flow case this is simply done by looking at the sign of the wave vector. The outgoing mode has a positive wave number and the incident is negative. However, in the flow case this is not as straightforward because the presence of the mean background flow shifts the eigenvalues. Here it can be advantageous to look at the sign of the intensity vector at the port to identify the propagation direction.

A plot of the modes at the source plane (z = 0) including the orientation of the axial intensity vector Iz is depicted in Figure 4 for the flow case (M = −0.5). The source plane is located below the duct, such that an incident mode has a positive intensity and vice versa. Note that the ordering of the modes is done based on the settings in the section of the study step. In this case the modes are sorted by ascending real part of the wave number. Filtering of the modes is also possible as done for the analysis of the terminal plane.

Figure 4: Mode pressures and axial intensity at the source plane (z = 0) for the case of a background flow with Mach number M = −0.5.

Figure 5 shows the resulting velocity-potential profile for all propagating modes at the source plane in the no-flow case (M = 0); and Figure 6 shows all the propagating modes in the flow case (M = −0.5). Notice that the mode profiles are not identical for the incident and outgoing modes (zoom in the graph in the model to see this clearly), indicating that the mean background flow is not uniform at this boundary.

Combining the information from all the plots, at the source plane, the mode identification necessary for the ports conditions can be done. In all cases the source is the first radial mode as defined in Ref. 1. See also the last two Evaluation Groups in the model where the solution index is shown for the various modes and flow conditions.

•

At the source plane in the no-flow case (M = 0), the outgoing mode has the wave number kn = 10.8. This is the second eigenvalue in solution list and thus as index 2 (necessary for referencing at the port). The incident mode (the source) has wave number kn = −10.8 (index 1).

•

At the source plane in the flow case (M = −0.5), the outgoing modes have wave numbers are kn = −0.9 (index 3) and kn = 5.8 (index 4). The first incident radial mode (the source) has wave number kn = −27.1 (index 1). In the port condition it is advantageous for clarity (but not necessary) that modes corresponding to the same mode shape should be defined together. This means that the first radial modes (index 1 and 4) are referenced in one port, and the second outgoing radial mode (index 3) is referenced in the second port (where no source is added).

Figure 5: Propagating modes at the source plane for the no-flow case (M = 0).

Figure 6: Propagating modes at the source plane for the flow case (M = −0.5).

Terminal-Plane Modes

A plot of the outgoing modes at the terminal plane (z = 1.86393) including the orientation of the axial intensity vector Iz is depicted in Figure 7 for the flow case (M = −0.5). The terminal plane is located above the duct, such that an incident mode has a negative intensity and vice versa. For this analysis filtering of the modes is defined such that only outgoing modes are stored in the output. The logic expression comp1.intop_tp(comp1.lpfbm2.Iz) is used in the section of the second study step. That is, only modes where the integral of the intensity is larger than 0 are stored.

Figure 8 shows the resulting velocity-potential profile for all propagating modes at the terminal plane in the no-flow case (M = 0); and Figure 9 shows all the propagating modes in the flow case (M = −0.5). Notice that the mode profiles are identical for the incident and outgoing mode, indicating that the mean background flow is uniform at this boundary (as expected).

Combining the information from all the plots, at the terminal plane, the mode identification necessary for the ports conditions can be done. See also the last two Evaluation Groups in the model where the solution index is shown for the various modes and flow conditions. No source is defined at the terminal plane.

•

At the terminal plane in the no-flow case (M = 0), the outgoing mode has the wave number kn = 9.6. This is the first and only eigenvalue in solution list and has solution index 1 (necessary for referencing at the port). No incident mode is necessary as no source is defined here.

•

At the terminal plane in the flow case (M = −0.5), the outgoing modes have wave numbers kn = 13.9 (index 1) and kn = 24.9 (index 2). No incident mode is necessary as no source is defined here.

Figure 7: Outgoing mode pressures and axial intensity at the terminal plane for the case of a background flow with Mach number M = −0.5.

Figure 8: Propagating modes at the terminal plane for the no-flow case (M = 0).

Figure 9: Propagating modes at the terminal plane for the flow case (M = −0.5).

The Aeroacoustic Field

The normalized pressure fields for the case without a no background mean flow (M = 0), shown in Figure 10, very closely match those for the corresponding finite element model (FEM) solutions presented in Figure 6 of Ref. 1. Similarly, the results for the attenuation between the source and inlet planes in the lined-wall case are in good agreement: 50.67 dB for the COMSOL Multiphysics solution versus 51.6 dB for the FEM solution, as shown in Table. 1 in Ref. 1.

Turning to the case with a mean flow (M = −0.5), the pressure field for the hard-wall as well as the lined wall (soft wall) cases in Figure 11 closely resembles the FEM solution obtained by Rienstra and Eversman in Ref. 1. This observation extends to the attenuation, for which the calculated value of 28.0 dB is in good agreement with the value of 27.2 dB obtained in Ref. 1.

Note that the source mode in the COMSOL Multiphysics calculation was derived for the case of a hard duct wall, whereas Rienstra and Eversman used a noise source adapted to the acoustic lining. However this fact does not seem to have a large influence on the solution for this particular problem. The propagating mode for the lined wall is actually a linear combination of the two hard-wall propagating modes.

Figure 10: Acoustic pressure field for the cases of hard (top) and lined (bottom) duct wall with no mean flow (M = 0); azimuthal mode number m = 10 and angular frequency ω = 16.

Figure 11: Acoustic pressure distribution for the cases of hard (top) and lined (bottom) duct wall with mean flow (M = −0.5); azimuthal mode number m = 10 and angular frequency ω = 16.

Notes About the COMSOL Implementation

Physics Interfaces

•

Compressible Potential Flow (cpf) — for modeling the background mean-flow velocity field as a potential flow (a lossless and irrotational flow).

•

Linearized Potential Flow, Boundary Mode (lpfbm) — for calculating the boundary eigenmode to be used by the port boundary conditions, defining outgoing and incident (the source) propagating acoustic field.

•

Linearized Potential Flow, Frequency Domain (lpff) — for modeling the time-harmonic acoustic field in the duct for the various excitation and flow condition.

Referencing the Acoustic Modes at the Ports

The modes computed at the source plane and the terminal plane with the boundary mode interfaces need to be referenced and used by the port conditions. This is achieved by using the withsol() operator. The operator is called with a solution tag, referencing which solution it should look at, here the tag used depends on if the model is solved at the source ('sol3') or the terminal planes ('sol7'). The input is the variable needed, for example, phi_sp for the source plane potential or lpfbm.kn for the source plane mode wave-number. Finally, the operator is called with two arguments using setval() to reference if the Mach number M is 0 or -0.5; and using setind() to set the solution index of the mode (of the eigenvalue lambda). For the index it is the number of the eigenvalue in the solution object. Note that lambda is always used as the internal eigenvalue variable in COMSOL Multiphysics.

Filtering and Sorting of Eigenmodes

It can be advantageous to control the sorting of the results of a mode analysis study. It can also be advantageous to control what modes are stored (filtered) in the output. Both can be achieved with the settings of the section of the study step.

References

1. S.W. Rienstra and W. Eversman, “A Numerical Comparison Between the Multiple-Scales and Finite-Element Solution for Sound Propagation in Lined Flow Ducts,” J. Fluid Mech., vol. 437, pp. 367–384, 2001.

2. M.K. Myers, “On the Acoustic Boundary Condition in the Presence of Flow,” J. Sound Vib., vol. 71, pp. 429–434, 1980.

3. W. Eversman, “The Boundary Condition at an Impedance Wall in a Non-Uniform Duct with Potential Mean Flow,” J. Sound Vib., vol. 246, pp. 63–69, 2001. Errata: ibid, vol. 258, pp. 791–792, 2002.

Acoustics_Module/Aeroacoustics_and_Noise/flow_duct_boundary_mode

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

In the tree, select > > .

Click .

In the text field, type phi_sp.

Here _sp stands for source plane.

In the tree, select > > .

Click .

In the text field, type phi_tp.

Here _tp stands for terminal plane.

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Root

In the window, click the root node.

In the root node’s window, locate the section.

From the list, choose .

This setting turns off all unit support in the model.

Global Definitions

Parameters 1

Load the parameters from a file. They define model, geometry, and physical properties including the liner impedance. Then proceed and create the geometry of the duct.

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file flow_duct_boundary_mode_parameters.txt.

Proceed and draw the geometry of the engine duct. Use the features to draw the shapes defined by the functions described in the main document.

Geometry 1

Parametric Curve 1 (pc1)

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type 1-0.18453*s^2+0.10158*(exp(-11*(1-s))-exp(-11))/(1-exp(-11)).

In the text field, type s*zi.

Parametric Curve 2 (pc2)

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type 0.7.

Locate the section. In the text field, type 0.64212-sqrt(0.04777+0.98234*s^2).

In the text field, type s*zi.

Line Segment 1 (ls1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the text field, type zi.

Union 1 (uni1)

In the toolbar, click

and choose .

Select the objects and only.

Line Segment 2 (ls2)

In the toolbar, click

and choose .

In the window for , locate the section.

Click to select the

toggle button for .

On the object , select Point 5 only.

Locate the section. Click to select the

toggle button for .

On the object , select Point 1 only.

Line Segment 3 (ls3)

In the toolbar, click

and choose .

On the object , select Point 4 only.

In the window for , locate the section.

Click to select the

toggle button for .

On the object , select Point 2 only.

In the toolbar, click

Delete Entities 1 (del1)

In the window, right-click and choose .

On the object , select Boundaries 1 and 3 only.

In the window for , click

Convert to Solid 1 (csol1)

In the toolbar, click

and choose .

Click in the window and then press Ctrl+D to clear all objects.

Click the

button in the toolbar.

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type ri.

Locate the section. In the text field, type zi.

Form Union (fin)

In the toolbar, click

Click the

button in the toolbar.

In the window, click .

Proceed and set up variables used for the results analysis. One is a normalized absolute pressure which uses a maximum operator over the domain. Define selections for the source, inlet and terminal planes. Finally, define an integration operator used to compute the power through the inlet plane.

Definitions

Variables 1

In the window, under right-click and choose .

In the window for , locate the section.

From the list, choose .

Select Domain 1 only.

Locate the section. In the table, enter the following settings:


Name	Expression	Unit	Description
Mz	-cpf.Vz		Axial Mach number
pabsn	abs(lpff.p)/comp1.maxop1(lpff.p)		Normalized pressure

Source Plane

In the toolbar, click

In the window for , type Source Plane in the text field.

Locate the section. From the list, choose .

Select Boundary 5 only.

Inlet Plane

In the toolbar, click

In the window for , type Inlet Plane in the text field.

Locate the section. From the list, choose .

Select Boundary 3 only.

Terminal Plane

In the toolbar, click

In the window for , type Terminal Plane in the text field.

Locate the section. From the list, choose .

Select Boundary 4 only.

Maximum 1 (maxop1)

In the toolbar, click

and choose .

Select Domain 1 only.

Integration 1 (intop1)

In the toolbar, click

and choose .

In the window for , type intop_ip in the text field.

Locate the section. From the list, choose .

From the list, choose .

Integration 2 (intop2)

In the toolbar, click

and choose .

In the window for , type intop_tp in the text field.

Locate the section. From the list, choose .

From the list, choose .

Now proceed and set up the physics for the Compressible Potential Flow as well as the two Boundary Mode interfaces. The latter two are used to compute the propagating modes at the source and the terminal planes. The modes are used in the Port boundary conditions used for the frequency domain analysis.

Compressible Potential Flow (cpf)

In the window, under click .

In the window for , locate the section.

In the pref text field, type cpf.rhoref^gamma/gamma.

In the ρref text field, type rho0.

In the vref text field, type M.

Compressible Potential Flow Model 1

In the window, under > click .

In the window for , locate the section.

From the γ list, choose . In the associated text field, type gamma.

Normal Flow 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Mass Flow 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Linearized Potential Flow, Boundary Mode (lpfbm)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Locate the section. In the m text field, type m.

Linearized Potential Flow Model 1

In the window, under > click .

In the window for , locate the section.

From the u0 list, choose .

Locate the section. From the ρ0 list, choose .

From the c0 list, choose .

Notice the orientation of the outward propagating normal, the red arrow seen in the window. For consistency it should point out of the domains. It is possible to reverse the orientation using the checkbox. This is not necessary here.

Linearized Potential Flow, Boundary Mode 2 (lpfbm2)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Locate the section. In the m text field, type m.

In the window, under > click .

In the window for , locate the section.

From the u0 list, choose .

Locate the section. From the ρ0 list, choose .

From the c0 list, choose .

Set up a fully user defined mesh for the computational domain.

Mesh 1

Free Triangular 1

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Domain 1 only.

Size 1

Right-click and choose .

In the window for , locate the section.

Click the button.

Locate the section. From the list, choose .

Select Points 4 and 7 only.

Locate the section.

Select the checkbox. In the associated text field, type 0.005.

Size

In the window, under > click .

In the window for , locate the section.

Click the button.

Locate the section. In the text field, type 0.015.

In the text field, type 0.001.

Mapped 1

In the toolbar, click

In the window for , click

Now, first solve the background flow and look at the results. The flow is solved for two Mach numbers using a Parametric Sweep.

Study 1 - Background Flow

In the window, click .

In the window for , type Study 1 - Background Flow in the text field.

Parametric Sweep

In the toolbar, click

In the window for , locate the section.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
M (Mean flow Mach number)	0 -0.5

In the toolbar, click

Results

Arrow Surface 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

In the toolbar, click

Cut Line 2D 1

In the toolbar, click

In the window for , locate the section.

In row , set to 0.8.

In row , set to zi.

Mean Flow: rho and Mz

In the toolbar, click

In the window for , type Mean Flow: rho and Mz in the text field.

Locate the section. From the list, choose .

From the list, choose .

Click to expand the section. From the list, choose .

Locate the section. From the list, choose .

Line Graph 1

Right-click and choose .

In the window for , locate the section.

In the text field, type rho.

Locate the section. From the list, choose .

In the text field, type z.

Click to expand the section. Select the checkbox.

Find the subsection. Select the checkbox.

Find the subsection. In the text field, type Ma = .

In the toolbar, click

Line Graph 2

In the window, right-click and choose .

In the window for , locate the section.

In the text field, type Mz.

Locate the section. From the list, choose .

In the text field, type z.

Locate the section. Select the checkbox.

Find the subsection. Select the checkbox.

Find the subsection. In the text field, type Ma = .

In the toolbar, click

Proceed to compute and analyze the mode shapes for the source plane. The mode search method will be used as it can filter out all the evanescent (nonpropagating) modes by setting a small range for the imaginary part of the out-of-plane wave number. The real part lies between -k0max_abs and +k0max_abs (the maximal absolute wave number) computed in the parameters list. Note that an additional small margin is used for the interval. A parametric sweep is also performed solving for the two Mach numbers.

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the table, clear the checkboxes for , , and .

Find the subsection. In the tree, select > .

Click the button in the window toolbar.

Study 2 - Source Plane Modes

In the window for , type Study 2 - Source Plane Modes in the text field.

Parametric Sweep

In the toolbar, click

In the window for , locate the section.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
M (Mean flow Mach number)	0 -0.5

Step 1: Mode Analysis

In the window, click .

In the window for , locate the section.

In the text field, type f.

From the list, choose .

In the text field, type 5.

In the text field, type 10.

Find the subsection. In the text field, type -1.1*k0max_abs.

In the text field, type 1.1*k0max_abs.

In the text field, type -0.1.

In the text field, type 0.1.

Click to expand the section. Find the subsection. From the list, choose .

From the list, choose .

Notice in the window the section called . This section contains settings for filtering the eigenmodes, as well as settings for sorting (or ordering) of the eigenmodes in the output. At the source plane both outgoing and incident modes are needed, so no filtering is performed. At the terminal plane however, only outgoing modes are needed (no source is present). The setup of the next study shows how such a filtering can be easily set up.

In the toolbar, click

Results

Source Plane: Acoustic Pressure and Axial Intensity

When setting up the Port boundary conditions it is necessary to know which modes are outgoing and incident. Simply knowing the sign of the computed out-of-plane wave number is not enough when a background flow is present and when higher order azimuthal modes are analyzed. In this case it is advantageous to plot the axial intensity vector to identify the propagation directions. This is done by modifying the next default plot, here at the source plane.

In the window for , type Source Plane: Acoustic Pressure and Axial Intensity in the text field.

Click to expand the section. From the list, choose .

From the list, choose .

Select the checkbox.

Click to expand the section. From the list, choose .

Arrow Line 1

Right-click and choose .

In the window for , click in the upper-right corner of the section. From the menu, choose > > > .

Locate the section. In the text field, type 0.

Locate the section. In the text field, type 30.

Locate the section. From the list, choose .

Source Plane: Acoustic Pressure and Axial Intensity

In the window, click .

In the toolbar, click

Click the

button in the toolbar.

In the window for , locate the section.

From the list, choose .

In the toolbar, click

From the list, choose .

In the toolbar, click

From the list, choose .

In the toolbar, click

Looking at the four out-of-plane wave numbers (for the M = -0.5 case) it can be seen that the first two are inward propagating and the last two are outward propagating.

Revolution 2D 1

In the window, under > click .

In the window for , click to expand the section.

From the list, choose .

The option can be used when evaluating the dependent variable (the velocity potential). However, here we will evaluate the pressure. The azimuthal component will be added manually in the plot, using the defined phi variable.

Source Plane: Acoustic Pressure, 3D (lpfbm)

In the window, under click .