John M Campbell
email: jmaxwellcampbell@gmail.com

π involving harmonic numbers

Keywords: Ramanujan-like series, Harmonic number, Pi formula, Complete elliptic integral Mathematics Subject Classification (2000) Primary 33C75, 33C20, Secondary 65B10

We introduce new classes of Ramanujan-like series for 1 π , by devising methods for evaluating harmonic sums involving squared central binomial coefficients, such as the Ramanujan-type series

introduced in this article. While the main technique used in this article is based on the evaluation of a parameter derivative of a beta-type integral, we also show how new integration results involving complete elliptic integrals may be used to evaluate Ramanujan-like series for 1 π containing harmonic numbers.

Introduction

The evaluation of summations containing binomial coefficients and harmonic numbers is an interesting topic, and a large amount of mathematical literature is devoted to this area. In this article, we consider the general problem of evaluating infinite series involving harmonic numbers and squared central binomial coefficients.

Many harmonic sums involving binomial coefficients may be evaluated by applying differential operators to known hypergeometric identities. In particular, the beautiful Ramanujan-like formula

∞ n=1 2n n 2 H n 16 n (2n -1) 2 =
12 -16 ln [START_REF] Chen | Interesting series associated with central binomial coefficients, Catalan numbers and harmonic numbers[END_REF] π which was noted in 2016 by Hongwei Chen may be proven using this method [START_REF] Chen | Interesting series associated with central binomial coefficients, Catalan numbers and harmonic numbers[END_REF]. The author of [START_REF] Chen | Interesting series associated with central binomial coefficients, Catalan numbers and harmonic numbers[END_REF] strongly encourages further exploration of properties of series involving powers of binomial coefficients and harmonic numbers.

Although certain types harmonic summations may be evaluated by applying differential operators to specific hypergeometric identities, it is not in general obvious as to how to evaluate a given a harmonic sum involving powers of binomial coefficients. Inspired in part by [START_REF] Chen | Interesting series associated with central binomial coefficients, Catalan numbers and harmonic numbers[END_REF], we introduce in our present article a very general method for evaluating series with a summand involving a factor of the form 2n n 2 H n . Using this method, we prove a variety of new Ramanujan-like formulas for 1 π . The main method introduced in this article for evaluating harmonic sums containing expressions of the form 2n n 2 is based on a parameter derivative of a beta-type integral. This technique is very useful and general, because this method may be used to evaluate series of the form

∞ n=1 g(n) 2n n 2 H n
for a given function g on N, whereas techniques given by differentiating specific classes of hypergeometric identities will only produce specific results. We also show how a generating function recently considered by Boyadzhiev in [START_REF] Boyadzhiev | Series with central binomial coefficients, Catalan numbers, and harmonic numbers[END_REF] may be used to evaluate certain harmonic sums with squared central binomial coefficients, and we also describe a general method for evaluating harmonic sums of the form

∞ n=1 g(n) 2n n 2 (H 2 n + H (2) n ).
The elegant Ramanujan-type formula

∞ n=0 (-1) n 1 2 3 n (1) 3 n (-2 + 3(4n + 1)H n ) = - 12 ln(2) π
is proven in [START_REF] Guillera | More hypergeometric identities related to Ramanujan-type series[END_REF], and several other Ramanujan-like formulas related to harmonic numbers are proven in [START_REF] Guillera | More hypergeometric identities related to Ramanujan-type series[END_REF]. Apart from [START_REF] Chen | Interesting series associated with central binomial coefficients, Catalan numbers and harmonic numbers[END_REF] and [START_REF] Guillera | More hypergeometric identities related to Ramanujan-type series[END_REF], there does not seem to be very much mathematical literature concerning harmonic sums for 1 π . It is thus natural to use a systematic approach towards the construction of new series of this form.

The formulas for 1 π introduced in this article can also be proven using some new formulas for definite integrals involving complete elliptic integrals, as we discuss in Section 4. We are interested in the general problem of evaluating integrals involving complete elliptic integrals in part because there are many applications related to elliptic functions and elliptic integrals in various fields connected to physics and engineering. As discussed in [START_REF] Kaplan | Multiple elliptic integrals[END_REF], many physical problems involve multiple integrals involving elliptic integrals such that once one of the integrations is completed, the resultant integrand contains an elliptic integral.

Harmonic sums involving squared central binomial coefficients

In this section, we describe some different techniques for evaluating Ramanujanlike series involving harmonic numbers and expressions of the form 2n n 2 , and in Section 3 we apply our main technique, which is described in Section 2.1, by introducing new harmonic series for 1 π .

First strategy

Our first strategy for constructing new Ramanujan-type series for 1 π is based on the evaluation of integrals of the form 1 0

x n ln(1-x 2 ) √ 1-x 2 dx.
Integrals of this form may be evaluated in a natural way in terms of harmonic numbers, as indicated in Lemma 1 below. We remark that a somewhat similar strategy based on the evaluation of integrals involving an expression of the form x ln 2m-1 x 1-x was recently used to evaluate certain classes of series involving harmonic numbers in [START_REF] Sofo | Integrals of logarithmic and hypergeometric functions[END_REF].

Lemma 1 Let f : N 0 → C be such that the series ∞ n=0 (-1) n 1 2 n f (n) x 2n ln 1 -x 2 √ 1 -x 2
is integrable on [0, 1], and write

g(n) = f (n) 16 n (2n-1) . Then ∞ n=0 g(n) 2n n 2 H n is equal to 2 π times the sum of 1 0 ∞ n=0 (-1) n 1 2 n f (n) x 2n ln 1 -x 2 √ 1 -x 2 dx and π ln(2) ∞ n=0 ( 2n n ) 2 f (n) 16 n (1-2n) .
Proof This follows immediately from the identity

1 0 x 2n ln 1 -x 2 √ 1 -x 2 dx = - π2 -2n nΓ (2n) (H n + 2 ln(2)) Γ (n + 1) 2 , which holds for n such that (n) > -3 2 .
The above lemma is especially useful because it is not in general otherwise obvious as to how to evaluate a series of the form

∞ n=0 g(n) 2n n 2 H n .
Since one of the key tools used in this paper is based on the evaluation of integrals of the form 1 0

x n ln(1-x 2 ) √ 1-x 2 dx, it is worthwhile to note that this integral is a derivative of the beta function with respect to a parameter. In particular, letting

B(x, y) = 1 0 t x-1 (1 -t) y-1 dt, we have that 1 2 ∂ ∂y B n + 1 2 , y + 1 = 1 0 u n 1 -u 2 y log 1 -u 2 du, so that 1 0 u n log 1 -u 2 √ 1 -u 2 du = √ πΓ 1+n 2 ψ (0) 1 2 2Γ 1 2 + 1+n 2 - √ πΓ 1+n 2 ψ (0) 1 2 + 1+n 2 2Γ 1 2 + 1+n 2
for n ∈ N 0 .

Boyadzhiev's generating function

Boyadzhiev recently evaluated the generating function for the sequence

2n n H n : n ∈ N 0
and similar sequences in [START_REF] Boyadzhiev | Series with central binomial coefficients, Catalan numbers, and harmonic numbers[END_REF], proving that

∞ n=0 2n n H n x n = 2 √ 1 -4x ln 1 + √ 1 -4x 2 √ 1 -4x
for |x| < 1 4 . Boyadzhiev's generating function may be used to construct certain types of Ramanujan-like series for 1 π involving squared central binomial coefficients and harmonic numbers. But there are certain series of this form which cannot feasibly be proven directly using Boyadzhiev's generating function or similar generating functions, and which can be proven easily using Lemma 1, as we later discuss.

To construct a Ramanujan-type series for 1 π using Boyadzhiev's generating function one may use the following integral identity:

2n n 4 n = 2 π ∞ 0 dx (x 2 + 1) n+1 .
To illustrate this idea, as well as the technique described in Section 2.1, we offer two corresponding proofs of the following new Ramanujan-like formula.

Theorem 1 ∞ n=1 ( 2n n ) 2 Hn 16 n (2n-1) = 8 ln(2)-4 π .
Proof #1: This follows immediately from Lemma 1 by letting f (n) = 1.

Proof #2: By Boyadzhiev's generating function, we have that

∞ n=0 2n n H n x 2n-2 4 n = 2 ln √ 1-x 2 +1 2 √ 1-x 2 x 2 √ 1 -x 2 .
Integrating both sides of this equality, we have that

∞ n=0 ( 2n n )Hny 2n-1 4 n (2n-1)
is equal to:

2 -1 -y 2 -1 -y 2 ln 1 √ 1-y 2 + 1 + 1 -y 2 ln(2) + 1 y .
Therefore,

∞ n=0 ( 2n n )Hn2 1 x 2 +1 n+1 4 n π(2n-1)
is equal to

4 √ x 2 + 1 - √ x 2 ln 1 x 2 x 2 +1 + 1 + √ x 2 (ln(2) -1) π (x 2 + 1) 3/2
.

Mathematica is able to evaluate the integral

∞ 0 4 √ x 2 + 1 - √ x 2 ln 1 x 2 x 2 +1 + 1 + √ x 2 (ln(2) -1) π (x 2 + 1)
3/2 dx as 8 ln(2)-4 π , thus proving the above evaluation of

∞ n=1 ( 2n n ) 2 Hn
16 n (2n-1) .

In Section 4, we offer another elegant proof of Theorem 1, through the use of integrals involving complete elliptic integrals.

Our main strategy introduced in Section 2.1 for evaluating Ramanujantype series for 1 π involving harmonic numbers is much more general and powerful compared to the alternative technique outlined in Section 2.2. For example, it seems that it would not be feasible to use Boyadzhiev's generating function to prove the following result due to Hongwei Chen, although Lemma 1 may be used to prove the following theorem.

Theorem 2 ∞ n=1 ( 2n n ) 2 Hn 16 n (2n-1) 2 = 12-16 ln(2) π [2].
Proof #1 : This follows immediately from Lemma 1 by letting f (n) = 1 2n-1 . By Boyadzhiev's generating function, we have that

∞ n=0 ( 2n n )Hny 2n-2 4 n (2n-1) is equal to 2 -1 -y 2 -1 -y 2 ln 1 √ 1-y 2 + 1 + 1 -y 2 ln(2) + 1 y 2 .
Integrating the above expression yields a very complicated expression involving the polylogarithm function. Mathematica is not able to evaluate the corresponding integral needed to evaluate the series

∞ n=1 ( 2n n ) 2 Hn
16 n (2n-1) 2 using the identity

( 2n n ) 4 n = 2 π ∞ 0 dx (x 2 +1) n+1
. This illustrates how the strategy introduced in Section 2.1 may be used to construct new Ramanujan-like series for 1 π which cannot be evaluated directly following the technique given in Section 2.2. We later offer an alternative proof of Theorem 2 using definite integrals involving complete elliptic integrals. We also offer an alternative proof of the following theorem using integrals containing complete elliptic integrals.

Theorem 3 ∞ n=1 ( 2n n ) 2 Hn
16 n (2n-3) = 120 ln(2)-68 27π .

Proof #1: This follows immediately from Lemma 1 by letting f (n) be equal to 2n-1 2n-3 .

Series involving an expression of the form

H 2 n + H (2) n
In Section 5, we show how the evaluation of integrals of the form

1 0 x n ln 2 1 -x 2 √ 1 -x 2 dx
may be used to prove new formulas for Ramanujan-like series, such as the elegant formula

