Thomas Hélie
email: thomas.helie@ircam.fr

Béatrice Laroche
email: beatrice.laroche@inra.fr

SIMULATION OF A DAMPED NONLINEAR BEAM BASED ON MODAL DECOMPOSITION AND VOLTERRA SERIES

Keywords: sound and vibration, partial differential equation, nonlinear systems, Volterra series, convergence

This paper addresses the model order reduction and the simulation of a damped nonlinear pinned beam excited by a distributed force. The model is based on : (H1) the assumption of Euler-Bernoulli kinematics (any cross-section before deformation remains straight after deformation); (H2) Von Karman's assumptions which couples the axial and the bending movements, introducing a nonlinearity in the model; (H3) some viscous and structural damping phenomena. The problem is first described and its linearized version is analyzed. This is used to build a reduced order model based on a standard modal decomposition. Then, the nonlinear system is examined in the framework of the regular perturbation theory. It is solved based on a Volterra series approach : the vibration is decomposed into a sum of nonlinear homogeneous contributions with respect to the excitation. The convergence of the series and the truncation error are characterized. Finally, numerical experiments are presented and discussed.

Model

Consider the Euler-Bernoulli model [START_REF] Géradin | Théorie des vibrations[END_REF][START_REF] Graff | Wave motion in elastic solids[END_REF] of a damped nonlinear pinned beam, initially at rest. The dimensionless model that governs the deflection waves w(z, t) are given by, for all z ∈ Ω =]0, 1[ and t ∈ T = R + ,

∂ 2 t w + 2 (a + b∂ 4 z ) ∂ t w + ∂ 4 z w -η 1 0 ∂ z w 2 dz ∂ 2 z w = f, (1)
with pinned-type boundary conditions at z ∈ {0, 1}, namely, w(z, t) = 0 (fixed extremities) and ∂ 2 z w(z, t) = 0 (no momentum), with zero initial conditions at t = 0, namely, w(z, 0) = 0 and ∂ t w(z, 0) = 0, and where f (z, t) is a distributed force. Coefficients a > 0 and b > 0 are fluid and structural damping parameters and η > 0 is the nonlinear coupling coefficient between the bending momentum and the displacement under the von Kármán assumption [START_REF] Nayfeh | Nonlinear Oscillations[END_REF].

2. Linearized problem (η = 0)

Modal decomposition and order reduction

For η = 0, the problem is linear and well-posed in an appropriate functional setting that is not detailed here (see e.g. [START_REF] Jacob | Analyticity and Riesz basis property of semigroups associated to damped vibrations[END_REF]). The solution w admits a modal decomposition w(z, t) = m≥1 w m (t)e m (z) on the basis of eigenfunctions (orthonormal in L 2 (Ω))

{e m } m≥1 with e m (z) = √ 2 sin(k m z) for z ∈ Ω and k m = mπ.
If only a finite number of modes E = [e 1 , . . . , e M ] T are excited by the distributed force, namely if

f (z, t) = E(z) T F (t) with F = [f 1 , . . . , f M ] T , then the exact solution w(z, t) = E(z) T W (t) with W = [w 1 , . . . , w M ] T is governed by ẅm + 2 (a + bk 4 m ) ẇm + k 4 m w m = f m , for 1 ≤ m ≤ M.
This problem can be restated as a state-space representation with input u = F , state x = [W T , Ẇ T ] T , where, denoting L = πdiag(1, . . . , M ),

ẋ(t) = A x(t) + B u(t) with A = 0 M ×M I M -L 4 -2(a I M + b L 4 ) , B = 0 M ×M I M , (2)
x(0) = 0 2M ×1 . (3)

Solution and properties

The solution x lin of this linear problem is the convolution of input u by the semi-group S, represented by the block diagram u → S → x lin and given by

x lin (t) = ∞ 0 S(t -τ )Bu(τ ) dτ with S(t) = exp(A t)1 t>0 . (4)
A consequence of the damping is that for all bounded input u ∈ U = L ∞ (T, U) living in U = R M , the trajectory x lin ∈ X = L ∞ (T, X) living in X = R 2M is also bounded. More precisely, S(t) L(X,X) is exponentially decreasing so that

x lin X ≤ γ B L(X,U) 1 u U where γ ≥ γ := T S(t) L(X,X) dt > 0. (5)
For musical applications, the damping coefficients are supposed to be such that: (i) the first mode (m = 1) has a damped oscillating dynamics, (ii) higher modes are more damped. By analyzing the complex-valued roots1 of the characteristic polynomials P m (λ) = λ 2 +2(a+bk 4 m )λ+k 4 m of (2.1), these conditions are found to be met if a π 2 + bπ 2 < 1 (i), and if 2b(abπ 4 ) ≤ 1 (ii). In this case, the real part of these eigenvalues are lower than -α < 0 and γ is found to be bounded by γ, with α = a + bπ 4 and γ = 1/α.

3. problem (η > 0)

Nonlinear state-space representation

Based on the same modal decomposition on E, a straightforward computation shows that (1) admits the following exact2 reduced order model

ẋ(t) = A x(t) + B u(t) + A 3 x(t), x(t), x(t) and x(0) = 0 2M ×1 , (7)
with

A 3 (a, b, c) = η a T L 2 0 M ×M 0 M ×M 0 M ×M b 0 M ×M 0 M ×M -L 2 0 M ×M c (multi-linear operator). ( 8
)
This multi-linear operator is bounded and an estimate of the bound (derivation not detailed here) yields

η 2 √ 10 ≥ A 3 = sup x 1 X =•••= x 3 X =1 A 3 (x 1 , x 2 , x 3 X . (9)

Solution based on a Volterra series expansion

The trajectory x of (7-8) can be represented by a Volterra series expansion [START_REF] Volterra | Sopra le funzioni che dipendono de altre funzioni[END_REF][START_REF] Volterra | Theory of Functionnals and of Integral and Integro-Differential Equations[END_REF] (see also [START_REF] Hélie | Sound synthesis of a nonlinear string using Volterra series[END_REF][START_REF] Roze | Introducing a Green-Volterra series formalism to solve weakly nonlinear boundary problems: application to the Kirchhoff's string[END_REF] for applications on nonlinear damped vibrating strings): it can be decomposed as a sum of contributions with homogeneous of order m w.r.t. to input u given by

x(t) = +∞ m=1 x m (t), (10)
where

x 1 (t) = x lin (t) = t 0 S(t -τ ) χ 1 (τ ) dτ, with χ 1 (τ ) = B u(τ ) (11)
and, for m ≥ 2,

x m (t) = t 0 S(t -τ ) χ m (τ ) dτ with χ m (τ ) = m k=2 p∈M 3 k A 3 x p 1 (τ ), . . . , x p k (τ ) , ( 12
)
where the index set M 3 m := p ∈ (N * ) 3 p 1 + p 2 + p 3 = m selects occurrences x p i such that their combination through the multi-linear operator A 3 is of homogeneous order m. Deriving the first contributions yields χ 2p = 0 for p ≥ 1 (only odd nonlinear orders are activated) and

χ 3 = A 3 (x 1 , x 1 , x 1 ), χ 5 = A 3 (x 1 , x 1 , x 3 ) + A 3 (x 1 , x 3 , x 1 ) + A 3 (x 3 , x 1 , x 1 ) = A 3 (x 1 , x 1 , x 3 ) + 2A 3 (x 1 , x 3 , x 1 ) (using the sym- metry A 3 (a, b, c) = A 3 (b, a, c)), χ 7 = A 3 (x 3 , x 3 , x 1 )+2A 3 (x 1 , x 3 , x 3 )+A 3 (x 1 , x 1 , x 5 )+2A 3 (x 1 , x 5 , x 1 ), etc.
The interest for simulation is that solution (10-12) only involves linear filtering (S) and static functions (A 3 ), to be combined through sums.

Convergence and error bound

According to [START_REF] Hélie | Computation of convergence bounds for Volterra series of linear analytic single input systems[END_REF][START_REF] Hélie | Computable convergence bounds of series expansions for infinite dimensional linearanalytic systems and application[END_REF], the convergence of the series (7-8) can be characterized by a radius ρ, a gain bound function Φ and an error bound function R M Φ as follows:

if x 1 X < ρ, then x X < Φ x 1 X , (convergence in norm) (13) and x - M m=1 x m X < R M Φ x 1 X . (truncation error bound) (14)
Based on [START_REF] Hélie | Computation of convergence bounds for Volterra series of linear analytic single input systems[END_REF][START_REF] Hélie | Computable convergence bounds of series expansions for infinite dimensional linearanalytic systems and application[END_REF] and using estimates ( 6) and ( 9), the radius ρ and the functions Φ and R M Φ are found to be (see also figure 1)

ρ = 2 4 √ 10 3 α η , ( 15
) ∀z ∈ [0, ρ[, Φ(z) = 3ρ cos π + arccos(z/ρ) 3 , (16)
R M Φ(z) = Φ(z) - M m=1 φ m z m , (Taylor expansion remainder of Φ) (17)
where the Taylor coefficients are φ 2p = 0 for all p ≥ 1 and

φ 1 = 1, φ 3 = 2 2 3 3 ρ 2 , φ 5 = 2 4 3 5 ρ 4 , φ 7 = 2 8 3 8 ρ 6 , φ 9 = 2 8 ×5×11
3 12 ρ 8 , etc.

10 onlinear system is simulated using a standard ODE solver (ode15s, These trajectories are compared to those obtained from truncated

-1 10
0 ϵ =10 -1 ϵ =10 -2 ϵ =10 -3 z/ρ 10 - 2
Φ/ρ M = 1 M = 3 M = 5 M = 7 M = 9 z/ρ R M Φ/ρ b ⃝ a ⃝ Function D ϵ Φ/ρ where D ϵ Φ(z) = Φ z(1 + ϵ) -Φ(z)

Numerical experiments and

We consider a modal excitation f (z, t) = M m=1 e m (z)u m (t) with a uniform impulsion u m = u 0 1 [0,T ] (t) with T > 0. Two configurations are examined.

Configuration 1

We consider a single mode (M = 1). In this case, the nonlinear problem corresponds to a Duffing oscillator. The damping is chosen to be of fluid type only (b = 0) and tuned to be close to the critical regime (a = 0.999π 2 ). Parameter η is chosen such that ρ = 1. The excitation duration is T = 3. Four amplitudes u 0 are tested such that x 1 X /ρ ∈ {0.8; 1; 1.2; 2}. Signals w(z = 0.5, t) = x 1 (z = 0.5, t) are displayed in figure 2 for several truncation orders (M = 1, 3, 5, 7) of [START_REF] Volterra | Sopra le funzioni che dipendono de altre funzioni[END_REF]. They are simulated by combining (exact) convolutions by S with (exact) operator A 3 , according to [START_REF] Volterra | Theory of Functionnals and of Integral and Integro-Differential Equations[END_REF](12) For x X = 0.8, a good approximation is obtained as soon as M ≥ 3. As guaranteed by (13), the convergence is numerically observed. For x X = 1, the approximations are significantly more accurate when increasing the order M (see also the zoom in figure 3). This is no longer true for x X = 1.2, for which the convergence seems to be lost. For x X = 2, the divergence is so fast (see the time range 1 ≤ t ≤ 3) that the best approximation is the linear approximation.

As a conclusion, the guaranteed bound ρ is close to the exact convergence radius.

x = 0.8 x X = 1 0 1 2 3 4 5

Configuration 2

The second configuration corresponds to an oscillating beam. The damping parameters a = 0.02 and b = 5 × 10 -5 are chosen to produce "wooden beam" sounds3 . As for configuration 1, η is chosen such that ρ = 1. The excitation is chosen such that M = 3 with T = 200. Numerical results are displayed in figure 4 for three amplitudes u 0 such that x 1 X /ρ ∈ {1; 3; 5}.

x X complete signal zoom For x 1 X = 1, the nonlinear dynamics is very close to the linear approximation and the convergence is fast. For x 1 X = 3, increasing the truncation order from M = 1 to M = 7 numerically improves the results (see the zoomed part). Although the convergence is no longer guaranteed, it can still be valid as ρ is a (guaranteed) lower bound of the maximal convergence radius. For x 1 X = 5, the convergence seems to be lost and nonlinear contributions x m≥3 build unrealistic high amplitudes signals in the first part with high magnitudes (see the zoomed part between t = 99 and t = 102).

For large input signals (outside the convergence domain), high order contributions x m have increasing amplitudes with m.

Outside the convergence domain, the high magnitudes observed on signals x m for orders m ≥ 3 correspond to the so-called phenomenon of "secular modes" [START_REF] Nayfeh | Nonlinear Oscillations[END_REF]. This phenomenon is well-known for the Duffing oscillator. It is due to the nonlinear effect type: in this system, the nonlinearity is mainly responsible for a frequency modulation of the input signal for large magnitudes. The regular perturbation method represents and approximates this frequency modulation with polynomial combinations of oscillating signals with fixed frequencies.

In the case of a conservative problem, this kind of approximation produces signal envelopes that increase as as t p-1 for orders m = 2p + 1. As a consequence, there is no possible convergence over T = R and ρ is zero. In the case of the damped beam, the damping makes these envelopes behaves as t p-1 exp e(λ n )α n t for each mode n, so that signals x m are all bounded and the convergence radius is nonzero. In other words, the convergence radius ρ appear to be a guaranteed bound for which high order contributions x m are not secular over T.

Conclusion and perspectives

The Volterra series approach is a regular perturbation method that provides external representation of nonlinear input-output systems, as transfer functions or impulse convolution do for linear time invariant problems. This has been successfully used to represents the vibrations of a damped nonlinear beam. The convergence analysis provides a computable domain, available for bounded excitation signals with any waveshape, for which (due to the damping) no divergence is activated by secular modes. In this domain, the series expansion provides accurate approximations of the solution: this is well-adapted to account for distortions and timbres modification of sounds and vibrations produced by any input waveshape. Outside this domain, secular modes appear and express that frequency modulations become a significant phenomena. A perspective is to propose extensions of this work to efficiently represent such modulations, based on other perturbation methods and approaches that are adapted to any input signal waveform.

ation 1 .

 1 for relative errors ϵ of and 10%. b ⃝ functions R M Φ/ρ (M = 1, . . . , 9) and Φ/ρ.erical resultssets of parameters are studied. We consider a single mode (n = N = 1), a fluid damping 0) tuned to be close to the critical regime (a = 0.999π 2 ). Parameter en such that ρ = 1, which corresponds to η = 13.857. The excitation is θ = 3. Four amplitudes are tested, chosen such that ∥x 1 ∥ X /ρ ∈ .2; 2}.

Figure 1 :

 1 Figure 1: Normalized gain bound function Φ/ρ and normalized remainder functions R M

6 -Figure 2 :

 62 Figure 2: Simulations for configuration 1 (the reference is computed by the standard ODE solver ode15s, Matlab).

Figure 3 :

 3 Figure 3: Zoom for the case x X = 1 in figure 2.

Figure 4 :

 4 Figure 4: Simulations for configuration 2 (oscillating wooden beam)

A similar analysis can be found in[START_REF] Hélie | Damping models for the sound synthesis of bar-like instruments[END_REF] (detailed analysis of poles in the Laplace domain), for the case of a beam excited by a force at one extremity and free boundary conditions.

In[START_REF] Géradin | Théorie des vibrations[END_REF], the nonlinear operator is a Laplace operator multiplied by a scalar (depending on an integral of (∂ z w) 2 ), that preserves the co-linearity of functions e m .

To listen to a result with a first mode at frequency f 0 on a sound card with sampling frequency f s , the simulated trajectories is sampled at T = 2πf 0 /(f s ms 1 ) where (s 1 , s 1 ) are the poles associated with the first mode (T ≈ 0.00292 for f 0 = 220Hz, f s =

48000Hz).

ICSV26, Montreal, 7-11 July 2019