http://www.diva-portal.org
This is the published version of a paper presented at ISMA2018-USD2018.
Citation for the original published paper:
Cuenca, J., Göransson, P., de Ryck, L., Lähivaara, T. (2018)
Inverse parameter estimation in resonant, coupled fluidstructure interaction problems In:
N.B. When citing this work, cite the original published paper.
Permanent link to this version:
http://urn.kb.se/resolve?urn=urn:nbn:se:kth:diva-240140
Inverse parameter estimation in resonant, coupled fluid- structure interaction problems
J. Cuenca1,2, P. G¨oransson2, L. De Ryck1, T. L¨ahivaara3
1Siemens Industry Software, Interleuvenlaan 68, B-3001 Leuven, Belgium email: jacques.cuenca@siemens.com
2 Aeronautical and Vehicle Engineering, KTH Royal Institute of Technology, Teknikringen 8, SE-10044 Stockholm, Sweden
3 Department of Applied Physics, University of Eastern Finland, Yliopistonranta 1, FIN-70211 Kuopio, Finland
Abstract
This work discusses the estimation of model parameters in resonant coupled systems over a wide frequency range. Such problems are known to be subjected to local minima, which represent a major obstacle in the field of parameter identification. This work evaluates various approaches to finding the global minimum, i.e. the true model parameters. The L2 norm between a target and a model frequency response function is used as the objective function. The estimation is performed through deterministic and statistical inversion frameworks, using both full-spectrum and step-wise approaches. The occurrence and the background of local minima are discussed and the performance of the different solution methods is evaluated. In particular, the correlation between different model parameters is analysed and observed to control the path towards the solution. Two illustrative examples are proposed, including the estimation of geometrical and material properties of an expansion chamber coupled to a limp porous material.
1 Introduction
Material and geometrical properties are observable individually by using direct methods. However, direct measurement requires dedicated setups for each individual parameter and such measurements are rarely practicable on a single material sample. An example of this is the case of porous media, where several separate setups are required for a direct observation of all acoustic and mechanical properties, e.g. [1–3].
In addition and most importantly, due to the mutual independence of individual parameter observations, the estimated material parameters are not guaranteed to satisfy a constitutive model.
As an alternative, inverse estimation methods rely on the ability for a model to describe an experimental observation. The material and geometrical properties are obtained indirectly as model parameters, as a result of a curve-fitting procedure. This guarantees that any set of trial parameters satisfies the constitutive model.
Typically the inversion procedure relies on optimisation algorithms, ensuring a low number of model evaluations. The most commonly used inversion frameworks rely on a gradient-based minimisation of an objective function. One of the major obstacles in such approaches is the presence of local minima, which can arise as a consequence of multiple factors. Recent work has shown that anisotropy in material properties leads to similarities in the behaviour of a structure for specific spatial orientations, which in turn results in local minima [4, 5]. In the case of inverse estimation of material or geometrical properties using frequency- dependent data, a major source of local minima is the presence of multiple resonances. Indeed, the inherent similarity among multiple resonances is a cause of wrong parameter estimations.
Other inversion paradigms do not suffer from local minima in the same manner, but generally require a large number of model evaluations. For instance, Bayesian inversion framework enables to investigate the inverse problem in a probabilistic sense. However, a significant number of model evaluations is required in order to provide a reliable estimate of the probability density function [6].
Effectively a trade-off exists between a low number of model evaluations and the assurance of finding the global minimum of the problem. Strategies to overcome the problem of multiple local minima in inverse estimation problems include, for instance, the use of multiple starting points with randomly distributed values [7], the use of neural networks or related paradigms [8], or the simultaneous use of different optimisation algorithms concurrently exploiting common evaluation points.
The present paper addresses the problem of inverse estimation of material and geometrical properties of structures exhibiting a non-monotonous behaviour with frequency. A strategy is here employed in order to overcome resulting local minima, based on physical principles governing the dynamic behaviour. The method, proposed in a recent paper [9], decomposes the broadband inverse estimation problem into a number of subsequent subproblems over increasingly broad frequency ranges. In the case of a coupled fluid-structure problem, a starting subproblem over a frequency band where the structure exhibits an asymptotic behaviour governed by a partial set of the unknowns has been shown to guide the solution towards the true solution.
In the present paper, the method is presented and applied to problems of reduced complexity, thus allowing for analytical derivations to support the observations. A gradient-based framework is used as the underlying inversion method in the proposed method. This is used in combination with a Bayesian framework, which provides an insight into the nature of the multiple local minima, including the true solution of the problem.
The combined use of deterministic and statistical approaches thus allows to benefit from a low number of model evaluations while gaining additional knowledge on the problem. The focus of the present paper is nonetheless on the physical aspects governing local minima rather than on the computational efficiency itself.
2 Inverse estimation problem
2.1 Observation model
This paragraph describes the general aspects of the inverse estimation problem, without a particular focus on a specific physical system. As described in the introduction, the goal is to measure material and/or geometrical properties of frequency-dependent systems. Accordingly, the system under consideration is characterised by an observation
g0(ωm), m ∈ [1, M ], (1)
where ωmis a sequence of discrete circular frequencies at which the observation is performed.
It is assumed that a known model describes the behaviour of the system as a function of its material and geometrical parameters, in the additive form
g0(ωm) = g(x0, ωm) + e, (2)
where x0is the set of unknown model parameters, considered as the true properties of the observed system, and e is a model of the measurement error (assumed to be Gaussian).
The primary goal is to find an estimate ˜x0of x0by means of a suitable inversion framework for the available model. The focus of the present paper is on dynamic problems and in particular those exhibiting resonances which frequently lead to multiple wrong estimations.
2.2 Gradient-based inversion
The objective function is here defined as the 2-norm distance between the modelled response and the observed response, which can be written as
d(x, M ) =
M
X
m=1
|g(x, ωm) − g0(ωm)|2. (3)
An estimate ˜x0of the unknown model parameters x0is thus obtained by minimising the objective function, in the form
˜
x0 = argmin (d(x, M )) . (4)
In the present paper, the Globally Convergent Method of Moving Asymptotes (GCMMA) [10] is used.
In addition to the global minimum of Eq. (4), resonant problems often exhibit a series of local minima, which arise as a consequence of the similarity of the system’s behaviour for different values of its properties.
A direct inversion thus proves to be ineffective, as multiple starting points will lead to multiple solutions.
Furthermore, in a general sense it is not guaranteed that the value of the objective function is strictly lower at the true solution, as shown in the first example below.
2.3 Bayesian inversion framework
In the Bayesian framework [11, 12], all quantities are modelled as random variables. The Bayesian inference is based on the Bayes’ rule and results in the conditional probability distribution of the unknowns x given the measurements g0. This distribution is commonly called as the posterior probability distribution of the unknowns x and can further be used for computing not only point estimates for x but also spread estimates, such as credible intervals.
For the sampling of the posterior densities, we use the Metropolis-Hastings algorithm with an adaptive proposal distribution scheme [13, 14], referred to as MCMC (Markov chain Monte Carlo) in the following sections. The point estimate for the parameters we use in this work is the maximum a posteriori and as an indication of the uncertainty in the point estimates we use the narrowest 95% credible interval. For more detailed discussion of the Bayesian framework used in this work, we refer to [9].
2.4 Incremental-frequency inverse estimation
The method proposed in this paragraph consists of solving the inversion problem as a sequence of sub- problems where the solution to a given sub-problem is used as the starting point for the next sub-problem.
By doing so, the aim is to guide the model parameters towards the true solution, without terminating at a local minimum.
The sequence of sub-problems is chosen following an a priori behaviour of dynamic systems, and as such is potentially suitable for a large number of problems. An initial sub-problem is defined by constraining Eq. (4) to a narrow frequency range in which the system can be expected to behave in an asymptotic manner, that is, with a dependence on a subset of the model variables and away from resonant behaviour. The solution to such sub-problem is used as the starting point for a second sub-problem, whose frequency range is wider.
The process is continued until the frequency range is that of the original full problem. The strength of such a process lies in the fact that the solution to the first sub-problem lies in a subspace that contains the true solution.
The sub-problems can be formally stated as
˜
x(n)0 = argmin (d(x, M (n))) , n ∈ [1, N ], (5) where M (n) is the number of frequencies in sub-problem n. Accordingly, ˜x(N )0 = ˜x0 ' x0 is the final solution, expected to provide a good estimate of the true model parameters.
3 Application examples
3.1 General framework
The two examples considered in this section rely on the knowledge of the sound transmission loss (STL) of the studied component, for instance as measured in an impedance tube setup. A transfer matrix approach [15]
is used, providing a relation between the acoustic pressure and particle flow at the inlet and outlet in the form
pin(ω) Sinvin(ω)
= T(ω)
pout(ω) Soutvout(ω)
, (6)
where p, v and S are respectively the acoustic pressure, particle velocity and tube cross-section area, at the inlet or outlet. By assuming that the studied component is composed of N equivalent homogeneous fluid layers. its transfer matrix may be written as the product
T(ω) =
N
Y
n=1
Tn(ω). (7)
The transfer matrix of layer n can be written in the form
Tn(ω) =
cos(knln) −Zn
iSn sin(knln) iSn
Zn sin(knln) cos(knln)
, (8)
where ln, Sn, kn and Zn respectively denote the layer’s length, cross-section area, wavenumber and characteristic impedance. It is worth noting that equivalent fluid models of complex media will exhibit a dependence of the wavenumber and impedance on intrinsic material properties. The STL of the full component is given by
STL = 20 log10 1 2
T11+ T12Sout
Zout + T21Zin
Sin + T22 Zin Zout
Sout Sin
. (9)
3.2 Expansion chamber 3.2.1 Model
As a first application example, a simple expansion chamber is considered, as illustrated in Fig. 1. The measured quantity is chosen as the sound transmission loss within a plane-wave approximation, which may be written as
STL(x) = 20 log10 1 2
2 cos(kl) + r2 R2 +R2
r2
i sin(kl)
, (10)
where R is the radius of the inlet and outlet ducts and l and r are respectively the length and radius of the expansion chamber, considered as the two unknowns of the problem, denoted as
x = {l, r}. (11)
The true values of the model parameters are denoted
x0 = {l0, r0}. (12)
l
R
r
Figure 1: Expansion chamber with unknown length and radius.
The aim of the following is to retrieve arbitrarily chosen values for l0 and r0 from the knowledge of the sound transmission loss of the system. Following Eq. (2), the observation is here modelled as
STL0(ωm) = STL(x0, ωm) + e, (13)
where e is a normally distributed random variable, representing an additive noise. The objective function is defined as
d(x) =X
m
|STL(x, ωm) − STL0(ωm)|2. (14)
First, the gradient-based inversion is applied in the full frequency spectrum, defined as f ∈ [1, 1200] Hz. The target values of the chamber length and radius are chosen as l0 = 0.5 m and r0= 0.12 m and the inlet/outlet duct radius is set to R = 0.04 m.
3.2.2 Full-spectrum gradient-based inversion
This paragraph discusses the gradient-based inverse estimation of model parameters l and r. As a convergence criterion, the GCMMA algorithm is set to stop after 5 consecutive iterations where the variation in the objective function and design variable values is less than 10−4relative to the previous iteration.
Figure 2 shows a map of the objective function as a function of (l, r) and the solutions found from running the gradient-based algorithm 1000 times using a random multi-start procedure. It shall be noted that e is chosen to be different for each run, in order to quantify the spread of the solution, rather than a solution bias.
The standard deviation of the additive noise is set to 2 dB in this case.
The figure shows that the gradient-based algorithm terminates at multiple local minima, including the true solution. This allows to conclude that the solution of this inverse problem on the full available spectrum is not viable.
r(m)
l (m) 0
0.05 0.1 0.15 0.2 0.25
0 0.2 0.4 0.6 0.8 1 1.2
f ∈[1,50] Hz
r(m)
l (m) 0
0.05 0.1 0.15 0.2 0.25
0 0.2 0.4 0.6 0.8 1 1.2
f ∈[1,50] Hz
1422.07e+06
(a)
0.118 0.12 0.122
0.498 0.5 0.502
r(m)
l (m) (b)
Figure 2: (a) Map of the full frequency range objective function (color) and local minima (dots); (b) closeup around the true solution.
An additional observation can be made on this simplified problem. In fact, Fig. 2 displays two sets of local minima, distributed horizontally. This is attributed to the dependence of the sound transmission loss on the expansion chamber radius, in the form r2/R2+ R2/r2. The expression is symmetric with respect to r = R, which implies that the objective function and its local minima within r ≤ R are a mirror image of those within r ≥ R. Therefore the observation of the sound transmission loss does not allow to distinguish an expansion from a constriction, unless a constraint such as r > R is imposed onto the inversion algorithm.
3.2.3 Incremental inverse estimation
The inverse estimation of the length and radius of the expansion chamber is here solved using the proposed incremental approach. As the starting sub-problem, a low-frequency narrow band where the system is non- resonant is chosen. An asymptotic expansion of the sound transmission loss for small kl yields
STL(x) ' 10 log10 1 + r2 R2 +R2
r2
2
(kl)2 4
!
. (15)
The behaviour of the system is thus monotonous in l and symmetric in r. The roots of the objective function, Eq. (14), at the limit kl 1, can be derived as functions of l and r and expressed as
r±= R s
r20 R2 +R2
r20
l0
l
±1/2
(16) as shown in Fig. 3. Therefore, by constraining the first sub-problem to a low frequency range, the corresponding estimate of x0is expected to fall onto one of the curves.
0 0.05 0.1 0.15 0.2 0.25
0 0.2 0.4 0.6 0.8 1 1.2
r(m)
l (m)
Figure 3: Roots of the objective function in the low frequency limit for the expansion chamber problem.
The sub-problem frequency ranges are chosen as
f ∈ [1 Hz − fmax], fmax= {50, 100, 200, 400, 700, 1200} Hz. (17) Each subproblem is solved using GCMMA with identical stabilisation criteria as for the full-range estimation.
An example of solution using the proposed approach is here presented in detail. Table 1 summarises the sub-problem frequency ranges, the final estimate per sub-problem and the deviation from the true solution for an arbitrary starting point. Figure 4 shows a map of the objective function for each sub-problem, together with the path of the design variables x = (l, r). It is worth noting that similar results are obtained by using different starting points. It has been observed for the present case that the procedure invariably converges to the true solution (l0, r0), or to its mirror image (l0, R2/r0).
Sub-problem Stable estimate (l, r) Absolute relative deviation f ∈[1,50] Hz (0.2663, 0.1601) m (46.74, 33.43) % f ∈[1,100] Hz (0.7709, 0.1011) m (54.18, 15.73) % f ∈[1,200] Hz (0.5088, 0.1202) m (1.75, 0.20) % f ∈[1,400] Hz (0.5003, 0.1202) m (0.06, 0.18) % f ∈[1,700] Hz (0.4992, 0.1196) m (0.16, 0.33) % f ∈[1,1200] Hz (0.4997, 0.1203) m (0.07, 0.27) %
Table 1: Sub-problem frequency ranges, final estimates and deviation from the target values x0 = (0.5, 0.12) m.
r(m)
l (m) 0
0.05 0.1 0.15 0.2 0.25
0 0.2 0.4 0.6 0.8 1 1.2
f ∈[1,50] Hz
r(m)
l (m) 0
0.05 0.1 0.15 0.2 0.25
0 0.2 0.4 0.6 0.8 1 1.2
f ∈[1,50] Hz
0.01432.3e+06
r(m)
l (m) 0
0.05 0.1 0.15 0.2 0.25
0 0.2 0.4 0.6 0.8 1 1.2
f ∈[1,100] Hz
r(m)
l (m) 0
0.05 0.1 0.15 0.2 0.25
0 0.2 0.4 0.6 0.8 1 1.2
f ∈[1,100] Hz
0.3242.3e+06
r(m)
l (m) 0
0.05 0.1 0.15 0.2 0.25
0 0.2 0.4 0.6 0.8 1 1.2
f ∈[1,200] Hz
r(m)
l (m) 0
0.05 0.1 0.15 0.2 0.25
0 0.2 0.4 0.6 0.8 1 1.2
f ∈[1,200] Hz
0.5421.99e+06
r(m)
l (m) 0
0.05 0.1 0.15 0.2 0.25
0 0.2 0.4 0.6 0.8 1 1.2
f ∈[1,400] Hz
r(m)
l (m) 0
0.05 0.1 0.15 0.2 0.25
0 0.2 0.4 0.6 0.8 1 1.2
f ∈[1,400] Hz
252.13e+06
r(m)
l (m) 0
0.05 0.1 0.15 0.2 0.25
0 0.2 0.4 0.6 0.8 1 1.2
f ∈[1,700] Hz
r(m)
l (m) 0
0.05 0.1 0.15 0.2 0.25
0 0.2 0.4 0.6 0.8 1 1.2
f ∈[1,700] Hz
48.42.08e+06
r(m)
l (m) 0
0.05 0.1 0.15 0.2 0.25
0 0.2 0.4 0.6 0.8 1 1.2
f ∈[1,1200] Hz
r(m)
l (m) 0
0.05 0.1 0.15 0.2 0.25
0 0.2 0.4 0.6 0.8 1 1.2
f ∈[1,1200] Hz
1422.07e+06
Figure 4: Objective function map and path of the design variables in the successive sub-problems.
—•
Path and final state of the design variables within a given sub-problem;×true solution. The frequency ranges of the sub-problems are indicated within the figures.3.2.4 Bayesian inversion
The MCMC algorithm is here applied to the present problem and the solution for each of the frequency ranges is given in Fig. 5. The posterior densities are consistent with the maps of Fig. 4 and moreover allow for point estimates and credible intervals to be evaluated. These are given in Tab. 2, providing an increasingly accurate and narrow estimation as the frequency range is widened. It is worth noting that, due to the shape of the posterior density the proposal density is relatively wide, enabling the chain to jump between separate branches, e.g. local minima, of the solution. The existence of two solutions is thus apparent in the first sub-problem, as predicted in Eq. (16).
Sub-problem MAP estimate (l, r) l95%(min, max) r95%(min, max) f ∈[1,50] Hz (0.1103, 0.2480) m (0.1025, 1.1631) m (0.0125, 0.2241) m f ∈[1,100] Hz (0.2648, 0.1599) m (0.1035, 0.5981) m (0.1033, 0.2361) m f ∈[1,200] Hz (0.4863, 0.1209) m (0.4268, 0.5352) m (0.1173, 0.1257) m f ∈[1,400] Hz (0.4995, 0.1194) m (0.4955, 0.5035) m (0.1176, 0.1213) m f ∈[1,700] Hz (0.4995, 0.1192) m (0.4973, 0.5017) m (0.1179, 0.1206) m f ∈[1,1200] Hz (0.5000, 0.1196) m (0.4990, 0.5009) m (0.1186, 0.1206) m Table 2: Sub-problem frequency ranges, MAP estimates, and 95% credible intervals.
Figure 5: 2D posterior densities for sub-problems. The black cross denotes the maximum a posteriori (MAP) estimate, the black triangle is the initial condition (based on Matlab’s lsqnonlin estimation), and the red circle is the true value. The dots represent a sample of the MCMC points, PCC is the Pearson correlation coefficient, and f indicates the frequency range of the corresponding sub-problem. The black curves on the bottom and left side of each plot are the 1D marginal densities. For the 1D densities, the light gray zone denotes the narrowest 95% credible interval (bounds are given in Table 2).
3.3 Expansion chamber coupled to a porous material 3.3.1 Model
As a second example, a coupled problem involving an expansion chamber and a porous material is considered, as illustrated in Fig. 6. The porous material is assumed to have a limp frame and is placed adjacent to the air expansion. Such a configuration is expected to ensure a strong coupling between the two components. The inverse estimation problem consists of 5 unknowns, namely the length and radius of the expansion chamber and the porosity, flow resistivity and tortuosity of the porous material. All other parameters are considered to be fixed.
r l
(φ, σ, τ )
Figure 6: Expansion chamber coupled to a porous material. The unknown properties are the length and radius of the expansion chamber and the porosity, flow resistivity and tortuosity of the porous material.
The porous material is modelled by means of the Johnson-Champoux-Allard model [16, 17], providing the equivalent density ˜ρ(ω) and bulk modulus ˜K(ω) as
˜
ρ(ω) = τ ρ0
φ 1 + σφ
iωρ0τ s
1 +iωρ04τ2η σ2Λ2φ2
!
(18) and
K(ω) =˜ γP0
φ
γ − (γ − 1)
"
1 + σφ
iωρ0τ Pr s
1 +iωρ04τ2ηPr σ2Λ02φ2
#−1
−1
, (19)
with the parameters of air defined as ρ0 the mass density, η the viscosity, Pr Prandtl’s number, P0 the atmospheric pressure and γ the heat capacity ratio, and the parameters of the porous material defined as φ the porosity, σ the flow resistivity, τ the tortuosity, Λ the viscous characteristic length and Λ0the thermal characteristic length.
Under the assumption of a limp frame [17, 18], the equivalent density of the porous material can be approximated as
˜
ρlimp(ω) = ρtρ(ω) − ρ˜ 20 ρt+ ˜ρ(ω) − 2ρ0
, (20)
where ρt = ρs+ φρ0 is the total apparent mass density of the porous material, ρsbeing the mass density of the frame. For the purposes of the present model, the equivalent characteristic impedance ˜Z(ω) and wavenumber ˜k(ω) are used, defined as
Z(ω) =˜ q
˜
ρlimp(ω) ˜K(ω), k(ω) = ω˜
sρ˜limp(ω)
K(ω)˜ . (21)
The unknowns of the problem are denoted
x = {l, r, φ, σ, τ } (22)
and the true values of the model parameters are denoted
x0 = {l0, r0, φ0, σ0, τ0}. (23) Table 3 summarises the model parameters and their values, as well as those considered as unknown for the inverse estimation problem.
The full-spectrum estimation and the proposed approach are applied below, using the sound transmission loss of the coupled system as the data to be reproduced by the inverse estimation, with Eq. (14) as the objective function. The numerical experiments are performed with identical settings to the previous example on the single expansion chamber, with the exception of the additive noise, here set to 0.5 dB.
3.3.2 Full-spectrum gradient-based inversion
The inverse estimation is here performed on the full frequency range of the problem. A multi-start scheme with 400 starting points with randomly generated values is used. Figure 7 shows the ensemble of solutions and the histogram of objective function residual values. The figure shows a high sparsity of the objective function residual values and a wide scattering of the model parameters, i.e. the appearance of multiple local minima. The model parameter values can be obtained as the solution exhibiting the lowest objective function residual value. However, such a full-spectrum multi-start procedure is only guaranteed to provide the true solution for a sufficiently large sample of starting points.
Property Symbol Target value Unknown Range Geometry
Expansion chamber length l0 0.5 m l [0.01, 1] m
Expansion chamber radius r0 0.12 m r [0.005, 0.25] m
Porous material length lp 0.15 m Porous material radius rp 0.09 m
Inlet/outlet radius R 0.04 m
Air
Mass density ρ0 1.2 kg·m−3
Viscosity η 1.81·10−5Pa·s
Prandtl’s number Pr 0.71
Atmospheric pressure P0 101325 Pa
Heat capacity ratio γ 1.4
Porous material
Frame density ρs 10 kg·m−3
Porosity φ0 0.98 φ [0, 1]
Flow resistivity σ0 104N·s·m−4 σ [103, 105] N·s·m−4
Tortuosity τ0 3 τ [1, 4]
Viscous characteristic length Λ 180·10−6m Thermal characteristic length Λ0 270·10−6m
Table 3: Values used for air properties and porous material properties for the Johnson-Champoux-Allard model within the limp frame limit. The ranges used for the inversion procedures are specified.
0 0.2 0.4 0.6 0.8 1
0 0.05 0.1 0.15 0.2
0 0.05 0.1 0.15 0.2
0.5 0.6 0.7 0.8 0.9 1
0.5 0.6 0.7 0.8 0.9 1
0 2 4 6 8 10
0 2 4 6 8 10
1 1.5 2 2.5 3 3.5 4 0
0.05 0.1 0.15 0.2 0.25
0 2 4 6 8 10 12 14
r (m) φ σ/10000 (Ns/m4) τ
l(m)r(m)φσ/10000(Ns/m4)
P(d)
d
Figure 7: Parameter values obtained for 400 runs of the full-spectrum gradient-based inversion and histogram of objective function residual values. true solution,solution with the lowest objective function residual.
3.3.3 Low frequency asymptotic behaviour
An expression for the low-frequency asymptotic behaviour of the coupled system may be obtained by considering the limp Johnson-Champoux-Allard model at the limit ω 1, which yields
Z(ω) ' ρ˜ tcp, k(ω) '˜ ω cp
, (24)
where cp =r P0 φρt
is the low-frequency limit of the speed of sound in the limp porous medium. This in turn leads to a low-frequency approximation of the sound transmission loss,
STL(x) ' 10 log10
1 + ω2
l
2cψ r2 R2
+ lp
2cpψ ρ0crp2 ρtcpR2
!!2
− lpl
cpcψ ρ0cr2p ρtcpr2
!
, (25) where
ψ(x) = x + x−1. (26)
The expression is comparable to Eq. (15) in that it is monotonous in l and presents a similar dependence in r. The dependence on the porosity of the foam is weak, only appearing through its total mass density ρt. Furthermore, this low-frequency asymptotic expression does not depend on the flow resistivity or tortuosity of the foam.
The roots of the objective function, Eq. (14), are given by
l±(r, φ) = − α(r)β(φ) − γ(r, φ) ±q
α(r)β(φ) − γ(r, φ)2
− α(r)2 β(φ)2− δ2
α(r)2 , (27)
where
α(r) = 1 2cψ r2
R2
, β(φ) = lp
2cp
ψ ρ0cr2p ρtcpR2
!
, γ(r, φ) = lp
2cpcψ ρ0cr2p ρtcpr2
! , δ2 = (α(r0)l0+ β(φ0))2− 2γ(r0, φ0))l0.
The root l+is positive and is shown in Fig. 8 for φ = 10−6and φ = 1. The figure shows that the behaviour of the coupled system in the low frequency range is governed by the parameters of the expansion chamber, while the dependence on the porosity is rather weak. The first step of the incremental estimation will thus provide the relation between l and r, the rest of the parameters being estimated in the subsequent steps.
0 0.05 0.1 0.15 0.2 0.25
0 0.2 0.4 0.6 0.8 1 1.2
r(m)
l (m)
Figure 8: Root of the objective function in the low frequency limit for the coupled expansion-foam problem.
—
φ = 10−6,—
φ = 1.3.3.4 Incremental inverse estimation
The proposed incremental approach is here applied to the case of a coupled expansion chamber and porous material. The inverse estimation is performed with identical incremental frequency ranges as in the previous example. Analogously to the full-spectrum estimation, a multi-start scheme with 400 randomly generated starting points is used here. Figure 9 shows the ensemble of solutions and the corresponding ensemble of values of the objective function, for the final substep of the method, f ∈ [1, 1200] Hz. The figure shows that the proposed incremental frequency approach invariably converges to the true solution. Furthermore, the histogram of objective function residual values is densely distributed, its dispersion being consistent with the noise model.
0.4994 0.4996 0.4998 0.5 0.5002 0.5004 0.5006 0.5008
0.1190.1195 0.12 0.12050.121
0.119 0.1195 0.12 0.1205 0.121
0.930.940.950.960.970.980.99 1
0.93 0.94 0.95 0.96 0.97 0.98 0.99 1
0.9 0.95 1 1.05 1.1
0.9 0.95 1 1.05 1.1
2.852.92.95 3 3.053.13.153.2 0
0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
0.220.230.240.250.260.270.28
r (m) φ σ/10000 (Ns/m4) τ
l(m)r(m)φσ/10000(Ns/m4)
P(d)
d
Figure 9: Parameter values obtained for 400 runs of the incremental-frequency gradient-based inversion and histogram of objective function residual values. true solution,solution with the lowest objective function residual, × mean value over the ensemble.
3.3.5 Bayesian inversion
The Bayesian inversion is here applied to the coupled problem. The solution for the first and last of the frequency ranges are given in Figs. 10 and 11. The marginal densities for the first frequency range exhibit a strong dependence on the length and radius of the expansion chamber, while the densities associated with the material properties of the porous material are close to the assumed uniform prior model. This result is consistent with the low-frequency asymptotic behaviour derived above, which predicts a strong dependence on the radius and length of the expansion, a weak dependence on the porosity of the foam, and no dependence on the flow resistivity and tortuosity. Furthermore, the point estimates and credible intervals evaluated in Fig. 11 are consistent with the solution using the incremental-frequency gradient-based inversion procedure.
Figure 10: 2D marginal densities for the substep f ≤ 50 Hz data. The black cross denotes the MAP estimate, the black triangle is the initial condition (based on Matlab’s lsqnonlin estimation), and the red circle is the true value. The dots represent a sample of the MCMC points and PCC is the Pearson correlation coefficient.
The black curves on the bottom and left side of each plot are the 1D marginal densities, for which the light gray zone denotes the narrowest 95% credible interval (bounds are given in the figure).
4 Conclusion
A method for the inverse estimation of material and geometrical properties of resonant coupled systems including duct components and porous media is proposed. The method consists in performing an inverse estimation procedure by steps over increasingly broader frequency ranges, in the present case starting from a low frequency limit. It has been shown, here and in a recent paper [9], that the proposed method allows to obtain the global solution of the inverse estimation problem in cases where the full spectrum inversion fails.
A crucial aspect is the availability of a frequency band where the system is non-resonant, as a starting point for the method. This ensures that the solution to the starting sub-problem belongs to a subspace containing the global solution. It is important to note that the method relies on the physical nature of resonant systems and therefore its success is in principle independent from the inversion algorithm used. In fact, both a
Figure 11: Results for the substep f ≤ 1200 Hz. Otherwise same caption as in Fig. 10.
deterministic and a Bayesian framework have been implemented and shown to benefit from the robustness of the procedure for the cases under consideration.
Future work will focus on evaluating the performance of the proposed incremental frequency approach to experimental cases.
The gradient-based implementation of the method presented herein is used in another paper in this conference [19] for the purpose of identifying the geometrical properties of muffler components.
Acknowledgements
This paper is based upon work from COST Action DENORMS CA15125, supported by COST (European Cooperation in Science and Technology). The work has been supported by the strategic funding of the University of Eastern Finland and by the Academy of Finland (Finnish Centre of Excellence of Inverse Modelling and Imaging), and the Centre for ECO2 Vehicle Design at KTH.
References
[1] R. Panneton and X. Olny. Acoustical determination of the parameters governing viscous dissipation in porous media. The Journal of the Acoustical Society of America, 119(4):2027–2040, 2006.
[2] L. Jaouen, A. Renault, and M. Deverge. Elastic and damping characterizations of acoustical porous materials: Available experimental methods and applications to a melamine foam. Applied Acoustics, 69(12):1129–1140, 2008.
[3] P. Leclaire, O. Umnova, K.V. Horoshenkov, and L. Maillet. Porosity measurement by comparison of air volumes. Review of Scientific Instruments, 74(3):1366–1370, 2003.
[4] J. Cuenca, C. Van der Kelen, and P. G¨oransson. A general methodology for inverse estimation of the elastic and anelastic properties of anisotropic open-cell porous materials — with application to a melamine foam. Journal of Applied Physics, 115(8):084904, 2014.
[5] J.P. Parra Mart´ınez, O. Dazel, P. G¨oransson, and J. Cuenca. Acoustic analysis of anisotropic poroelastic multilayered systems. Journal of Applied Physics, 119(8):084907, 2016.
[6] M. Niskanen, J.-P. Groby, A. Duclos, O. Dazel, J.C. Le Roux, N. Poulain, T. Huttunen, and T. L¨ahivaara. Deterministic and statistical characterization of rigid frame porous materials from impedance tube measurements. The Journal of the Acoustical Society of America, 142(4):2407–2418, 2017.
[7] W. Tu and R.W. Mayne. Studies of multi-start clustering for global optimization. International Journal for Numerical Methods in Engineering, 53(9):2239–2252, 2002.
[8] T. L¨ahivaara, L. K¨arkk¨ainen, J.M.J. Huttunen, and J.S. Hesthaven. Deep convolutional neural networks for estimating porous material parameters with ultrasound tomography. The Journal of the Acoustical Society of America, 143(2):1148–1158, 2018.
[9] P. G¨oransson, J. Cuenca, and T. L¨ahivaara. Parameter estimation in modelling frequency response of coupled systems using a stepwise approach. ArXiv:1808.04985 [physics.comp–ph], 2018.
[10] K. Svanberg. A class of globally convergent optimization methods based on conservative convex separable approximations. SIAM Journal on Optimization, 12(2):555, 2002.
[11] J.P. Kaipio and E. Somersalo. Statistical and Computational Inverse Problems. Springer-Verlag, 2005.
[12] D. Calvetti and E. Somersalo. Introduction to Bayesian Scientific Computing: Ten Lectures on Subjective Computing (Surveys and Tutorials in the Applied Mathematical Sciences). Springer-Verlag New York, Inc., 2007.
[13] H. Haario, E. Saksman, and J. Tamminen. Adaptive proposal distribution for random walk Metropolis algorithm. Computational Statistics, 14:375–396, 1999.
[14] H. Haario, E. Saksman, and J. Tamminen. An adaptive Metropolis algorithm. Bernoulli, 7:223–242, 2001.
[15] M.L. Munjal. Acoustics of Ducts and Mufflers With Application to Exhaust and Ventilation System Design. 1987.
[16] Y. Champoux and J.-F. Allard. Dynamic tortuosity and bulk modulus in air-saturated porous media.
Journal of Applied Physics, 70(4):1975–1979, 1991.
[17] J.F. Allard and N. Atalla. Propagation of Sound in Porous Media: Modelling Sound Absorbing Materials 2e. Wiley, 2009.
[18] R. Panneton. Comments on the limp frame equivalent fluid model for porous media. The Journal of the Acoustical Society of America, 122(6):EL217–EL222, 2007.
[19] M.H. Alkmim, J. Cuenca, L. De Ryck, and P. G¨oransson. Model-based acoustic characterisation of duct components and extrapolation to inhomogeneous thermal conditions. In ISMA/USD, 2018/9/17-19.