∞ n=1 2n n 2 H 2 n + H (2) n 16 n (2n -1) = 4π 3 + 16 2 ln(2) -2 ln 2 (2) -1 π
introduced in this article.

3 New Ramanujan-type series for 1 π involving harmonic numbers

We begin by presenting new Ramaujan-like series for 1 π with a summand of the form ( 2n n )

2 Hn 16 n pn , letting p n denote a polynomial with integer coefficients. Lemma 1 may be used to evaluate summations of this form in the case whereby p n is a polynomial with two different factors, but this case can be reduced to simpler cases by decomposing the appropriate fraction into simpler fractions. So, we omit consideration of summands of the form

( 2n n ) 2 Hn
16 n pn in the case whereby p n has two different factors.

The Ramanujan-type formulas given in this section are new in the sense that state-of-the-art computer algebra systems are not able to directly evaluate the series given in this section. For example, Mathematica 11 is only able to evaluate the series

∞ n=1 ( 2n n ) 2 Hn
16 n (2n-1) as

2 F 1 (0,0,1,0) - 1 2 , 1 2 , 1, 1 .
The Ramanujan-like formulas given in this section are also new in the sense that these results have apparently not appeared in any mathematical literature concerning summations involving harmonic numbers and binomial coefficients, with the exception of the formula which is given in Theorem 2. The following new Ramanujan-like formulas are easily verified using Lemma 1.

∞
∞ n=1 2n n 2 H n 16 n (2n -1) = 8 ln(2) -4 π ∞ n=1 2n n 2 H n 16 n (2n -3) = 120 ln(2) -68 27π ∞ n=1 2n n 2 H n 16 n (2n -5) = -6508 + 10680 log(2) 3375π ∞ n=1 2n n 2 H n 16 n (2n -1) 2 = 12 -16 ln(2) π ∞ n=1 2n n 2 H n 16 n (2n -3) 2 = 164 -176 ln(2) 27π ∞ n=1 2n n 2 H n 16 n (2n -5) 2 = 70724 -67760 log(2) 16875π .
We remark that analogues of Lemma 1 may be used to evaluate series involving expressions of the form H 2n . For example, an analogue of Lemma 1 may be used to prove that

∞ n=1 2 -6n 2n n 4n 2n H 2n 2n -1
is equal to the following expression involving the dilogarithm:

1 4π 2Li 2 17 -12 √ 2 -4 √ 2+20 √ 2 ln(2)+ln 19601 -13860 √ 2 +cos -1 (3) 2 + cosh -1 (17) 2 + 2 ln 3 -2 √ 2 sinh -1 (1) -π 12 .
We remark that the Ramanujan-like series introduced in this section may be rewritten as double series using identities such as the following series identity:

∞ n=1 1 n(m + n) = ψ (0) (m + 1) + γ m .
For example, it is easily seen that the formula

Definite integrals involving complete elliptic integrals

It is convenient for our purposes to define the complete elliptic integral of the first kind so that K(k) = 1 2 π 2 F 1 1 2 , 1 2 ; 1; k . Similarly we define the complete elliptic integral of the second kind so that E(k) = 1 2 π 2 F 1 1 2 , -1 2 ; 1; k . The series for 1 π given in the previous section are closely related to properties concerning complete elliptic integrals. To illustrate this idea, we begin by showing how the known integral formula Proof #3 of Theorem 1: Using integration by parts with respect to the integral

1 0 ln(1 -x)K(x) dx = 8 ln(2) -8, by letting u = ln(1 -x) and v = K(x), with v(x) = 2xK(x) -2K(x) + 2E(x) -2,
we have that 8 log(2) -8 is equal to

1 0 -2 + 2E(x) -2K(x) + 2xK(x) 1 -x dx.
Therefore,

8 log(2) -8 = -2 1 0 1 -E(x) 1 -x dx + 1 0 -2K(x) + 2xK(x) 1 -x dx.
So, we have that

2 -4 ln(2) = 1 0 1 -E(x) 1 -x dx (1)
So, since

1 -E(x) = 1 2 π ∞ n=0 1 2 n -1 2 n (1 -x n ) (1) n n! ,
we have that

1 -E(x) 1 -x = 1 2 π ∞ n=0 1 2 n -1 2 n 1-x n 1-x (1) n n! .
By interating both sides of this equality, we find that

2 -4 log(2) = 1 2 π ∞ n=0 1 2 n -1 2 n H n (1) n n! ,
as desired.

The equality given in (1) turns out to be unexpectedly useful, since (1) may be used to construct an alternative proof of Theorem 2, as shown below. The formula in (1) can also be used to construct another proof of Theorem 3. We remark that state-of-the-art computer algebra systems are not able to evaluate the definite integral given in [START_REF] Boyadzhiev | Series with central binomial coefficients, Catalan numbers, and harmonic numbers[END_REF], and that this integral does not seem to appear in any mathematical literature concerning definite integrals.

Proof #2 of Theorem 2: Rewriting the equality

1 0 K(x) dx = 2 as 1 0 - K(x) 2(-1 + x) + xK(x) 2(-1 + x) dx = 1, from (1)
, we have that

1 0 -1 + E(x) -1 + x - K(x) 2(-1 + x) + xK(x) 2(-1 + x) dx = 3 -4 ln(2).
Rewrite this equality as follows:

1 0 π 2(x-1)K(x) π + 4E(x) π -4 π(x -1) dx = 12 -16 ln(2) π .
This equality may be rewritten as follows:

1 0 -4 + π 2 F 1 -1 2 , -1 2 ; 1; x π(-1 + x) dx = 12 -16 ln(2) π .
But since

∞ n=0 16 -n (1 -x n ) 2n n 2 (-1 + 2n) 2 (1 -x) = π 2 F 1 -1 2 , -1 2 ; 1; x -4 π(x -1)
, by integrating both sides of the above equality appropriately we thus have that ∞ n=1

( 2n n ) 2 Hn
16 n (2n-1) 2 = 12-16 ln(2) π as desired.

Proof #2 of Theorem 3: Begin by evaluating

∞ n=0 ( 2n n ) 2 ( 1-x n 1-x ) 16 n (2n-3) as 10 -3π 3 F 2 -3 2 , 1 2 , 1 2 ; -1 2 , 1; x 9π(-1 + x) ,
which is equal to

10 -3π 4(1-x)K(x) 3π + 2(4x+1)E(x) 3π 9π(x -1) .
So, the infinite series given in Theorem 3 is equal to the following integral:

1 0 10 -3π 2(1+4x)E(x) 3π + 4(1-x)K(x) 3π 9π(-1 + x) dx.
Therefore, the series given in Theorem 3 is equal to:

8 9π + 1 0 10 9π(-1 + x) - 2E(x) 9π(-1 + x) - 8xE(x) 9π(-1 + x) dx,
which is equal to:

8 9π - 2 9π 1 0 1 -E(x) 1 -x dx - 8 9π 1 0 1 -xE(x) 1 -x dx.
We can find a closed-form evaluation of the above expression using (1), since

1-xE(x) 1-x = E(x) + 1-E(x) 1-x .
Similar approaches may be used to construct alternative proofs for the other formulas for 1 π given in Section 3. We remark that new formulas for integrals containing complete elliptic integrals, such as the formula

1 0 ln 2 (1 -x)dE(x) = 8 - 2π 2 3 + 16(ln(2) -1) ln(2),
may be proven by rewriting the harmonic numbers in the series given in Section 3 and Section 5 using integrals such as 1 0 x m ln(1 -x) dx. Many integrals of this form can also be proven using integration by parts and by applying known results concerning complete elliptic integrals, so we omit a full investigation of these integrals.

Summations with sqaured harmonic numbers

Given the variety of integration results and Ramanujan-type series related to Lemma 1, it is natural to consider generalizations of the strategy introduced in Section 2.1. We begin by considering the integral

1 0 x n ln 2 1 -x 2 √ 1 -x 2 dx,
which is equal to the following expression for (n) > -5.

1 4nΓ n 2 +1 Γ n 2 √ πΓ n+1 2 4Γ n 2 +1 ψ (0) n 2 +1 2 -ψ (1) n 2 +1 +nΓ n 2 4(γ+ ln(4))ψ (0) n 2 + 1 + π 2 + 2γ 2 + 8 ln 2 (2) + 8γ ln(2) . Theorem 4 ∞ n=1 ( 2n n ) 2 (H 2 n +H (2)
n )

16 n (2n-1)

= 4π 3 + 16(2 ln(2)-2 ln 2 (2)-1) π . Proof Since ∞ n=0 (-1) n x 2n 1 2 n ln 2 1 -x 2 √ 1 -x 2 = ln 2 1 -x 2 ,
we have that

∞ n=0 (-1) n 1 2 n 1 0 x 2n ln 2 1 -x 2 √ 1 -x 2 dx = 8 - π 2 3 + 4 ln 2 (2) -8 ln(2).
Expanding the above summand using the above evaluation of the expression 1 0

x n ln 2 (1-x 2 ) √ 1-x 2
dx together with Theorem 1 may be used to evaluate

∞ n=1 ( 2n n ) 2 ((Hn) 2 +H (2) n ) 16 n (2n-1)
.

The following propositions may be proven by analogy with our proof of Theorem 4.

Proposition 1 ∞ n=1 ( 2n n ) 2 (H 2 n +H (2) n ) 16 n (2n-1) 2 = 64+64 ln 2 (2)-96 ln(2) π -8π 3 . Proposition 2 ∞ n=1 ( 2n n ) 2 (H 2 n +H (2)
n )

16 n (2n-3)

= 20π 27 + 16(-53-90 ln2 (2)+102 ln(2)) 81π

.

Proposition 3 ∞ n=1 2 -4n ( 2n n ) 2 (H 2 n +H (2) n ) (2n-3) 2 = 32(268+9 ln(2)(22 ln(2)-41)) 243π -88π 81 .
It is often useful to translate summations such as the series given in the above propositions using integral identities such as

1 0 nx n-1 ln 2 (1 -x) dx = (H n ) 2 + H (2)
n . For example, we have that

1 0 K(x) ln 2 (1 -x) dx = 48 - 4π 2 3 + 32(ln(2) -2) ln(2) (2)
from Proposition 1, together with the aforementioned integral identity. We remark that such integration results often may be used to prove elegant formulas for double series involving harmonic numbers. For example, it is easily seen that m,n≥0

16 -m 2m m 2 H n (n + 1)(m + n + 2) = 48 + 32(ln(2) -2) ln(2) π - 4π 3
using (2) together with the identity ∞ n=0

x n+1 Hn n+1

= 1 2 ln 2 (1 -x).

Variations

Recall that the main techniques introduced in this article for constructing series for1 π involving harmonic numbers are based upon the evaluation of the following integrals:

1 0 x n ln(1 -x 2 ) √ 1 -x 2 dx, 1 0 x n ln 2 1 -x 2 √ 1 -x 2 dx.
Intuitively, our approach towards constructing new series for 1 π involved the integration of "variations" of the binomial expansion of √ 1 -x 2 using the above integrals (see Lemma 1 and Section 5), so that the expression Example 1 The interesting formula

π -tanh(π) √ π = ∞ n=0 2n n (-i -n) -1 2 + (i -n) -1 2 4 n
is easily seen to hold by by evaluating

∞ n=0 (-1) n x n -1 2 n sin(ln(x)) √ 1 -x = - sin(ln(x)) x -1 ,
and integrating both sides of this equation by evaluating

1 0 x n sin(ln(x)) √ 1-x dx. Example 2 It is easily seen that ∞ n=0 2n n (H -1-i-n -H -1+i-n ) 4 n (2n -1)
is equal to

4 cosh(π) Γ (1 + i)Γ 3 2 -i -Γ (1 -i)Γ 3 2 + i 5 √ π by evaluating ∞ n=0 (-1) n x n 1 2 n sin(ln(x)) -1 + x = - sin(ln(x)) √ 1 -x ,
and integrating both sides of this equality, by evaluating 1 0

x n sin(ln(x)) -1+x

dx.

Example 3 The equality and integrating both sides of this equation.

- 4i coth(π) 5 √ π = ∞ n=0 2n n 4 n (2n -1) 1 (-i -n) 1 2 - 1 (i -n) 1
We currently leave it as an open problem to further investigate variations of the techniques given in Section 2.

Conclusion

The evaluation of integrals involving complete elliptic integrals based on the evaluation of harmonic sums given by analogues of Lemma 1 may be an interesting area to explore. Also, we currently leave the natural problem of evaluating Ramanujan-like series involving higher powers of central binomial coefficients and harmonic numbers as an open problem.

Recall that one of the main tools used in our paper is based on the following parameter derivative for a beta-type integral: It would be interesting to further explore the use of other kinds of parameter derivatives of beta-type integrals in the construction of Ramanujan-like series involving 1 π and definite integrals involving complete elliptic integrals. While there are many new formulas which can be obtained by applying parameter derivatives to hypergeometric identities, the main technique indicated in Section 2.1 can be applied in a very general way to compute series involving products of the form 2n n 2 H n , as discussed in Section 1. However, we encourage the exploration of new formulas for 1 π which can be obtained from parameter derivatives of hypergeometric identities.

2 16

 2 m (2m -1)(n + 1)(m + n + 2) .

1 0

 1 ln(1 -x)K(x) dx = 8 ln(2) -8 may be used to construct yet another proof of Theorem 1.

2(- 1 ) n x 2n 1 2 n

 12 is easily seen to hold by evaluating∞ n=0 sin 2 (ln(x)) √ 1 -x 2 = sin 2 (ln(x))

u n 1 -

 1 u 2 y log 1 -u 2 du.

√ 1-x

would in some sense "cancel" with an expression similar to √ 1 -x 2 , thus yielding a relatively simple logarithmic integral being equal to π times a Ramanujan-like series involving harmonic numbers. It is natural to consider simple variations of this strategy, based on the evaluation definite integrals of the form 1 0x n e(x) √ 1-x 2 dx and similar integrals, letting e(x) denote an elementary function.

Acknowledgements The author would like to thank Dr. Jonathan Sondow for a useful discussion concerning Ramanujan-like formulas for 1 π . The author would also like to thank two anonymous reviewers for many useful comments.


Antonio Marcos Gonçalves De Lima

Noureddine Bouhaddi
email: bnoureddine.bouhaddi@univ-fcomte.fr

Alves Domingos
email: cdomingos@ufu.brdmarcelobelonsi@bol.com.br

Marcelo Rade

Marcos Antônio

Gonçalves De Lima

Domingos Alves

Marcelo Belonsi

A time-domain finite element model reduction method for viscoelastic linear and nonlinear systems

Keywords: Passive vibration control, viscoelastic damping, fractional derivative model, condensation, nonlinear vibrations

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Many authors have shown that the effective design of viscoelastic systems can be conveniently carried out by using modern mathematical models to represent the frequency-and temperaturedependent behavior of viscoelastic materials. However, in the quest for design procedures of real-word engineering structures, the large number of exact evaluations of the dynamic responses during iterative procedures, combined with the typically high dimensions of large finite element models, makes the numerical analysis very costly, sometimes unfeasible. It is especially true when the viscoelastic materials are used to reduce vibrations of nonlinear systems. As a matter of fact, which the resolution of the resulting nonlinear equations of motion with frequency-and temperature-dependent viscoelastic damping forces is an interesting, but hard-to-solve problem. Those difficulties motivate the present study, in which a time-domain condensation strategy of viscoelastic systems is addressed, where the viscoelastic behavior is modeled by using a four parameter fractional derivative model. After the discussion of various theoretical aspects, the exact and reduced time responses are calculated for a three-layer sandwich plate by considering nonlinear boundary conditions.

A time-domain finite element model reduction method for viscoelastic linear and nonlinear systems

INTRODUCTION

The use of viscoelastic materials has been regarded as a convenient strategy in many types of industrial applications, where these materials can be applied either as discrete devices or surface treatments at a relatively low cost [START_REF] Nashif | Vibration damping[END_REF][START_REF] Rao | Recent applications of viscoelastic damping for noise control in automobiles and commercial airplanes[END_REF][START_REF] Samali | Use of viscoelastic dampers in reducing wind and earthquake induced motion of building structures[END_REF]. A typical example is the application of viscoelastic constrained layers in compressors systems [START_REF] De Lima | Vibration attenuation in refrigerators using constrained viscoelastic layers: preliminary finite element modeling and experimental assessment[END_REF]. However, the behavior of viscoelastic materials presents some inherent com-plexities such as the influence of operational and environmental factors (frequency, temperature, pre-loads, etc.) [START_REF] Nashif | Vibration damping[END_REF]. Thus, several approaches have been developed for performing dynamic responses of engineering structures containing viscoelastic damping devices [START_REF] Bagley | A generalized derivative model for an elastomer damper[END_REF][START_REF] Golla | Dynamics of viscoelastic structures -A time domain finite element formulation[END_REF][START_REF] Mctavish | Modeling of linear viscoelastic space structures[END_REF][START_REF] Lesieutre | Time domain modeling of linear viscoelasticity using anelastic displacement fields[END_REF]. However, many authors have found-out that the proposed mathematical models used for performing the frequency-and temperature-dependent behavior of viscoelastic materials, based on the concept of internal variables lead to global systems of equations of motion whose numbers of degrees of freedom (DOFs) largely exceed the order of the associated undamped structures [START_REF] Balmès | Tools for viscoelastic damping treatment design: Application to an automotive floor panel[END_REF][START_REF] De Lima | An efficient modeling methodology of structural systems containing viscoelastic dampers based on frequency response function substructuring[END_REF][START_REF] De Lima | Component mode synthesis combining robust enriched Ritz approach for viscoelastically damped structures[END_REF]. As a result, if such response evaluations are made based on computations performed on the full finite element (FE) models of structures treated with viscoelastic materials composed by many thousands of DOFs, which is not rarely the case, the computational cost necessary to perform exact response evaluations during iterative processes such as structural optimization and uncertainty propagation, can become prohibitive, even unfeasible [START_REF] De Lima | Component mode synthesis combining robust enriched Ritz approach for viscoelastically damped structures[END_REF][START_REF] Zghal | Reduced-order model for non-linear dynamic analysis of viscoelastic sandwich structures in time domain[END_REF]. Moreover, the increased dimension of the analytical models can preclude their use for active control and nonlinear vibration computations, in which time-domain analyses must be performed.

These drawbacks can be circumvented by combining directly the so-called model-reduction techniques in an attempt to reduce the order of the FE model while preserving its capability to represent the dynamic behavior of viscoelastic systems. However, it must be reminded that an inherent difficulty of any model reduction procedure applied to viscoelastic systems lies on the fact that the viscoelastic stiffness matrix depends on frequency and temperature. Another difficulty is to choose the projection bases in order to reduce truncation effects [START_REF] De Lima | Component mode synthesis combining robust enriched Ritz approach for viscoelastically damped structures[END_REF].

One possible strategy is to generate a set of Ritz vectors capable of representing with satisfactory accuracy the structural behaviour under structural modifications. For example, [START_REF] Balmès | Tools for viscoelastic damping treatment design: Application to an automotive floor panel[END_REF] studied the possibility of using a constant Ritz basis to create parametric families of reduced models developed for viscoelastic structures to be applied in the frequency-domain. [START_REF] Bouazzouni | Selecting a Ritz basis for reanalysis of frequency response functions of modified structures[END_REF] developed an optimal method to construct additional vectors by using the dynamic behavior of undamped structures before modifications. De [START_REF] De Lima | An efficient modeling methodology of structural systems containing viscoelastic dampers based on frequency response function substructuring[END_REF] have proposed the use of a frequency-and temperature-independent reduction basis of viscoelastic systems based on the construction of residual static vectors taking into account the external loads and the viscoelastic damping forces.

The disadvantage of such approaches lies in the fact that a basis of reduction composed by a large number of residual static vectors is obtained, thus augmenting the computation effort involved in the condensation. Also, the proposed condensation strategies are restricted to the frequency-domain analysis, since the complex modulus approach is used to represent the viscoelastic dynamic features. More recently, [START_REF] Zghal | Reduced-order model for non-linear dynamic analysis of viscoelastic sandwich structures in time domain[END_REF] have proposed a model reduction method of nonlinear dynamic analysis of viscoelastic systems in time domain using the well-known Golla-Hughes-McTavish (GHM) model [START_REF] Golla | Dynamics of viscoelastic structures -A time domain finite element formulation[END_REF][START_REF] Mctavish | Modeling of linear viscoelastic space structures[END_REF]. However, it leads to augmented global systems of equations of motion whose numbers of DOFs largely exceed the order of the associated undamped structures. This fact motivates the proposition of a time-domain condensation method based on the use of constant enriched reduction bases, which take into account a priori information of the viscoelastic damping forces and the local perturba-tions to be applied in the transient analysis of viscoelastic systems. The perturbations can be applied related to either structural viscoelastic modifications or local nonlinearities.

As for the viscoelastic models, among the widely used mathematical representations accounting for the typical dependence of the viscoelastic properties with respect to frequency and temperature [START_REF] Bagley | A generalized derivative model for an elastomer damper[END_REF][START_REF] Golla | Dynamics of viscoelastic structures -A time domain finite element formulation[END_REF][START_REF] Mctavish | Modeling of linear viscoelastic space structures[END_REF][START_REF] Lesieutre | Time domain modeling of linear viscoelasticity using anelastic displacement fields[END_REF], in the present study, the proposed time-domain condensation approach is specially adapted for viscoelastic systems in which the viscoelastic behavior is modeled by using the four-parameter fractional derivative model (FDM) originally proposed by [START_REF] Bagley | A generalized derivative model for an elastomer damper[END_REF] and modified by [START_REF] Galucio | Finite element formulation of viscoelastic sandwich beams using fractional derivative operators[END_REF]. Although the association of this viscoelastic model in a FE discretization has been largely used in previous studies devoted to linear vibrating systems, the main contribution intended for the present study is the proposition of a straightforward time-domain condensation strategy to reduce linear and nonlinear viscoelastic systems, which requires specially adapted time-domain analytical and numerical resolution procedures.

In the remainder, after the presentation of various aspects related to the theoretical foundations, the description of numerical applications composed by a three-layer sandwich plate structure with embedded viscoelastic materials, subjected to linear and nonlinear structural modifications, demonstrates the effectiveness of the proposed condensation strategy.

REVIEW OF THE FE MODELING OF PLATES WITH VISCOELASTIC CONSTRAINING LAYERS

In this section, the model of a moderately thin three-layer sandwich plate FE, which can be frequently found, for example, in aerospace systems, is summarized based on the original developments made by [START_REF] Khatua | Bending and Vibration of Multilayer Sandwich Beams and Plates[END_REF] and implemented by [START_REF] De Lima | Sensitivity analysis of viscoelastic structures[END_REF]. The inclusion of the frequency-and temperature-dependent behavior of the viscoelastic material is made by using the so-called Elastic-Viscoelastic Correspondence Principle [START_REF] Nashif | Vibration damping[END_REF], according to which, the structural matrices are first generated for specific types of finite elements (rods, beams, plates, etc.) assuming that the longitudinal modulus and/or the shear modulus (according to the stress state assumed) are constant (independent on frequency and temperature). Then, after the finite element matrices are constructed, the frequency-temperature dependency of those moduli can be introduced according to the complex modulus approach combined with the Frequency-Temperature Superposition Principle [START_REF] De Lima | An efficient modeling methodology of structural systems containing viscoelastic dampers based on frequency response function substructuring[END_REF].

Figure 1 depicts a rectangular element formed by an elastic base-plate (1), a viscoelastic core (2) and an elastic constraining layer (3). This element contains four nodes and seven DOFs per node, representing the in-plane displacements in the middle plane of the base-plate in directions x and y (denoted by 1 u and 1 v , respectively), the in-plane displacements of the middle plane of the constraining layer in directions x and y (denoted by 3 u and 3 v , respectively), the transverse displacements, w , and the cross-section rotations about x and y , denoted by x  and y  , respectively.

In the development of the theory, the following assumptions are adopted: (i) all the materials involved are homogeneous and isotropic and present linear behavior; (ii) normal strains in direction z are null for the three layers; (iii) the elastic layers (1) and (3) are modeled according to Kirchhoff's theory; (iv) for the viscoelastic core, Mindlin's theory is adopted (transverse shear is included); (v) the cross-section rotations, x  and y  , are the same for the elastic layers; (vi) the transverse displacement, w , is the same for all the layers. These assumptions have been considered by many authors as being adequate for the modeling of thin panels, as is the case of the structures addressed herein [START_REF] Austin | Variations on modeling of constrained-layer damping treatments[END_REF]. Moreover, previous studies carried-out by the authors demonstrated satisfactory correlation between model predictions and experimental results [START_REF] De Lima | Finite element modelling and experimental characterization of beams and plates treated with constraining damping layers[END_REF].

The discretization of the displacement fields within the element is made by using bilinear interpolation functions for the displacements in the middle plane of the elastic layers in directions x and y , and a cubic interpolation function for the transverse displacement, according to the relation,

  ( , , ) ( , ) ( )
e x y t x y t  u N u
, where ( , )

x y N is the matrix formed by the shape functions, and

( ) ( ) e t  u 1 1 3 3 T i i i i i i i x y u v u v w       
 with i  1 to 4 designates the vector containing the nodal variables as functions of time. The strain-displacement relations in vector form,

   ε x, y, z,t       e t B u x, y, z
, are used and the resulting strains for the elastic layers, k  ε

( ) ( ) ( ) T k k k x y xy            1, 3 k 
, and for the viscoelastic core, 2  ε

(2) (2) (2) (2) (2) T x y xy xz yz        
 , are generated. Thus, based on the stress-states assumed for each layer and the corresponding stress-strain relations, following standard analytical developments based on variational approaches, the strain and kinetic energies of the three-layer sandwich plate FE are formulated and the elementary mass and stiffnesses matrices are obtained as follows:

      3 1 0 0 , , d d b a e T k k k x y h x y x y y x         Μ N N (1.a)       1,3 0 0 0 , , , , d d d k h b a e T e k k k k z x y x y z x y z y x z          Κ B C B (1.b)         2 2 2 2 0 0 0 , , , , d d d v h b a e T z x y x y z x y z y x z          Κ B C B (1.c)
where

k C   1, 3 k 
is the isotropic elastic material properties matrix, and

2 v    C
C contains the frequency-and temperature-dependent material properties for the viscoelastic core. Matrix

 

x, y, z B is formed by differential operators appearing in the strain-displacement relations for each layer, as detailed in de [START_REF] De Lima | Sensitivity analysis of viscoelastic structures[END_REF]. Also, b and a designate, respectively, the dimen-sions of the plate element in directions x and y, and k h and k  are the thickness and the mass density of the k -th layer, respectively.

According to the theory of the sandwich plate finite element summarized above and assuming that the Poisson ratio is independent from frequency and temperature in such a way that the frequency-and temperature-dependent longitudinal and transverse moduli are related to each other through the relation,

  2 1 v v E G     
, the design parameters of mass and stiffness of each layer can be factored-out of the elementary matrices by uncoupling membrane, bending and shear effects as:

                      1 1 1 1,1 1 1 1,1 2 2 2 2 2,0 2 2 1 2,1 2 2 1 2,2 2 2 2,0 3 3 3 3,1 3 3 3,1 e e e
h h d h d h h h                 M M M M M M M M M M M
(2.a)

                                    3 1 1 1 1,1 1 1 1,1 2 2 2 2 2,0 1 2,1 1 2,2 2,0 2 2,1 2 2,2 2 2 2,0 3 2,1 3 2 2,2 2 3 3 3 3 3,1 3 3 3,1
E h E h E h d d d d E d d d h E h E h                   K K K K K K K K K K K K K K K K (2.b)
where

3 1 2 3 2 d h h h    and k E ( 1, 3 k 
) represents the longitudinal moduli of the elastic layers. Also, subscripts m, b and s indicate, respectively, the membrane, bending and shear effects in the structural matrices.

It should be noted that, in Eqs.

(2), the matrices appearing in the right-hand side are constant sparse matrices which are independent from the design parameters. It is clearly perceived how the use of such expressions facilitates a variety of finite element computations, since it enables to account for structural modifications and/or parametric uncertainties in the values of the design parameters in a straightforward way during iterative processes. Also, it facilitates, to a large extent, the evaluation of the sensitivities of the responses with respect to the design parameters.

From the elementary matrices computed for each element of the FE mesh, and neglecting other forms of damping, the elementary equations of motion are given as: (1)(2)

                  e e e e e
e v e t t t

           M u K K u f (3)
where the stiffness matrices,

      1 3 e e e e   Κ Κ Κ and     2 e e v    Κ Κ
give, respectively, the contributions of the purely elastic and viscoelastic parts of the damped structure, and

        1 2 3 e e e e    M M M M
includes the contribution of the base-plate, the viscoelastic core and the constraining layer to the inertia matrix. Finally,     e t u and     e t f are, respectively, the vectors of the generalized displacements and external loads.

Cleary, Eq. ( 3) must be used to construct the global equations of motion, accounting for the node connectivity of the discretization mesh, using standard FE assembling procedures.

Fractional derivative model incorporated into FE matrices

The dependency of the viscoelastic stiffness matrix on frequency and temperature is a consequence of the dependency of the material moduli with respect to these factors. A variety of mathematical models has been developed to represent the viscoelastic behavior and have been shown to be suitable to be used in combination with the FE discretization.

In this paper, as the interest is confined to a proposition of a time-domain reduction procedure of viscoelastic systems, the Fractional Derivative Model (FDM) initially proposed by [START_REF] Bagley | A generalized derivative model for an elastomer damper[END_REF] is used in combination with the Grünwald discretization technique [START_REF] Galucio | Finite element formulation of viscoelastic sandwich beams using fractional derivative operators[END_REF] to approximate the fractional operator appearing in the following one-dimensional stress state constitutive equation assumed for the viscoelastic material:

      0 d d t E E t t E t             (4)
where 0 E and E  are, respectively, the relaxed or static modulus, and non-relaxed or highfrequency limit value of the modulus,  is the relaxation time, and  represents the fractional order of the time derivative  

0 1    .
One important aspect regarding the use of the following complex modulus function,

        0 1 v E E E i i         
, obtained by applying the Fourier transform to Eq. ( 4) is the identification of the four parameters ( 0E , E  ,  ,  ) from experimental data-sheets provided by the manufactures containing the material storage modulus and loss factor as functions of frequency and temperature. Thus, from the complex modulus definition, the determination of the values of the material parameters can be carried out by formulating an optimization problem in which the objective function represents the difference between the experimental data points and the corresponding model predictions in the frequency band of interest for a fixed temperature value [START_REF] Galucio | Finite element formulation of viscoelastic sandwich beams using fractional derivative operators[END_REF]. As a result, the designer can obtain the storage modulus and the loss factor at any given temperature into a frequency band of interest, as illustrated in Fig. 2 for the particular material considered in the present study. Also, it must be emphasize that Eq. ( 4) was obtained by introducing the following internal variable as a strain function,      

t t t E      
, in the standard one-dimensional stress state constitutive equation for linear viscoelastic materials. As a result, Eq. ( 4) contains only one fractional derivative term,

  d t d t   
, instead of two presented in the classical constitutive equation, and can be approximated by the following Grünwald definition,

  d t d t       1 0 p n j j t A t j t         
, noting that 1 1 A  , to generate the following discretized form of the anelastic strain:

        0 1 1 1 p n j j E E t c t c A t j t E               (5)
where

  c t         , t t n   is the time step, p n n 
is the number of discretization points, and

1 j
A  represents the Grünwald coefficients given by the recurrence formula,

1 j A     1 j j A j    . At this point, the term       e v e t  K u
appearing in equations of motion (3) can be expressed as,

      e v e t   K u     , , , , , d d d T z y x x y z x y z t x y z    σ B ,
where for a general stress-strain relation, the vector  

, , , x y z t σ can be written in terms of the anelastic strain (5) taking into account the strain relation,      

t t t E      
, as follows:

      0 0 1 1 0 0 , , , 1 , , , , , , p n j
j E E E x y z t E c x y z t c A x y z t j t E E                                  σ ε ε (6)
By combining Eq. ( 6) with the elastic and anelastic strain-displacement relations,

        e t  ε B u x, y, z,t x, y, z and         e t  B ε u x, y, z,t
x, y, z , the resulting modified elementary equation of motion takes the following form:

                e e e v e v e v t t t j t        K u K u f (7)
where

    0 0 1 e e v v E E c E                  K K and           1 1 0 p n e e v v j e j E t j t c A t j t E           f K u .
Upon introduction of Eq. ( 7) into Eq. ( 3) and after matrix assembling, the global equations of motion of the viscoelastic system containing N DOFs can be written as:

          e v v t t t t j t          M u K K u f f (8)
where

 

v t j t    f is a load vector, which depends on the anelastic displacements,   t u , and the Grünwald coefficients. These later account for the fading memory of viscoelastic materials.

It is interesting to highlight that Eq. ( 8) represents a time-domain model of structures containing viscoelastic materials, which can be interpreted as the result of incorporating modifications in the original elastic stiffness of the base-structure associated to viscoelastic materials, which is, in fact, independent of time, and by adding viscoelastic time-varying damping forces to the original external loads.

By comparing this modeling approach to other alternatives that have been proposed based on the adoption of particular representations for the frequency-and temperature-dependent behavior of the viscoelastic materials [5][6][7][8]12], the proposed strategy does not require transformation of the equations of motion into a state-space form to generate frequency-and temperature-independent state matrices.

Based on Eq. ( 8), a time-domain condensation method based on the use of a frequency-and temperature-independent reduction bases is addressed next.

TIME-DOMAIN CONDENSATION OF VISCOELASTIC SYSTEMS FOR DESIGN PROCEDURE

In the case of complex structures of industrial interest, FE models are usually constituted by a large number of DOFs (hundreds of thousand or even millions). In such cases, it becomes practically impossible to compute the time-domain response directly from Eq. ( 8), owing to high computation times and storage memory required. This fact motivates the use of model reduction procedures, which aim at reducing the model dimensions (and the associated computational burden), while keeping a reasonable predictive capability of the numerical models. This can be done based on the assumption that the exact response, given by the resolution of Eq. ( 8), can be approached by a projection on a reduced vector basis as:

    0 t t  u T u (9)
where 0

N NR C   T
is the transformation matrix formed column-wise by vectors,   ˆNR t C  u is the vector of generalized coordinates, and NR N  is the number of vectors in the basis. Once associated to the transformation expressed in Eq. ( 9), the equations of motion (8) can be rewritten as follows:

          ê v v t t t t j t          M u K K u f f (10)
where

0 0 ˆT e e  K T K T , 0 0 ˆT v v    K T K T , 0 0 ˆT  M T MT ,     0 ˆT v v t j t t j t        f T f and   ˆt  f   0 T t T f
are the reduced matrices and vectors.

For models of structures containing viscous or structural damping, it is relatively common to use a constant projection basis formed by a subset of eigenvectors of the associated conservative structure, as the mass and stiffness matrices are invariant [START_REF] Nashif | Vibration damping[END_REF]. However, for systems containing viscoelastic materials, the selection of the reduction basis is more delicate as this condition does not hold. Owing to the dependence of the stiffness matrix with respect to frequency and temperature, the reduction basis should be able to represent the changes of the dynamic behavior as frequency and temperature vary [START_REF] De Lima | Component mode synthesis combining robust enriched Ritz approach for viscoelastically damped structures[END_REF]. This motivates the proposition of a special time-domain condensation method based on the use of a frequency-and temperature-independent enriched reduction basis which takes account a priori information of the viscoelastic damping forces. These static responses are computed by using the tangent viscoelastic stiffness matrix computed in the low frequency range by considering the real value of the modulus, 0 E , as follows:

0 0 v E    K K K e (11)
Hence, the nominal basis can be obtained by the resolution of the eigenvalue problem:

    0 0 1 2 0 1 0 1, , , , , i i NR
NR i N diag                  ϕ ϕ ϕ ϕ ϕ Λ K M (12)
The basis, 0 ϕ , is further enriched by introducing the following residues formed by the static displacements associated, respectively, to viscoelastic damping forces and external excitations:

0 1 0 0 v v    ϕ R K K (13.a) 1 0   R K b (13.b)
where

N f R   b
is a Boolean matrix which enables to select, among the DOFs, those in which the unit excitation forces are applied.

The residues (13.a) are interpreted as the columns of the flexibility matrix of the associated undamped system, related to the coordinates of application of the viscoelastic damping forces, which can be better understood by examining Eq. ( 8), noting that the term involving the viscoelastic behavior can be moved to the right-hand side, where it plays the role of additional forces applied to the associated conservative structure.

Thus, the enriched basis of reduction for the viscoelastic system is given as: (13)

0 0 0 v        ϕ T R R (14)

Static residual vectors accounting for local nonlinearities

Local nonlinear behavior can be frequently found in a number of real-word engineering systems such as those having mechanical connections and joints. Thus, the objective of this section is to extend the proposed reduction method to viscoelastically damped systems supported by lumped nonlinear springs. In this case, the inclusion of nonlinear effects into the global equations of motion (8) can be easily done by the concept of dyadic structural modifications [START_REF] Maia | Theoretical and experimental modal analysis[END_REF], in which the loading system exerted by nonlinear springs on the viscoelastic system at the node of attachment, nd , can be expressed in the following vector form, by assuming cubic nonlinearity:

      nd nd nd nl l nd nl nd t t t    f K u K u (15)
where  

1 1 3 3 T nd x y t u v u v w          u and   3 3 3 3 3 3 3 1 1 3 3 T nd x y t u v u v w           u
are the vectors containing the DOFs associated to the node of attachment of the linear and nonlinear contributions of the nonlinear springs in the global coordinate system.

Hence, the global system of equations of motion for the nonlinear viscoelastic system can be expressed under the following form:

            e l v v t t t t j t t             M u K K K u f f f (16)
where

T
    nl t t    f K u
. For small nonlinearities, Eq. ( 16) shows that the nonlinear behavior can be understood as a perturbation indicated by,

  t  f
, of the linear behavior, and can be view as a modification introduced on the linear viscoelastic system. As a result, the reduction method presented in the previous section can be extended to the case of nonlinear viscoelastic systems by a further enrichment of the reduction basis including normal modes related to the linear conservative associated system perturbed by the nonlinear springs. These modes can be obtained by the resolution of the eigenvalue problem:

  0 0 0 1, , i i i N         ϕ K M (17)
where

0 0 l    K K K .
Based on the fact that the vibration modes of the perturbed system does not differ appreciably from that of the non-perturbed nominal system, it is possible to assume that, i   ϕ

0 i i     ϕ ϕ
, where, 0 i  ϕ , designates the modes of the nominal conservative associated system obtained by the resolution of Eq. ( 17), and i   ϕ is the perturbation due to the local nonlinearities which can be interpreted as an external load applied on the linear model, 0 i  ϕ . Thus, a set of static responses can be generated:

1 0        ϕ K b (18)
where

N m R     b
is also a Boolean matrix which enables to select the nonlinear DOFs in which the unit forces are applied and static responses are computed, and m is the number of nonlinear DOFs.

Hence, the final condensation basis taking into account a priori knowledge of the perturbations can be expressed as follows:

0 1 m             ϕ ϕ T T (19)
The residues,   ϕ , are not necessarily of maximum rank. Thus, with the aim of obtaining a limited number of independent residual vectors, it is appropriate to select the dominating directions of this basis, which can be done by performing the Singular Value Decomposition (SVD) of T to identify its dominant singular values.

By considering the basis ( 19), the equations of motion ( 16) can be rewritten as:

          0 ˆv t t t t j t t           M u K u f f f (20)
where

ˆT  M T MT , 0 0 ˆT    K T K T ,     ˆT t t  f T f ,     ˆT v v t j t t j t        f T f and   ˆt   f   T t  T f
are the reduced matrices and vectors. As will be demonstrated later on, the reduction basis given in Eq. ( 19) is robust enough to provide accurate predictions of the reduced model under the influence of the viscoelastic treatment, external excitations and nonlinear supports.

ITERATIVE RESOLUTION PROCEDURE

Based on the modeling features of the viscoelastic behavior and the presence of local nonlinearities, the resulting nonlinear differential equation of motion (20) were numerically solved using a variant of the Newmark integration scheme [START_REF] Bathe | Finite element procedures[END_REF] performed in the MATLAB environment.

Figure 3 shows the main steps of the numerical procedure which can be summarized as follows: (i) at the beginning of the process, 0 t , after defining the operating temperature of the viscoelastic material, v T , the Newmark integration parameters,   ,   , and the variables   0 , , , E E    for FDM model obtained from curve-fitting of experimental data [START_REF] Galucio | Finite element formulation of viscoelastic sandwich beams using fractional derivative operators[END_REF], it is possible to compute the real and anelastic displacements, velocities and accelerations fields,  

0 t u ,   0 t u ,   0 t  u and   0 t  u
, respectively, and the full mass, M , and stiffnesses matrices  

, , ,

e v l nl  K K K K
taking into account the mechanical boundary conditions; (ii) next, an eigenvalue problem is performed according to Eq. ( 17) in order to extract the eigenvalue modes, 0  ϕ , of the linear conservative associated system to generate the first basis of reduction to be further enriched by the residual static vectors associated to the external excitations,

1 0    R K b , the viscoelastic damping forces, 0 1 0 0 v v      R K K ϕ , the modifications introduced on the viscoelastic zones, 1 0        K b ϕ , or by the local nonlinearities, 1 0       K b ϕ .
With the enriched basis of reduction it is possible to reduce the system in order to generate the reduced matrices, M and 0  K and vectors, û , û , f , ˆ f and  f ; (iii) next, an iteration is initiated and the loading vectors due to the external loads, f , the nonlinear perturbation, ˆ f , and the viscoelastic damping forces, ˆv  f , are evaluated at each iteration process. Also, a new set of reduced displacement fields,

APPLICATIONS TO PLATE TREATED WITH CONSTRAINED VISCOELASTIC LAYER

In this first application, numerical tests were performed using the representative system shown in Fig. 4 composed by a clamped-clamped-free-free rectangular plate made of aluminum, fully treated with a layer of 3M ISD112™ (2013) viscoelastic material constrained between the baseplate and an outer aluminum sheet. The FE model is composed of 140 elements and the viscoelastic treatment is modeled according to the three-layer sandwich plate element developed according to Section 2. The geometric dimensions in x and y directions are depicted in the same figure, and the thicknesses of the base-plate, constraining layer, and viscoelastic core are, respectively, 1.0 mm, 0.5 mm and 0.254 mm. The material properties for the base-plate and constraining layer are: Young modulus E  70×10 9 N/m 2 ; mass density   2700 kg/m 3 ; Poisson ratio,   0.3. For the 3M™ ISD112,   950kg/m 3 ,   0.49, and the parameters of the FDM model can be obtained by performing a curve-fitting procedure taking into account the viscoelastic material data provided in Fig. 2 in terms of storage modulus and loss factor, for a given temperature of the viscoelastic material.

The plate is subjected to a harmonic force of the form,    

0 0 sin 2 F t F f t  
, where 0 F  10 N is the amplitude of the excitation and 0 f is the excitation frequency, and to a transverse triangular impulse loading applied to point I indicated in Fig. 4 The interest now is to evaluate the enriched basis of reduction by using the static residues associated to the external excitations and viscoelastic damping forces. To verify the direct condensation procedure, one considers the following nominal basis:

1 0        T R ϕ ( 19
eigenvectors, plus one residual vector computed by the Eq. (13.b));

0 2 0 v        T R R ϕ ( 19
eigenvectors, one residual computed by the Eq. (13.b), 14 residual vectors computed according to the Eq. (13.a) after SVD filtering). The residues 0 v R were computed based on the largest singular values, for which the relation, max i    1x10 6 , for i  1 to 14 holds. The interest in examining these situations is to quantify the improvement entailed by the residual vectors associated to the viscoelastic damping forces and to demonstrate the capability of the time-domain reduction modeling procedure to accommodate such design variants. Figure 6 shows the normalized time responses obtained by using the two nominal bases, as compared to the reference displacements computed by using the exact system. It can be clearly seen that the accuracy is improved upon enrichment of the reduction basis by the inclusion of residual vectors accounting for the static residues associated to the viscoelastic damping forces for both input signals, demonstrating that the use of these residues are sufficient to represent with accuracy the time-domain behavior of viscoelastically damped structures. As a complementary demonstration of the utility of the proposed condensation strategy in the analysis of viscoelastic systems, Fig. 7 shows the amplitudes of the frequency response functions (FRFs) of the sandwich plate computed by using the two nominal bases for the impulse loading condition, as compared to the amplitudes of the FRF computed for the exact system. Again, Fig. 7(b) confirms that the use of first order residues associated with viscoelastic damping forces are sufficient to represent with satisfactory accuracy the dynamic behavior of the viscoelastic system. Moreover, these results provide a sense of the efficiency of the surface viscoelastic treatment in mitigating the transverse vibrations of the plate. T and 2 T .

Structural modifications

In this second application the utility of the time-domain reduction strategy in the analysis of modified viscoelastic systems is shown. The interest is to evaluate the robustness of the nominal basis further enriched to account for small perturbations introduced into the nominal model, according to Eq. ( 19). The modifications considered consist in increasing the thicknesses of the viscoelastic and constraining layers of the nominal system in 80%, as indicated in Fig. 4, and the exact time-domain responses and the amplitudes of the FRFs of the perturbed system was computed using the modified FE model, as shown in Figure 8. It can be noted that the dynamic behavior of the perturbed system does not differ appreciably from that of the nominal system. In Fig. 9, the responses of the perturbed exact system in the time and frequency domains are compared to the counterparts computed by using the nominal reduction basis 2 T containing 34 vectors, without further enrichment accounting for structural modifications. The observed differences for both domains lead to conclude that this basis is not capable of representing accurately enough the changes of the dynamic behavior induced by the structural viscoelastic modifications, since the basis 2 T represents with accuracy the behavior of the nominal system, as shown in Fig. 6. Thus, in a number of real-word engineering applications of passive constraining damping layers, this reduction may be found to be insufficient to provide the necessary representation of the dynamic behavior of viscoelastic materials associated to structural variations. In the same figure, it is possible to compare the dynamic responses of the perturbed exact system to the counterparts computed by using the basis 3

2         T R
Τ containing 41 vectors (including 7 SVD-filtered residual vectors) associated to the structural modifications. This time, one can observe a very satisfactory agreement between the amplitudes of the dynamic responses of both models. This leads to conclude that the reduction basis 3 T is robust enough to represent the response of the perturbed viscoelastic system.

Modifications due to local nonlinearities

In this section, it is considered the same system shown in Fig. 4 subjected to non-linear boundary conditions through the application of transverse non-linear springs on the right edge of the plate located at x  0.30 m. The non-linear behavior of the springs are computed by the Eq. ( 15), where the linear, l l l k   K K , and non-linear, nl nl nl k   K K , stiffnesses are evaluated by assuming the values l k  100 N/m and nl k  1x10 9 N/m 3 , where l K and nl K are diagonal identity matrices. For the harmonic loading condition, 0 F  30 N is the amplitude of the excitation and the excitation frequency, 0 f  24 Hz, is close to the first flexural mode of the conservative associated system in which the stiffness is computed by the relation, 0

0 v l E      e K K K
K . The interest is to evaluate the enriched basis of reduction by using the static residues associated to the non-linear perturbations by considering the following basis:

0 2 0 v        T R R ϕ ( 19
eigenvectors, one residual vector computed by the Eq. (13.b), 14 residues computed according to the Eq. (13.b) after SVD filtering);

3 2          T T
ϕ (the basis 2 T and 11 vectors computed by the Eq. ( 18)).

Figure 10 shows that a reasonable agreement between the dynamic responses of both models is obtained by using the basis, 3

T . The agreement is less satisfactory for the reduction basis, 2 T , leading to conclude that the static residues associated to the local non-linearities is robust enough to represent the response of the non-linear viscoelastic system. Clearly, it should be reminded that this kind of reasoning is strictly valid for small non-linear damped behaviors as, depending on the amplitude of the excitation, the use of a constant enriched basis can only provide rough approximations of exact non-linear behaviors, as shown in Fig. 11. It can be noted the degree of influence of the non-linear behavior on the prediction: the larger the non-linear behavior with respect to a given excitation, the larger the influence of it on the accuracy of the reduced model. In this case, the results could possibly be improved by updating the reduction basis, 3 T , at each iteration cycle to account for the non-linear variations. In general, it results in a prohibitive computational time, since a set of updated eigenproblems must be solved. Another strategy is the so-called Combined Approximations Approach (CA) [START_REF] Kirsch | Nonlinear and dynamic structural analysis using combined approximations[END_REF] developed originally for linear and non-linear undamped reanalysis. These aspects should be considered in future works.

CONCLUDING REMARKS

A robust time-domain condensation procedure intended to be used for dealing with linear and non-linear viscoelastically damped systems was suggested and evaluated. The original aspects of the procedure lies in the adaptation of the concept of robust condensation, initially developed for linear undamped structures in the frequency-domain, for systems containing viscoelastic materials in the time-domain, and the extension of the proposed methodology to non-linear structures incorporating viscoelastic materials with the aim of vibration attenuation. This approach allows design parameters of a superelement to be modified, for example in the context of linear and nonlinear reanalysis such as an optimization and/or model updating processes and structural damage analysis, without the necessity to perform a complete superelement analysis at each point in parameter space. Thus, the improved model reduction transformation can be prepared in advance and then used directly during the iteration processes to avoid the exact recalculation of the modified zones.

At the present time, the proposed time-condensation process is not constrained solely to the dynamic applications discussed herein, and can already be used to approximate reanalysis of linear and non-linear behaviors to viscoelastically damped structures in control technology and optimization and structural damage analyses where the iterative solutions consist of repeated analyses followed by redesign steps. However, it must be reminded that depending on the structural local perturbations introduced on the viscoelastic zones and the non-linear behavior to be considered, the proposed time-condensation are no longer valid and the reduction basis must be adapted to account for the local perturbations during the iteration cycles.

Finally, it is important to pointed-out that that the efficiency of the proposed timecondensation compared with the complete analysis of the modified viscoelastic systems can be measured by the accuracy the dynamic responses, as demonstrated by the obtained results, and computational effort. It was found that the CPU effort for the dynamic reanalysis of the modified linear viscoelastic systems is reduced by more than 64%, compared with the complete analyses of the modified design. For the viscoelastic systems incorporating local non-linearities the total computational effort can be reduced by more than 96% for both cases investigated herein.

Figure 1 :

 1 Figure 1: Schematic representation of three-layer sandwich plate element.

Figure 2 :

 2 Figure 2: (a) Storage modulus and (b) loss factor for different temperatures for the 3M™ ISD112.

  finite element mass and stiffness matrices, and symbol  indicates matrix assembling.   t u is the vector of global DOFs and   t f is the vector of the generalized external loads.

  are, respectively, the diagonal stiffnesses matrices corresponding to the linear and nonlinear contributions, th columns of the identity matrices of order N , and

R

  , performed according to the Newmark integration scheme, and taking into ac-count the reduced matrices, M and 0  K ; (iv) in the next step, a new set of reduced velocities, , are generated. The iterative process is stopped when the maximum number of iterations is reached with a tolerance criterion supported by the Newmark integration scheme. Also, the iterative process is the same for linear viscoelastic systems, where 0

Figure 3 :

 3 Figure 3: Flowchart of the numerical procedure used to compute the system response.

  . The numerical tests were performed in the time interval, an arbitrarily chosen constant time step, t   0.01 ms, and the computations consisted in obtaining the normalized transverse displacements,

Figure 4 :

 4 Figure 4: Sketch of the sandwich plate with the mechanical boundary conditions.

Figure 5 :

 5 Figure 5: Transverse displacements of the plate with and without surface viscoelastic treatment due to the applied harmonic (left) and impulse (right) loadings.

Figure 6 :

 6 Figure 6: Transverse displacements of the reference and reduced systems by using the nominal basis T 1 and T 2 .

Figure 7 :

 7 Figure 7: Amplitudes of the FRFs of the reference and reduced systems by using the nominal basis 1 T and 2 T .

Figure 8 :

 8 Figure 8: Time-domain responses and FRFs amplitudes for the nominal and perturbed exact systems.

Figure 9 :

 9 Figure 9: Time-domain responses and FRFs amplitudes for the exact and reduced systems by using the basis T 2 and T 3 .

Figure 10 :

 10 Figure 10: Time-domain responses of the reference and reduced systems by using the basis 2 T and 3 T for 30 N.

Figure 11 :

 11 Figure 11: Time-domain responses of the reference and reduced systems by using the basis 2 T and 3 T for 50 N and 100 N.

Acknowledgements

The authors are grateful to the Minas Gerais State Agency FAPEMIG for the financial support to their research activities and the Brazilian Research Council -CNPq for the continued support to their research work, especially through research projects 308310/2010-1 (A.M.G. de Lima) and 482963/2009-3 (D.A. Rade). The authors express their acknowledgements to the National Institute of Science and Technology of Smart Structures in Engineering, jointly funded by CNPq, CAPES and FAPEMIG.


Khaoula Chikhaoui
email: khaoula.chikhaoui@femto-st.fr

Diala Bitar
email: diala.bitar@femto-st.fr

Najib Kacem
email: najib.kacem@femto-st.fr

Noureddine Bouhaddi
email: noureddine.bouhaddi@femto-st.fr

ROBUSTNESS ANALYSIS OF THE COLLECTIVE DYNAMICS OF NONLINEAR PERIODIC STRUCTURES UNDER PARAMETRIC UNCERTAINTY

Keywords: uncertainty, nonlinearity, periodic structure, collective dynamics

In order to ensure more realistic design of nonlinear periodic structures, the collective dynamics of a coupled pendulums system is investigated under parametric uncertainties. A generic discrete analytical model combining the multiple scales method, the perturbation theory and a standing-wave decomposition is proposed and adapted to the presence of uncertainties. These uncertainties are taken into account through a probabilistic modeling implying that the stochastic parameters vary according to random variables of chosen probability density functions. The proposed model leads to a set of coupled complex algebraic equations written according to the number and positions of the uncertainties in the structure and numerically solved using the time integration Runge-Kutta method. The uncertainty propagation through the established model is finally ensured using the Latin Hypercube Sampling method. The analysis of the dispersion, in term of variability of the frequency and amplitude intervals of the multistability domain shows the effects of uncertainties on the stability and nonlinearity of a three coupled pendulums structure. The nonlinear aspect is strengthened, the multistability domain is wider, more stable branches are obtained and thus the multimode solutions are enhanced.

INTRODUCTION

Under the hypothesis of perfect periodicity, many engineering researches provided interesting insights in the dynamic analysis of nonlinear periodic structures. For instance, a straight-forward perturbation analysis was applied by Boechler et al. [1] to investigate amplitude-dependent dispersion of a discrete one dimensional nonlinear periodic chain with Hertzian contact. Moreover, the multiple scales approach was combined with the secular perturbation theory, in a discrete model, by Lifshitz and Cross [2] to calculate the responses of a coupled array of nonlinear oscillators under parametric excitation and recently by Bitar et al. [3] to investigate the collective dynamics of a periodic structure of coupled nonlinear Duffing-Van Der Pol oscillators under simultaneous parametric and external excitations. Multiple scales perturbation approach was also used by Gutschmidt and Gottlied [4] in a continuum-based model to investigate the dynamic behavior of an array of N nonlinearly coupled microbeams and by Manktelow et al. [5] to study the wave interactions in diatomic chain with two degrees of freedom per unit cell. Recently, Andreassen et al. [6] studied the wave interactions in a periodically perforated plate through the two-1 dimensional dispersion characteristics, group velocities and internal resonances investigation.

Nevertheless, far from the mathematical idealization due to the perfect periodicity hypothesis, imperfections, which can be due to material defects, structural damage, ageing, fatigue, manufacturing processes, etc., and which reflect the nature of real engineering systems, can perturb the perfect arrangement of cells in a structure and change significantly the dynamic behavior from the predictions done under perfect periodicity. In the literature, Kissel [7], for instance, investigated the effects of the disorder in one-dimensional periodic structure theoretically and numerically using Monte Carlo (MC) simulations and proved that the disorder causes wave attenuation and pronounced spatial localization of normal modes at frequencies near the stopbands of the perfectly periodic associated structure. Statistical investigation of the effect of disorder on the dynamics of multidimensional coupled periodic structures, using the MC method, was carried out by Cha et al. [START_REF] Cha | A statistical investigation of the forced response of finite\ nearly periodic assemblies[END_REF]. Moreover, to study the effects of the randomness of flexible joints on the free vibrations of simply-supported periodic large space trusses, Koch [START_REF] Koch | Structural Dynamics of Large Space Structures Having Random Parametric Uncertainties[END_REF] combined an extended Timoshenko beam continuum model, MC simulations and first-order perturbation method.

Nonlinear coupled pendulums arrays are widely used in mechanical engineering and have been purpose of several researches in the literature. For instance, Jallouli et al. [START_REF] Jallouli | Nonlinear dynamics of a 2D array of coupled pendulums under parametric excitation[END_REF] investigated the nonlinear dynamics of a two-dimensional array of coupled pendulums under parametric excitation. Moreover, in the context of imperfections, Tjavaras and Triantafyllou [START_REF] Tjavaras | Non-linear response of two disordered pendula[END_REF] investigated numerically the effect of nonlinearities on the forced response of two disordered pendula coupled through a weak linear spring. Hai-Quing and Yi [START_REF] Hai-Quing | Parametrically driven solitons in a chain of nonlinear coupled pendula with an impurity[END_REF] developed a discrete theoretical model based on the envelope function approach to study analytically and numerically the effect of mass impurity on nonlinear localized modes in a chain of parametrically driven coupled pendula. The influence of impurities on the envelope waves in a driven nonlinear pendulum chain has been investigated numerically under a continuum-limit approximation in [START_REF] Hu | Interactions between impurities and nonlinear waves in a driven nonlinear pendulum chain[END_REF] and experimentally in [START_REF] Zhu | Experiments on the interactions between impurities and solitary waves in lattice model[END_REF].

Extending the framework to the dispersion analysis of nonlinear periodic structures dynamics under parametric uncertainty is a complex challenge and the main aim of this work is to deal with. Therefore, based on the work performed by Diala et al. [START_REF] Bitar | Multi-mode solutions in a periodic array of coupled nonlinear pendulums under primary resonance[END_REF] and under the hypothesis of small imperfections which implies that the collective dynamics and the localization phenomena due to the weak coupling of components are not destroyed, the present paper is focused on the collective dynamics analysis of a stochastic periodic nonlinear structure composed of coupled pendulums. A generic discrete analytical model combining the multiple scales method, the perturbation theory and a standing-wave decomposition is proposed and adapted to the presence of uncertainties. The latter are taken into account, in a probabilistic framework, as parametric uncertainties modeled by random variables. The obtained analytical model leads to a set of coupled complex algebraic equations written according to the number and positions of the uncertainties in the structure and then numerically solved using the time integration Runge-Kutta method. To propagate uncertainties through the established model, the statistical Latin Hypercube Sampling method [START_REF] Helton | Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems[END_REF] is used. To illustrate the stochastic effect on the dynamic behavior of a three coupled pendulums structure, when the length of the first pendulum is supposed to be uncertain, a statistical analysis of the frequency responses dispersions, both in modal and physical coordinates, is carried out. Figure 1 illustrates a generic structure of coupled pendulums of identical mass and viscous damping coefficient generated by the dissipative force acting on the supporting point of each one. The pendulums are coupled by linear springs of stiffness and are subject to an external excitation cos Ω each one. The inclination angle from the equilibrium position quantifies the rotational displacement of the pendulum. Applied boundary conditions are such as the pendulums labeled 0 and 1 are fixed so that 0. The periodicity of the structure is broken by the presence of pendulums containing parametric uncertainties, which can, for instance, be the pendulum length as illustrated in Fig. 1.

MODEL

Equation of motion

Applying the Lagrange approach leads to the equation of motion of the pendulum of the form

, cos Ω (1)
where , , , when the pendulum is deterministic and , , ,

1 1 f cos(Ω t) l l s l k k k m m m 2
when, for instance, the length of the pendulum, of stochastic displacement , is uncertain. Since the linear coupling between the pendulums is supposed to be very weak and under the hypothesis of small imperfections, each frequency is supposed to be equal to the eigenfrequency ( ).

The linear coupling term , depends on the position of the stochastic pendulums in the chain. In fact, considering the case where the stochastic pendulums are not adjacent, four different configurations are distinguished: a. If the concerned pendulum is deterministic as well as its neighbors, , 2 ; b. If the concerned pendulum is deterministic but the previous one is stochastic, , 2 ; c. If the concerned pendulum is deterministic but the following is stochastic, , 2 ; d. If the stochastic pendulum is concerned, the deterministic displacement is replaced by the stochastic one, , in Eq. (1) such as

, cos Ω (2)
and , 2 , in this case.

Numerical Solving

The multiple scales method consists on replacing the single time variable by an infinite sequence of independent time scales ( ), where is a dimensionless parameter assumed to be small, and eliminating secular terms in the fast time variable . The scaled equation of motion of the pendulum, with respect to each configuration mentioned above, is assumed as follows

, cos Ω ( 3 )
where the excitation frequency Ω is expressed as Ω with being the detuning parameter. The solution of Eq. (3) can generally be given by a formal power series expansion ∑ . Up to the first order, one obtains an equation of the form

2 , (4)
For simplicity, Eqs. (3-4) are illustrated with respect to Eq. (1). The same methodology is applied to Eq. (2). The displacement can be expanded as a sum of standing wave modes with slowly varying amplitudes. Taking into account the boundary conditions 0, the standing wave modes are

u sin n ; , 1 … (5)
The displacement of the pendulums is thus expressed as

∑ sin n . . ( 6
)
when it is deterministic and

sin n . . (7)
if it contains uncertainty. Substituting Eq. (6) into Eq. ( 4) leads to equations of the form [START_REF] Cha | A statistical investigation of the forced response of finite\ nearly periodic assemblies[END_REF] where the secular terms should be equated to zero to satisfy the condition of solvability of the multiple scales method. An equation similar to Eq. ( 8) in terms of is obtained when substituting Eq. (7) into Eq. (4).

Projecting the response on the standing-wave modes implies multiplying all terms by sin and summing over which leads to the generic complex equation of the amplitude

2 2 cos * * sin n cos n sin ∑ * , , Δ , exp ∑ sin n 0 (9)
where Δ , is the delta function defined in terms of the Kronecker deltas as demonstrated in [2, 3]. The introduced functions and are chosen according to the number and the positions of the stochastic pendulums in the chain, which is derived from the linear coupling. For purely deterministic case according to which the concerned pendulum is deterministic as well as its neighbors (case a, section 2.1), 2 and 0. If the deterministic pendulum is adjacent to a stochastic one, and 1 or 1 if the stochastic pendulum is, respectively, the previous (case b) or the next (case c).

refers to the stochastic amplitude. When the concerned pendulum is stochastic (case d) and his neighbors are deterministic, replaces the deterministic amplitude in Eq. ( 9), 2 and 0. The number of generated complex equations depends on the number and positions of stochastic pendulums. In fact, the equations obtained in a deterministic case becomes in a stochastic one, where is the number of stochastic pendulums, the number of deterministic pendulums neighbors of the stochastic ones and 1 if the structure contains deterministic pendulums having deterministic neighbors and 0 otherwise.

To solve the obtained complex equations, a transformation of the amplitude to Cartesian form is needed [START_REF] Jallouli | Nonlinear dynamics of a 2D array of coupled pendulums under parametric excitation[END_REF] Substituting Eq. ( 10) into Eq. ( 9) and simplifying by , one can obtain two generic equations for the real and imaginary parts of each amplitude of the form

2 cos * sin n cos n sin ∑ , , Δ , (11)
2 cos * sin n ∑ cos n sin ∑ , , Δ , ∑ sin n ( 1 2 )
These obtained coupled algebraic equations are then numerically solved using the time integration Runge-Kutta method.

The stochastic analysis requires applying an uncertainty propagation method. In this paper, the statistical LHS method [START_REF] Helton | Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems[END_REF] is used. The choice of this method is due to its simple implementation consisting on a generation of a succession of deterministic computations according to a set of random variables and its reduced computing time with respect to the very time-consuming MC method. In fact, a great number of samples are typically required in MC to achieve good accuracy while partitioning the variability space into regions of equal probability and picking up one sampling point in each region allow reducing the number of samples needed for the LHS method implementation.

NUMERICAL SIMULATIONS

The uncertainty propagation using the LHS method through the generic discrete model proposed in this paper is numerically implemented in order to illustrate the effects of uncertainties on the behavior of a structure composed of three coupled pendulums (

3) for which, the design parameters are listed in Tab. 1. The length of the first pendulum is supposed to be uncertain and varies such as

1 ( 1 3 )
where and are respectively chosen dispersion level and Gaussian random variable.

The problem is numerically time-consuming when using the time integration Runge-Kutta method, in each LHS samples, to solve the obtained algebraic equations due to the necessity to sufficiently refine frequency step and vary initial conditions to obtain more stable solutions. This constraint imposes minimizing as possible the number of samples for the LHS method. 200 samples are therefore generated.

A statistical analysis of the frequency responses dispersion, both in generalized (modal) and physical coordinates, is carried out in order to investigate the robustness of the multimode branches against uncertainties in the multistability domain. Indeed, the results are presented at first in term of envelope of the displacement amplitude of each pendulum (Figs. 345). Then, a comparative study with respect to the deterministic solutions, presented by Fig. 2, permits to show the effect of uncertainties in the multistability domain through the quantification of the amplitude and frequency ranges variation (Tab. 2 and Fig. 7).

Deterministic frequency responses in modal and physical coordinates are shown in Fig. 2. In the multistability domain, different stable solutions exist and are generated by the modal interactions between the responses which are driven by the collective dynamics. Six stable solutions could thus be obtained in the multistability domain. The obtained stable branches can be classified as Single (SM), Double (DM) and Triple mode (TM) solutions. In fact, the only SM branch corresponds to the null trivial solution of the second amplitude and is presented by black curve in Fig. 2 (A). The DM branches, which are presented by red curves, correspond to the solutions for the first and third amplitudes, whereas the TM solutions (blue curves) result from the coupling between the three pendulums. A correspondence in term of bifurcation points between the amplitudes with respect to each branch type is observed. It is generated by the bifurcation topology transfer.

(A) (B) FIGURE 2. DETERMINISTIC RESPONSE AMPLITUDES OF THE PENDULUMS IN (A) GENERALIZED AND (B) PHYSICAL

COORDINATES.

In order to enable comparative study of the obtained results with respect to the deterministic ones, the samples of the LHS method are generated according to the initial conditions fixed for the deterministic case. The shown discontinuity in the envelope curves is due to the "lack" of initial conditions. However, it is important to note that varying more the initial conditions makes computation prohibitive and it still difficult to cover all possible solutions and so obtain smooth curves.

To illustrate the stochastic effect on the structure behavior when uncertainty is taken into account, two dispersion levels are considered: 2% and 5%. For 2%, the envelopes of the response amplitudes in modal coordinates are shown in Fig. 3, while Fig. 4 illustrates the amplitude envelopes for 5%. Moving to physical coordinates, response envelopes are presented by Fig. 5 for different dispersion levels.

It is clearly shown that increasing the dispersion level affects more the variability of the response amplitudes, both in modal and physical coordinates. Large variability is obtained for 5% and the overlap of envelopes makes detecting the extreme statistics of adjacent branches more difficult.

The effects of uncertainties on the dispersion of the response amplitudes depend on the position of each concerned pendulum with respect to the stochastic one. In fact, as shown in Figs. 34, the envelope width differs when moving from one pendulum to another. The greatest impact is on the response of the first pendulum since it contains the uncertain parameter. The second one is less affected while the behavior of the third pendulum remains deterministic.

When increasing the dispersion level, all branches are larger in terms of amplitude and frequency ranges. Consequently, the multimode solutions are enhanced and the multistability domain is wider.

In addition to the curves presented in Fig. 5 showing the response envelopes in physical coordinates, Tab. 2 and Fig. 7 illustrate the effects of uncertainties through calculating the intervals of variation of the amplitude (Δ ) and frequency (Δ ) ranges of the multistability domain, between extreme bifurcation points, according to the example showed in Fig. 6. As mentioned above and proved through Fig. 7, the first pendulum response is more affected by uncertainties than the second since its length is stochastic. The third pendulum has deterministic behavior because no coupling terms relate its equation with the one associated to the stochastic pendulum. 5 More important variation is detected for the frequency ranges (Δ ) than the amplitude ones (Δ ). When computing the ratio Δ Δ ⁄ , one can show the impact of uncertainties on the nonlinearity of the structure. Indeed, dispersion and nonlinearity are proportional. Consequently, the nonlinear aspect of the structure is enhanced.

σ -3 -2 -1 0 1 2 |A 1 |² 0 0.1 0.2 0.3 0.4 0.5 0.6 σ -3 -2 -1 0 1 2 |A 3 |² 0 0.05 0.1 0.15 0.2 σ -3 -2 -1 0 1 2 |A 2 |² 0 0.

CONCLUSION

Under the hypothesis of small imperfections, a generic discrete analytical model combining the multiple scales method, the perturbation theory and a standing-wave decomposition is proposed in order to investigate the collective dynamics of a nonlinear periodic structure composed of three linearly coupled pendulums under parametric uncertainties. The statistical LHS method permits to propagate uncertainty through the proposed model.

The analysis of the dispersion, in terms of variability of the frequency and amplitude intervals of the multistability domain according to different dispersion levels of uncertain input parameters, shows the effects of uncertainties on the stability and nonlinearity of the structure. The obtained results prove that the presence of uncertainties in a periodic structure enhances the multimode solutions introduced by the collective dynamics and interactions between pendulums by obtaining additional stable branches and widening the multistability domain. The uncertainties enhance also the nonlinear aspect of the structure. Further works concern the extension of the proposed method to large-size periodic structures.

FIGURE 1 .

 1 FIGURE 1. PERIODIC NONLINEAR COUPLED PENDULUMS STRUCTURE WITH PARAMETRICUNCERTAINTY.

FIGURE 3 .FIGURE 4 .FIGURE 5 .FIGURE 6 .

 3456 FIGURE 3. ENVELOPES OF RESPONSE AMPLITUDES OF THE (A) FIRST, (B) SECOND AND (C) THIRD PENDULUM INGENERALIZED COORDINATES FOR 2%.

  [1]. Boechler N., Daraio C., Narisetti R. K., Ruzzene M., Leamy M.J., 2010, "Analytical and experimental analysis of bandgaps in nonlinear one dimensional periodic structures", In Tsung-Tsong Wu and Chien-Ching Ma, editors, IUTAM Symposium on Recent Advances of Acoustic Waves in Solids, Vol 26 of IUTAM Bookseries, pp. 209-219. Springer Netherlands. [2]. Lifshitz R., Cross M. C., 2003, "Response of parametrically driven nonlinear coupled oscillators with application to micromechanical and nanomechanical resonator arrays", Physical Review B, 67, 134302. [3]. Bitar D., Kacem N., Bouhaddi N., Collet M., 2015, "Collective dynamics of periodic nonlinear oscillators under simultaneous parametric and external excitations", Nonlinear Dynamics, 82, pp. 749-766. [4]. Gutschmidt S., Gottlieb O., 2012, "Nonlinear dynamic behavior of a microbeam array subject to parametric actuation at low, medium and large dc-voltages", Nonlinear Dynamics, 67(1), pp. 1-36. [5]. Manktelow K., Leamy M. J., Ruzzene M., 2014, "Weakly nonlinear wave interactions in multi-degree of freedom periodic structures", Wave Motion 51, pp. 886-904. [6]. Andreassen E., Manktelow K., Ruzzene M., 2015, "Directional bending wave propagation in periodically perforated plates", Journal of Sound and Vibration, 335, pp. 187-203. [7]. Kissel G. J., 1988, "Localization in disordered periodic structures", Ph.D. thesis, Massachusetts Institute of Technology.

TABLE 1 .

 1 DESIGN PARAMETERS OF THE THREE COUPLED PENDULUMS STRUCTURE.

		.	. .		.
	0.25 0.062	9.10 -4	0.01	0.01	12.58

TABLE 2 .

 2 STOCHASTIC EVOLUTION OF THE FREQUENCY AND AMPLITUDE RANGES OF THE MULTISTABILITY DOMAIN WITH RESPECT TO DETERMINISTIC ONES ACCORDING TO THE DM AND TM.

				AMPLITUDE AND FREQUENCY
	RANGES OF THE MULTISTABILITY DOMAIN, ACCORDING TO
	THE DM AND TM, WITH RESPECT TO THE DETERMINISTIC
				RANGES.
			dof	DM = 2 %	= 5 %	= 2 %	TM	= 5 %
			1	172.86	272.86	165.06		287.95
	Δ %	2	110.47	113.95	143.18		165.91
			3	100	100	100		100
			1	104.76	114.15	113.29		126.44
	Δ %	2	100.94	103.10	100.83		97.39
			3	100	100	100		100
	Δ Δ	%	1 2 3	165.88 109.48 100	240.00 110.34 100	145.61 140.98 100		227.19 168.85 100


Diala Bitar

Najib Kacem

Noureddine Bouhaddi

Basins of attraction of coupled nonlinear resonators in periodic lattices

Collective dynamics in periodic lattices of coupled nonlinear Duffing-Van Der Pol oscillators is modeled and investigated under simultaneous external and parametric resonances. The resonators are coupled with linear and nonlinear springs. Numerical simulations have been performed in the case of two coupled oscillators for which the basins of attraction have been analyzed in the multistability domain in order to check the efficiency of the multimode branches.

Introduction

Interest in the nonlinear dynamics of periodic nonlinear lattices has grown rapidly over the last few years. Actually, it exists a practical need to understand nonlinearities and functionalize them in order to efficiently exploit the collective nonlinear dynamics of smart structures. For instance, Lifshitz et al. [START_REF] Lifshitz | Nonlinear dynamics of nanomechanical and micromechanical resonators[END_REF] investigated the dynamic behavior of an array of N coupled micro-beams using a discrete model. Manktelow et al. [START_REF] Manktelow | Multiple scales analysis of wave interactionns in a cubically nonlinear monoatomic chain[END_REF] focused on the interaction of wave propagation in a cubically nonlinear mono-atomic chain, while Bitar et al. [START_REF] Bitar | Energy transfer in externally driven periodic nonlinear structures[END_REF] investigated the bifurcation and enery transfers in periodic lattices of coupled nonlinear Duffing-Van Der Pol oscillators under an external excitation. Particularly, the basins of attraction can be used for qualitative as well as quantitative analysis of the collective dynamics robustness. In a nonlinear nanomechanical resonator, Kozinsky et al. [START_REF] Kozinsky | Basins of attraction of a nonlinear nanomechanical resonator[END_REF] experimentally probe the basins of attraction of two fixed points. Moreover, Sliwa et al. [START_REF] Śliwa | Periodic orbits, basins of attraction and chaotic beats in two coupled Kerr oscillators[END_REF] investigated the basins of attraction of two coupled Kerr oscillators. Furthermore, Ruzziconi et al. [START_REF] Ruzziconi | Multistability in an electrically actuated carbon nanotube a dynamical integrity perspective[END_REF] studied frequency response curves, behavior charts and attractor-basins phase portraits of a considered NEMS constituted by an electrically actuated carbon nanotube. In this context, an array of coupled Duffing-Van Der Pol oscillators under simultaneous primary and parametric resonances is developed and the distribution of the basins of attractions is analyzed for multistable solutions including single and multi-modes branches.

Model and basin of attraction analysis

The proposed model presents an array of a finite coupled Duffing-Van Der Pol oscillators, under simultaneous external and parametric excitations. The scaled equation of motion (EOM) governing the behavior of the n th resonator can be written as:

ün + ε ω 0 Q un + ω 2 0 u n + εh cos[2(ω 0 + εΩ)t]u n + 1 2 εd(-u n+1 + 2u n -u n-1 ) +α[(u n -u n+1 ) 3 + (u n -u n-1 ) 3 ] + ηu 2 n un = ε 3 2 g cos[(ω 0 + εΩ)t], (1)
where u n is the displacement of the n th oscillator, with fixed boundary conditions u 0 = u N +1 = 0, ω 0 and Ω are respectively the natural frequency and the detuning parameter. Q is the quality fractor, d represents the linear coupling, α is the cubic spring constant and η represents the Van Der Pol damping coefficient, h and g are parametric and external excitation amplitudes respectively, and ε is a small dimensionless parameter. The method of multiple time scales was used to solve the coupled EOM analytically. With the expectation that the motion of the resonator far from its equilibrium will be on the order of ε 1 2 , we try a solution of the form:

u n (t) = ε 1 2 N m=1 (A m (t) sin( nmπ N + 1 )e iω0t + c.c.) + ε 3 2 u (1) n (t) + • • • n = 1, ..., N, (2)
where T = εt is a slow time variable. The slowly varying amplitudes A m (T ) = (a m (T ) + ib m (T ))e iω0t obeys to 2N differential equations. In Figure 1, we display the response amplitude of the first oscillator, in function of the detuning parameter of two coupled oscillators, for specific design parameters. Remarkably, there are frequency bands where four stable solutions can exist. The multivaluedness of the response curves due to the nonlinearity has a significance from the physical point of view because it leads to jump phenomena wich are localized at the bifurcation points. The basins of attraction are numerically plotted to investigate the trajectories of the system response and the probability for which the system follows either the resonant or non-resonant, for single or double mode branches. Although, they are usually plotted in the phase plane (u n , un ), we chose to represent them in the Nyquist plane. Several numerical integrations of differential equations have been performed for a specific domain of initial conditions, in order to localize the maximum amplitude |A 1 | 2 in the steady-state domain. As shown in Figure 1, |A 1 | 2 takes one of four values, depending in the chosen initial conditions which allows for the representation of the corresponding basins of attraction.

Conclusion

The collective nonlinear dynamics of periodic nonlinear lattices was modeled for specific discrete systems of coupled Duffing-VDP oscillators under simultaneous primary and parametric excitations. The case of two coupled oscillators was investigated for a specific design parameters for which, the basins of attraction have been analyzed in the multistability domain for two coupled nonlinear oscillators to quantitatively asses the efficiency and reliabiliy of additional branches

Figure 1 :

 1 Figure 1: Response intensity as a function of the detuning parameter Ω for the first oscillator, where SMRB1 and SMRB2 are Single Mode Resonant Branches due respectively to primary and parametric resonances, SMNRB and DMNRB are respectively Single and Double Mode Non Resonant Branches.

Figure 2 :

 2 Figure 2: Variation of the basins of attraction in the Nyquist plane (a1, b1), with respect to the inital condition a2 for fixed detuning parameter Ω = 36 and b2 = 0. Magenta, purple, blue and cyan colors indicate respectively SMRB1, SMRB2, SMNRB and DMNRB.

Acknowledgments

This project has been performed in cooperation with the Labex ACTION program (contract ANR-11-LABX-01-01).