Loic Le Gratiet
email: loic.le-gratiet@edf.fr

Reformulation of the GCQC3 rule for three correlated excitation components with different spectra shapes

Keywords: response spectrum method, GCQC3, CQC3, directional correlation, direction combination, interaction hyper-ellipsoid

The response spectrum method is commonly used to design multi-degree-of-freedom structures subject to random vibrations. However, the common direction combination rules (SRSS, Newmark 30-100-100, Newmark 40-100-100) consider uncorrelated excitation components. Theoretically, it is always possible to rotate the directions in order to obtain principal directions along which the excitation components are not correlated. This problem is handle by the combinations CQC3 and GCQC3. However, these rules need simplifications to be applicable in an industrial framework. For specific applications, this simplifications are not consistent (e.g. for substructures anchored to a structure subject to plane crash). This article addresses the problem of reformulating the GCQC3 rule when vibrations along each direction are correlated, spectra have different shapes and without any simplifications. This new formulation avoids lot of calculations and makes the GCQC3 method industrially applicable in its most general form. In particular, it integrates the instantaneous correlations between excitation components. Furthermore, this formulation is propagated to the elliptical response envelopealso called interaction ellipsoidmethod.

Introduction

The response spectrum method (RSM) is commonly used to design multi-degree-of-freedom structures subject to random vibrations. This method works with linear structures and it uses a modal superposition procedure where each eigenmode corresponds to a single-degree-offreedom oscillator. Then, the peak responses of modal oscillators under each excitation component are evaluated and combined. The combinations give a probable upper bound for the temporal response of the structure [START_REF] Kiureghian | On response of structures to stationary excitation[END_REF][START_REF] Kiureghian | A response spectrum method for random vibrations[END_REF][START_REF] Erlicher | Seismic design by the response spectrum method: a new interpretation of elliptical response envelopes and a novel equivalent static method based on probable linear combinations of modes[END_REF]. Both modal and direction combinations have to be performed. A well-known rule for the modal combination is the Complete Quadratic Combination (CQC) [START_REF] Kiureghian | On response of structures to stationary excitation[END_REF][START_REF] Kiureghian | A response spectrum method for random vibrations[END_REF]. It allows for taking into account the correlation between the eigenmodes. For the direction combinations, the SRSS, the 30% [START_REF] Rosenblueth | Approximate design for multicomponent earthquakes[END_REF] and the 40% [START_REF] Newmark | Seismic design criteria for structures and facilities Trans-Alaska pipeline system[END_REF] rules are the most common. They all assume that the excitation components are uncorrelated. However, this assumption is frequently not verified in practical applications [START_REF] Penzien | Characteristic of 3-dimensional earthquake ground motions[END_REF][START_REF] Chen | Definition of statistically independent time histories[END_REF][START_REF] Hadjian | On the correlation of the components of strong ground motion[END_REF][START_REF] Marsan | Local amplification effects recorded by a local strong motion network during the 1997 Umbria-Marche Earthquake[END_REF]. In that case, the excitation components can be rotated to the principal orthogonal directions along which their correlations are null [START_REF] Penzien | Characteristic of 3-dimensional earthquake ground motions[END_REF][START_REF] Clough | Dynamics of Structures[END_REF][START_REF] Erlicher | Seismic design by the response spectrum method: a new interpretation of elliptical response envelopes and a novel equivalent static method based on probable linear combinations of modes[END_REF]. Then, a question of interest for seismic design is to establish the most critical orientation of the principal axes. Since it is often convenient to work with the structure axes, many works have been investigated the generalization of the CQC rule when the principal components are rotated with respect to the structure axes [1, 8, 10, 11 13, 14, 1]. In that framework, the GCQC3 rule has been developed [START_REF] Hernandez | Response to three-component seismic motion of arbitrary direction[END_REF][START_REF] Anastassiadis | Directions sismiques défavorables et combinaisons défavorables des efforts[END_REF]. However, this methods in its most general case could be difficult to handle for practical applications. For that reason, many works deal with this problem by suggesting simplifications. For example, [START_REF] Hernandez | Response to three-component seismic motion of arbitrary direction[END_REF] suggest a simplification by imposing restriction on the maximum inclination or intensity of a principal seismic component. Furthermore, most of them consider rotations only for horizontal directions [START_REF] Wilson | Three-dimensional dynamic analysis for multicomponent[END_REF][START_REF] Wilson | A clarification of the orthogonal effects in a three-dimensional seismic analysis[END_REF][START_REF] Smeby | Modal Combination rules for multicomponent earthquake excitation[END_REF][START_REF] Lopez | The critical angle of seismic incidence and the maximum structural response[END_REF][START_REF] Menun | A replacement for the 30%, 40%, and SRSS Rules for Multicomponent Seismic Analysis[END_REF][START_REF] Lopez | Critical response of structures to multicomponent earthquake excitation[END_REF]. In particular, in the work [START_REF] Wilson | Three-dimensional dynamic analysis for multicomponent[END_REF] the horizontal spectra are supposed to be identical and the solicitation components are uncorrelated. Then, a rule where the spectra have the same shape (i.e. they are equal up to a multiplicative constant) is presented in [START_REF] Wilson | A clarification of the orthogonal effects in a three-dimensional seismic analysis[END_REF]. The correlation are taken into account in [START_REF] Smeby | Modal Combination rules for multicomponent earthquake excitation[END_REF] where the spectra are supposed to have the same shape. A general rule with different spectra shapes and considering the correlations is given in [START_REF] Lopez | The critical angle of seismic incidence and the maximum structural response[END_REF]. Finally, the equations of [START_REF] Lopez | The critical angle of seismic incidence and the maximum structural response[END_REF] are reused in [START_REF] Menun | A replacement for the 30%, 40%, and SRSS Rules for Multicomponent Seismic Analysis[END_REF] with identical spectra shapes, the resulting formulation is called CQC3. This last work has been reformulating in [START_REF] Lopez | Critical response of structures to multicomponent earthquake excitation[END_REF] in order to avoid the evaluation of the critical orientation which implies lot of calculations.

The works cited previously essentially focus on structures under seismic solicitations. This article concerns the specific application of a substructure, such as equipment, which is anchored to a main structure subject to a vibrational solicitation. This solicitation could be seismic or induced by a plane crash. This application implies that the solicitations at the anchorages have different spectra shapes along the three directions and can present significant correlations between all directions (especially for vibrations induced by plane crash). Therefore, the suggested simplification of the GCQC3 rule [START_REF] Hernandez | Response to three-component seismic motion of arbitrary direction[END_REF][START_REF] Anastassiadis | Directions sismiques défavorables et combinaisons défavorables des efforts[END_REF] are not appropriate. Another particularity of this application is that the solicitations are determined with a time-history numerical model of the main structure. This point allows for reformulating the GCQC3 rule without performing simplifications by introducing the instantaneous correlations between the pairs of directions. Indeed, in practical application, it is straightforward to obtain these correlations by post-processing. The main result of this work is the establishment of an alternative formulation of the GCQC3 rule depending only on these instantaneous correlations. As a consequence, the resulting rule is quite easy to apply which is mandatory in the nuclear industry. For example, for sizing at plane crash, several thousand tests based on RSM have to be performed on substructures anchored to the main structures. Finally, the suggested reformulation allows to easily adapt the elliptical response envelopes formulae. The reader is referred to [START_REF] Erlicher | Seismic design by the response spectrum method: a new interpretation of elliptical response envelopes and a novel equivalent static method based on probable linear combinations of modes[END_REF][START_REF] Chu | Spectral treatment of actions of three earthquake components on structures[END_REF][START_REF] Gupta | Design of column sections subjected to three components of earthquake[END_REF][START_REF] Menun | Envelopes for seismic response vectors: I. theory[END_REF][START_REF] Menun | Envelopes for seismic response vectors: I I. Application[END_REF] for more detail about this approach. It allows to work with signed simultaneous responses for the design. More precisely, under the assumptions considered for the response spectrum method, the responses of interest belong to an interaction hyper-ellipsoid which has a closed-form expression. Then, by bounding the hyper-ellipsoid with a polytope, one can deduce concomitant responses which enfolds with high probability all possible responses. As for the Square Root of Sum of Squares (SRSS), the 40% and the 30% rules, the elliptical envelopes make the assumption of uncorrelated excitation components (by using a rotation of directions if necessary). This article generalizes this method by introducing the instantaneous correlations between the components.

3.

Response to three-dimensional excitations

Temporal response of a linear structures

The equations of relative motion for a linear multi-degree-of-freedom structure subject to a three-dimensional excitation 𝑺 ̈ can be written as:

𝑴𝑼 ̈+ 𝑪𝑼 ̇+ 𝑲𝑼 = -𝑴𝑶𝑨,
where 𝑴, 𝑪 and 𝑲 respectively are the mass, damping and stiffness matrices, 𝑼 is the vector of nodal displacements, 𝑨 = (𝐴 𝑥 (𝑡) 𝐴 𝑦 (𝑡) 𝐴 𝑧 (𝑡)) 𝑇 is the vector of accelerations exciting the structure and 𝑶 = ( 𝑶 𝒙 𝑶 𝒚 𝑶 𝒛 ) 𝑇 is the matrix coupling the degrees of freedom of the structure to the excitation motion. Let us introduce the matrix 𝝓 containing the eigenmodes of the eigenvalue equation (𝑲 -𝑴𝜔 2 )𝜙 = 0 and let us suppose that (𝝓 𝑻 𝑴𝝓) -1 𝝓 𝑻 𝑲𝝓 is diagonal. Then, using the transformation 𝑼 = 𝝓𝒒, the 𝑖th uncoupled eigenmode motion equation is:

𝑞 𝑖 ̈+ 2𝜉 𝑖 𝜔 𝑖 𝑞 𝑖 ̇+ 𝜔 𝑖 2 = -𝒑 𝒊 𝑨,
where 𝒑 𝒊 = 𝝓 𝑖 𝑇 𝑴𝑶/(𝝓 𝑖 𝑇 𝑴𝝓 𝑖 ) is the vector of participation factors for mode 𝑖, 𝝓 𝑖 is the 𝑖th column of 𝝓, 𝜔 𝑖 is its natural pulsation of mode 𝑖 and 𝜉 𝑖 is its damping ratio.

Consider a response of interest 𝑅(𝑡) that can be written as a linear combination of the nodal displacements, i.e. it has the form 𝑅(𝑡) = 𝒅 𝑇 𝑼. For a linear structure, responses can be forces, strains, stresses, nodal forces or loads. Then, using a mode superposition approach, 𝑅(𝑡) can be written as:

𝑅(𝑡) = ∑ ∑ 𝑆 𝑖,𝑘 (𝑡)𝑅 𝑖,𝑘 , 𝑛 𝑖 𝑘=𝑥,𝑦,𝑧
where 𝑛 is the number of eigenmodes, 𝑅 𝑖,𝑘 = (𝝓 𝑖 𝑇 𝑴𝑶 𝒌 /(𝝓 𝑖 𝑇 𝑴𝝓 𝑖 ))𝒅 𝑇 𝝓 𝑖 is the modal response weighted by its participation factor and 𝑆 𝑖,𝑘 (𝑡) is the response of an 1dof oscillator of pulsation 𝜔 𝑖 and damping ratio 𝜉 𝑖 to the excitation 𝐴 𝑘 (𝑡).

Response Power Spectral Density for correlated excitation components

Similarly to [START_REF] Kiureghian | On response of structures to stationary excitation[END_REF][START_REF] Kiureghian | A response spectrum method for random vibrations[END_REF] the excitation 𝑨 is supposed to be a zero-mean stationary random process. Nevertheless, the components of 𝑨 are not anymore considered as independent. Therefore, one cannot consider them separately as in [START_REF] Kiureghian | On response of structures to stationary excitation[END_REF][START_REF] Kiureghian | A response spectrum method for random vibrations[END_REF]. According to the multivariate Bochner's Theorem [START_REF] Gikhman | The theory of stochastic processes I[END_REF], the Power Spectral Density (PSD) of a multivariate stationary random process has the following form:

𝐺 𝑨 (𝜔) = ( 𝐺 𝐴 𝑥𝑥 (𝜔) 𝐺 𝐴 𝑥𝑦 (𝜔) 𝐺 𝐴 𝑥𝑧 (𝜔) 𝐺 𝐴 𝑥𝑦 (𝜔) 𝐺 𝐴 𝑦𝑦 (𝜔) 𝐺 𝐴 𝑦𝑧 (𝜔) 𝐺 𝐴 𝑥𝑧 (𝜔) 𝐺 𝐴 𝑦𝑧 (𝜔) 𝐺 𝐴 𝑧𝑧 (𝜔) ),
where ∀𝜔, 𝐺 𝐴 (𝜔) is nonnegative definite. In practice, this property needs to be verified a posteriori. Then, by the means of linearity, the power spectral density of the response is: where * stands for the complex conjugate and 𝐻 𝑖 (𝜔) is the complex frequency response function of a 1dof oscillator with damping ratio 𝜉 𝑖 and pulsation 𝜔 𝑖 :

𝐺 𝑅 (𝜔) = ∑ ∑ ∑ ∑
𝐻 𝑖 (𝜔) = 1 𝜔 𝑖 2 -𝜔 2 + 2𝑖𝜉 𝑖 𝜔 𝑖 𝜔 .
It can be emphasized that 𝐺 𝑅 (𝜔) is a natural generalization of the single excitation component case presented in [START_REF] Kiureghian | On response of structures to stationary excitation[END_REF][START_REF] Kiureghian | A response spectrum method for random vibrations[END_REF]. The 𝑚-order spectral moment of the response is: The term 𝜆 0,𝑖𝑗,𝑘𝑙 represents the covariance between 𝑆 𝑖,𝑘 (𝑡) and 𝑆 𝑗,𝑙 (𝑡) and 𝜆 2,𝑖𝑗,𝑘𝑙 represents the covariance between 𝑆 ̇𝑖,𝑘 (𝑡) and 𝑆 ̇𝑗,𝑙 (𝑡) [START_REF] Kiureghian | On response of structures to stationary excitation[END_REF][START_REF] Kiureghian | A response spectrum method for random vibrations[END_REF]. Therefore, the correlation between 𝑆 𝑖,𝑘 (𝑡) and 𝑆 𝑗,𝑙 (𝑡) writes:

𝜆 𝑚 = ∫ 𝜔 𝑚 𝐺 𝑅 (𝜔)𝑑𝜔 = ∑ ∑ ∑ ∑ 𝑅 𝑖,
𝜌 𝑖𝑗,𝑘𝑙 = 𝜆 0,𝑖𝑗,𝑘𝑙 √ 𝜆 0,𝑖𝑖,𝑘𝑘 𝜆 0,𝑗𝑗,𝑙𝑙 .

Then, the mean square of the response 𝑅(𝑡) can be written as follows: According to [START_REF] Kiureghian | A response spectrum method for random vibrations[END_REF] the mean of the peak response over a duration 𝜏 is given by:

𝜆 0 = ∑ ∑ ∑ ∑
𝐹 ̅ 𝑚𝑎𝑥 = 𝑝√𝜆 0 ,
where 𝑝 is a peak factor depending on the duration 𝜏, the mean zero-crossing rate and a shape factor for the response PSD [START_REF] Kiureghian | A response spectrum method for random vibrations[END_REF]. An estimate value of 𝑝 is given in Equation ( 16) in [START_REF] Kiureghian | A response spectrum method for random vibrations[END_REF]. Similarly to the response 𝑅(𝑡), the mean value of the maximum absolute response of a 1dof oscillator of pulsation 𝜔 𝑖 and damping ratio 𝜉 𝑖 to the excitation 𝐴 𝑘 (𝑡) writes:

𝑆 ̅ 𝑖,𝑘 = 𝐸 [max t |𝑆 𝑖,𝑘 (𝑡)|] = 𝑝 𝑖,𝑘 √𝜆 0,𝑖𝑖,𝑘𝑘 ,
where 𝑝 𝑖,𝑘 is given in Equation ( 16) in [START_REF] Kiureghian | A response spectrum method for random vibrations[END_REF] To obtain an estimate of 𝐹 ̅ 𝑚𝑎𝑥 , the expected maximum value of 𝑆 ̅ 𝑖,𝑘 is substituted by its observed maximal value max 𝑡 | 𝑆 𝑖,𝑘 (𝑡)|. Moreover, it is shown in [START_REF] Kiureghian | A response spectrum method for random vibrations[END_REF] that the ratio 𝑝 𝑝 𝑖,𝑘 ⁄ is around unity. Their argument is based on the slight dependency of the peak factor to the shape factor and on the fact that the response mean zero-crossing rate is a weighted root-mean-square of the mean zero-crossing rates 𝜈 𝑖,𝑘 . From these points, the mean of the peak response value can be approximated by:

𝐹 𝑚𝑎𝑥 = √ ∑ ∑ ∑ ∑ 𝐹 𝑖,𝑘 𝐹 𝑗,𝑙 𝜌 𝑖𝑗,𝑘𝑙 𝑛 𝑗=1 𝑛 𝑖=1 , 𝑙=𝑥,𝑦,𝑧 𝑘=𝑥,𝑦,𝑧
where 𝐹 𝑖,𝑘 is the peak modal response :

𝐹 𝑖,𝑘 = max 𝑡 | 𝑆 𝑖,𝑘 (𝑡)| × 𝑅 𝑖,𝑘 .
We emphasize that the approximation 𝑝 𝑝 𝑖,𝑘 ⁄ ≈ 1 can be avoided by using Equations ( 14)-( 18), ( 21) and ( 22) in [START_REF] Kiureghian | A response spectrum method for random vibrations[END_REF].

4.

Generalized CQC rule for correlated excitation components

Three-components stationary random excitation

A general closed-form expression is given in Section 3 for the response PSD of a structure enduring a multivariate stationary random excitation. By analogy with [START_REF] Kiureghian | On response of structures to stationary excitation[END_REF][START_REF] Kiureghian | A response spectrum method for random vibrations[END_REF], we suppose a particular probabilistic form for 𝑨. In particular, the following expression is considered for its PSD: ),

𝐺 𝑨 (𝜔) = 𝐺(𝜔
where 𝐺(𝜔) is the PSD a stationary random process, 𝜌 ̃𝑘,𝑙 is the instantaneous correlation between 𝐴 𝑘 (𝑡) and 𝐴 𝑙 (𝑡) and 𝛾 𝑘 2 is the second-order moment of 𝐴 𝑘 (𝑡). From a practical point of view, 𝛾 𝑘 2 is related to the magnitude of the excitation component 𝑘 and 𝜌 ̃𝑘,𝑙 has to be estimated statistically (instantaneous correlation means that no lag is introduced between the excitation components). The nonnegativity of 𝐺 𝑨 (𝜔) can be verified by evaluating the determinant of the right-hand term matrix. The response DSP is:

𝐺 𝑅 (𝜔) = ∑ ∑ 𝜌 ̃𝑘,

Generalized CQC rule (GCQC3)

In the expression of 𝐺 𝑨 (𝜔), we consider a unique PSD 𝐺(𝜔) for all margins. It means that the excitation components have the same probabilistic form. As in [START_REF] Kiureghian | On response of structures to stationary excitation[END_REF][START_REF] Kiureghian | A response spectrum method for random vibrations[END_REF], we consider a white noise excitation. It corresponds to the following DSP:

𝐺(𝜔) = 𝐺 0 ,
where 𝐺 0 is a constant. Since the magnitude of the DSP is already taken into account through the coefficient 𝛾 𝑘 , we set 𝐺 0 = 1. In this framework, the covariance between 𝑆 𝑖,𝑘 (𝑡) and 𝑆 𝑗,𝑙 (𝑡) has the following expression [START_REF] Kiureghian | On response of structures to stationary excitation[END_REF][START_REF] Kiureghian | A response spectrum method for random vibrations[END_REF]: ). The expression of 𝐹 𝑚𝑎𝑥 provides a generalization of the CQC ruledenoted GCQC3 -when excitation components are correlated.

𝜆 0,

Discussions about the GCQC3 rule

As presented in Subsection 4.2, one can obtain the peak response value of a quantity of interest by considering the following combination rule:

√ ∑ ∑ ∑ 𝜌 𝑖,𝑗 𝜌 ̃𝑘,𝑙 𝐹 𝑖,𝑘 𝐹 𝑗,𝑙 𝑛 𝑖,𝑗=1
.

𝑙=𝑥,𝑦,𝑧 𝑘=𝑥,𝑦,𝑧 Contrary to the SRSS, the Newmark 30-100-100 and the Newmark 40-100-100 rules, it requires to perform the modal and the direction combinations simultaneously. However, the peak modal responses (𝐹 𝑖,𝑘 ) 𝑘=𝑥,𝑦,𝑧 𝑖=1,…,𝑛 are those already evaluated when using classical modal and direction combinations. Therefore, this generalized rule only requires in addition the evaluation of the correlation coefficients 𝜌 ̃𝑘,𝑙 . Moreover, the independent case corresponds to 𝜌 ̃𝑘,𝑙 = 𝛿 𝑘=𝑙 where 𝛿 is the delta of Kroenecker. This assumption leads to the following rule: Then, using the triangle inequality the following inequality holds: 𝐹 𝑚𝑎𝑥 ≤ √𝐹 𝑥,𝑥 + 𝐹 𝑦,𝑦 + 𝐹 𝑧,𝑧 + 2|𝜌 ̃𝑥,𝑦 |√𝐹 𝑥,𝑥 √𝐹 𝑦,𝑦 + 2|𝜌 ̃𝑦,𝑧 |√𝐹 𝑦,𝑦 √𝐹 𝑧,𝑧 + 2|𝜌 ̃𝑥,𝑧 |√𝐹 𝑥,𝑥 √𝐹 𝑧,𝑧 .

√ ∑ ∑ ∑ 𝜌 𝑖,
As a consequence, when the instantaneous correlations are unknown, an upper bound for the response can be provided.

5.

GCQC3 method for elliptical response envelopes

As for the SRSS, Newmark 30-100-100 and Newmark 40-100-100 rules, the GCQC3 is valid for a single response [START_REF] Nie | On the correct application of the 100-40-40 rule for combining responses due to three directions of earthquake loading[END_REF]. As state in [START_REF] Nie | On the correct application of the 100-40-40 rule for combining responses due to three directions of earthquake loading[END_REF], to deal with multiple response, SRSS, 30%, 40% and GCQC3 rules can be performed for each response and then all combinations of maximum responses have to be considered. This can lead to strong overestimates. The elliptical response envelopes method can be used to avoid them. It allows to work with concomitant responses by considering their cross-correlation [START_REF] Erlicher | Seismic design by the response spectrum method: a new interpretation of elliptical response envelopes and a novel equivalent static method based on probable linear combinations of modes[END_REF]. Like the classical rules, it considers uncorrelated excitation components by arguing that there are principal directions leading to un-correlation. However, it could be difficult in practice to work with rotated directions. To address this issue, the elliptical response envelopes is re-written by integrating excitation component correlations.

First of all, let us consider a multiple-response 𝑹(𝑡) = (𝑅 1 (𝑡) ⋯ 𝑅 𝑝 (𝑡)) 𝑇 with 𝑝 components and 𝑅 𝑚 (𝑡) = 𝒅 𝒎 𝑇 𝑼. It can be written [START_REF] Erlicher | Seismic design by the response spectrum method: a new interpretation of elliptical response envelopes and a novel equivalent static method based on probable linear combinations of modes[END_REF]: For 𝑝 = 1, this is equivalent to the equation of 𝑅(𝑡) in Subsection 3.1. From the GCQC3 method, we have the flowing inequality: As demonstrated in [START_REF] Erlicher | Seismic design by the response spectrum method: a new interpretation of elliptical response envelopes and a novel equivalent static method based on probable linear combinations of modes[END_REF], 𝜶(𝒕) must fulfills the following inequality:

𝑹(𝑡) = ∑ ∑ 𝛼 𝑖,𝑘 (
|𝑅 𝑚 (𝑡)| = |𝑭 𝑚 𝑇 𝜶(𝑡)| ≤ 𝐹 𝑚 𝑚𝑎𝑥 = √ ∑ ∑ ∑ 𝜌 𝑖,
𝜶(𝑡) 𝑇 𝑯 -1 𝜶(𝑡) ≤ 1,
which means that 𝜶(𝒕) belong to an hyper-ellipsoid of dimension 3𝑛. Therefore, 𝛼(𝑡) verifies:

{ 𝑹(𝑡) = 𝑭 𝑇 𝜶(𝑡) 𝜶(𝑡) 𝑇 𝑯 -1 𝜶(𝑡) ≤ 1 .

To obtain a unique value for 𝜶(𝒕), the one minimizing the quadratic form 𝜶(𝒕) 𝑻 𝑯 -1 𝜶(𝒕) is considered [START_REF] Erlicher | Seismic design by the response spectrum method: a new interpretation of elliptical response envelopes and a novel equivalent static method based on probable linear combinations of modes[END_REF]. It leads to the following closed-form expression:

𝜶 𝑚𝑖𝑛 (𝑡) = 𝑯𝑭(𝑭 𝑇 𝑯𝑭) -1 𝑹(𝑡).

Integrating the value of 𝜶 𝑚𝑖𝑛 (𝑡) into 𝜶(𝑡) 𝑇 𝑯 -1 𝜶(𝑡) ≤ 1, we obtain the inequality:

𝑹 𝑇 (𝑡)(𝑭 𝑇 𝑯𝑭) -1 𝑹(𝑡) ≤ 1,
which implies that the multiple response 𝑹(𝑡) belongs to an hyper-ellipsoid of dimension 𝑝. It is so-called the elliptical response envelopes. From it, we can deduce a set of concomitant responses. In practice, they are deduced from the vertices of a polytope bounding the elliptical response envelopes.

GCQC3 with rigid mode

To take into account the missing mass due to the projection 𝑼 = 𝝓𝒒, an additional mode called "rigid mode" or "pseudo-mode" is often considered. It assumes that the dynamic amplification of the modes with frequency larger than the one of the Zero-Period Acceleration 𝜉 𝑍𝑃𝐴 is negligible. Then, the response becomes:

Example

Let us consider an embedded structure composed of three beams and four nodes. The structure and its eigenmodes are illustrated in Figure 1. The node N3 in Figure 1 has a punctual mass of 10kg and the masses of the others nodes and beams are null. The structure is embedded at node N0, the beam lengths equal 1m and they have square sections with length 3cm and thickness 5mm. The Young's modulus of the beams is 200GPa and their Poisson coefficient is 0.3. Furthermore, the degrees of freedom along the Y axis are blocked as well as the rotations around the X and Z axes. The structure has two modes presented in Figure 1 that take into account the whole mass. The frequencies, the participation factors and the modal masses are presented in Table 1.

For the participation factors, the eigenmodes are normalized such that the largest translation value equals one. Moreover, the modal dumping is the same for the two modes and equals 𝜉 = 5%. The excitations along the X and Z directions have the following form:

𝐴(𝑡) = Γ(𝑡)𝑌(𝑡),
where Γ(𝑡) is the Gamma modulation function defined by:

Γ(𝑡) = 𝛼 𝑡 𝛽-1 exp(𝛾𝑡),
and 𝑌(𝑡) is a zero-mean stationary Gaussian process with an exponential covariance function:

𝑘(ℎ) = σ 2 exp (- |ℎ| 𝜃 ).
The coefficients of the excitation model are given in Table 2.

Table 2. Model coefficients for the excitation.

8.

Conclusion

This article presents a reformulation of the GCQC3 rule when vibrations along each direction are correlated, spectra have different shapes and without any simplifications. Such a rule is necessary in the nuclear industry, for example to size substructures subject to vibrations induced by a plane crash. Furthermore, the reformulation is also required since the general GCQC3 rule without simplifications implies lot of calculations which are difficult to achieve in an industrial context. The obtained GCQC3 rule only requires the evaluation of the instantaneous correlations between excitation components and the peak modal responses. Therefore, its implementation is as easy as the one of the classic rules (SRSS, 40%, 30%). Furthermore, this rule can be adapted for the elliptical response envelopes method which allows to consider concomitant responses. Finally, the efficiency of the GCQC3 rule has been illustrated in a numerical example. It emphasizes that it allows to fit the maximum responses of a time-history analysis even in presence of correlations. Furthermore, a good consistency between the time-history responses and the elliptical response envelopes is also observed.

9.

  𝑘 𝑅 𝑗,𝑙 𝜆 𝑚,𝑖𝑗,𝜆 𝑚,𝑖𝑗,𝑘𝑙 = 𝑅𝑒 [∫ 𝜔 𝑚 𝐺 𝐴 𝑘𝑙 (𝜔) ∞ 0 𝐻 𝑖 (𝜔)𝐻 𝑗 * (𝜔)𝑑𝜔].

  𝑗 𝛿 𝑘=𝑙 𝐹 𝑖,𝑘 𝐹 𝑗,𝑙 𝑛 𝑖,𝑗=1 𝑙=𝑥,𝑦,𝑧 𝑘=𝑥,𝑦,𝑧 = √ ∑ ∑ 𝜌 𝑖,𝑗 𝐹 𝑖,𝑘 𝐹 𝑗,𝑙 𝑛 𝑖,𝑗=1 𝑘=𝑥,𝑦,𝑧 , which is exactly the classical SRSS rule. Introducing the notation: 𝐹 𝑘,𝑙 = ∑ 𝜌 𝑖,𝑗 𝐹 𝑖,𝑘 𝐹 𝑗,𝑙 𝑛 𝑖,𝑗=1 , we can explicit the GCQC3: (𝐹 𝑚𝑎𝑥 ) 2 = ∑ ∑ 𝜌 ̃𝑘,𝑙 𝐹 𝑘,𝑙 𝑙=𝑥,𝑦,𝑧 𝑘=𝑥,𝑦,𝑧 = ∑ 𝐹 𝑘,𝑘 𝑘=𝑥,𝑦,𝑧 + 2𝜌 ̃𝑥,𝑦 𝐹 𝑥,𝑦 + 2𝜌 ̃𝑦,𝑧 𝐹 𝑦,𝑧 + 2𝜌 ̃𝑥,𝑧 𝐹 𝑥,𝑧 , where 𝐹 𝑘,𝑘 is the square of the classical direction response deduced from a CQC rule. Finally, if only the classical CQC responses are available, one can deduce an upper bound to 𝐹 𝑚𝑎𝑥 . Indeed, the Cauchy-Schwarz inequality implies: ∑ ∑ |𝜌 ̃𝑘,𝑙 𝐹 𝑘,𝑙 | 𝑙=𝑥,𝑦,𝑧 𝑘=𝑥,𝑦,𝑧 ≤ ∑ ∑ |𝜌 ̃𝑘,𝑙 |√ ∑ 𝜌 𝑖,𝑗 𝐹 𝑖,𝑘 𝐹 𝑗,𝑘 𝑛 𝑖,𝑗=1 𝑙=𝑥,𝑦,𝑧 𝑘=𝑥,𝑦,𝑧 √ ∑ 𝜌 𝑖,𝑗 𝐹 𝑖,𝑙 𝐹 𝑗,𝑙 𝑛 𝑖,𝑗=1.

  𝑡)𝑭 𝑖,𝑘 𝑛 𝑖 , 𝑘=𝑥,𝑦,𝑧 where -1 ≤ 𝛼 𝑖,𝑘 (𝑡) ≤ 1, 𝑭 𝑖,𝑘 = ( 𝐹 1,𝑖,𝑘 = max 𝑡 | 𝑆 𝑖,𝑘 (𝑡)| × 𝑅 1,𝑖,𝑘 ⋮ 𝐹 𝑝,𝑖,𝑘 = max 𝑡 | 𝑆 𝑖,𝑘 (𝑡)| × 𝑅 𝑝,𝑖,𝑘 ), and 𝑅 𝑚,𝑖,𝑘 = (𝝓 𝑖 𝑇 𝑴𝑶 𝒌 /(𝝓 𝑖 𝑇 𝑴𝝓 𝑖 ))𝒅 𝒎 𝑇 𝝓 𝑖 . The multiple-response can be written in the following form: 𝑹(𝑡) = 𝑭 𝑇 𝜶(𝑡), with 𝑭 = (𝑭 1 ⋯ 𝑭 𝑝 ), 𝑭 𝑚 = ([𝐹 𝑚,𝑖,𝑥 ] 𝑇 𝑖=1,…,𝑛 [𝐹 𝑚,𝑖,𝑦 ] 𝑇 𝑖=1,…,𝑛 [𝐹 𝑚,𝑖,𝑧 ] 𝑇 𝑖=1,…,𝑛 ) 𝑇 , and 𝜶(𝑡) = ([𝛼 𝑖,𝑥 (𝑡)] 𝑇 𝑖=1,…,𝑛 [𝛼 𝑖,𝑦 (𝑡)] 𝑇 𝑖=1,…,𝑛 [𝛼 𝑖,𝑧 (𝑡)] 𝑇 𝑖=1,…,𝑛 ) 𝑇 .

  𝑅(𝑡) = ∑ 𝐴 𝑘 (𝑡)𝑅 𝑛+1,𝑘 + ∑ 𝑆 𝑖,𝑘 (𝑡)𝑅 𝑖,𝑘 𝑛 𝑖 𝑘=𝑥,𝑦,𝑧 = ∑ 𝑆 𝑖,𝑘 (𝑡)𝑅 𝑖,𝑘 𝑛+1 𝑖 , where 𝑅 𝑛+1,𝑘 = 𝑲 -1 𝑴𝝓 𝑅,𝑘 , 𝝓 𝑅,𝑘 = 𝑶 𝒌 -∑ 𝑅 𝑖,𝑘 𝑛 𝑖=1is the rigid mode and 𝑆 𝑛+1,𝑘 (𝑡) = 𝐴 𝑘 (𝑡) is the ground acceleration for direction 𝑘. Then, considering the following correlations: 𝑐𝑜𝑟 (𝑆 𝑛+1,𝑘 (𝑡), 𝑆 𝑗,𝑘 (𝑡)) = 𝛿 𝑗=𝑛+1 , 𝑐𝑜𝑟 (𝑆 𝑛+1,𝑘 (𝑡), 𝑆 𝑛+1,𝑙 (𝑡)) = 𝜌 ̃𝑘,𝑙 , and using the notation 𝜌 𝑛+1,𝑗 = 𝛿 𝑗=𝑛+1 , taking into account the rigid mode is straightforward : 𝐹 𝑚𝑎𝑥 = √ ∑ ∑ ∑ 𝜌 𝑖,𝑗 𝜌 ̃𝑘,𝑙 𝐹 𝑖,𝑘 𝐹 𝑗,

Figure 1 .

 1 Figure 1. Eigenmodes of the structure.

Table 1 .

 1 Modal frequencies, participation factors and masses.

𝛼 1/ 6 Figure 3 .

 63 Figure 3. Maximum displacements with respect to the excitation component correlations.Finally, we present in the series of figures in Figure4, the results of the elliptical response envelopes method with the GCQC3 rule for displacements along the X and Z directions. The black points in Figure4represent the values of the time-history simulations. For each correlation value, the solid lines represent the elliptical envelopes for the 200 excitations and the dashed line is the mean of these envelopes. Figure4highlights the good behavior of the GCQC3 rule since the elliptical envelopes fit the cloud of concomitant responses obtained with the time-history analysis.

Figure 4 .

 4 Figure 4. GCQC3 elliptical response envelopes for different component correlations.

  𝐺 𝐴 𝑘𝑙 (𝜔)𝑅 𝑖,𝑘 𝑅 𝑗,𝑙

			𝑛	𝑛
				𝐻 𝑖 (𝜔)𝐻 𝑗 * (𝜔),
	𝑘=𝑥,𝑦,𝑧	𝑙=𝑥,𝑦,𝑧	𝑖=1	𝑗=1

  𝑅 𝑖,𝑘 𝑅 𝑗,𝑙 𝜌 𝑖𝑗,𝑘𝑙 √𝜆 0,𝑖𝑖,𝑘𝑘 𝜆 0,𝑗𝑗,𝑙𝑙

			𝑛	𝑛
				,
	𝑘=𝑥,𝑦,𝑧	𝑙=𝑥,𝑦,𝑧	𝑖=1	𝑗=1

  with the mean zero-crossing rate 𝜈 𝑖,𝑘 = √

				𝜆 2,𝑖𝑖,𝑘𝑘
				𝜆 0,𝑖𝑖,𝑘𝑘
	𝐹 ̅ 𝑚𝑎𝑥 = √ ∑ ∑ ∑ ∑ 𝑛 𝑛 𝑗=1 𝑖=1 𝑙=𝑥,𝑦,𝑧 𝑘=𝑥,𝑦,𝑧	𝑝 2 𝑝 𝑖,𝑘 𝑝 𝑗,𝑙	𝑅 𝑖,𝑘 𝑅 𝑗,𝑙 𝜌 𝑖𝑗,𝑘𝑙 𝑆 ̅ 𝑖,𝑘 𝑆 ̅ 𝑗,𝑙	.

/𝜋 and the response PSD shape factor 𝛿 𝑖,𝑘 = √1 -

𝜆 1,𝑖𝑖,𝑘𝑘 2 𝜆 2,𝑖𝑖,𝑘𝑘 𝜆 0,𝑖𝑖,𝑘𝑘

. Using the expression of the expected maximum value 𝑆 ̅ 𝑖,𝑘 , the mean of the peak response can be expressed as follows:

  𝑙 𝛾 𝑘 𝛾 𝑙 ∑ ∑ 𝐺(𝜔)𝑅 𝑖,𝑘 𝑅 𝑗,𝑙 the covariance between 𝑆 𝑖,𝑘 (𝑡) and 𝑆 𝑗,𝑙 (𝑡) writes: 𝜆 0,𝑖𝑗,𝑘𝑙 = 𝛾 𝑘 𝛾 𝑙 𝜌 ̃𝑘,𝑙 𝑅𝑒 [∫ 𝐺(𝜔)

			𝑛	𝑛
				𝐻 𝑖 (𝜔)𝐻 𝑗 * (𝜔),
	𝑘=𝑥,𝑦,𝑧	𝑙=𝑥,𝑦,𝑧	𝑖=1	𝑗=1
			∞ 0	𝐻 𝑖 (𝜔)𝐻 𝑗

and * (𝜔)𝑑𝜔].

  𝜉 𝑗 𝜉 𝑖 𝜔 𝑖 𝜔 𝑗 (𝜉 𝑖 𝜔 𝑖 + 𝜉 𝑗 𝜔 𝑗 )𝜔 𝑖 𝜔 𝑗

		𝑖𝑗,𝑘𝑙 = 𝛾 𝑘 𝛾 𝑙 𝜌 ̃𝑘,𝑙	2𝜋 2 + 𝜔 𝑗 2 + 4𝜉 𝑖 𝜉 𝑗 (𝜔 𝑖 2 ) 2 -𝜔 𝑗 (𝜔 𝑖 2 ) + 4(𝜉 𝑖 2 + 𝜉 𝑗 2 )𝜔 𝑖 2 𝜔 𝑗 2	(𝜉 𝑖 𝜔 𝑖 + 𝜉 𝑗 𝜔 𝑗 ),
	and their correlation is:
		𝜌 𝑖𝑗,𝑘𝑙 = 𝜌 ̃𝑘,𝑙	8 √ (𝜔 𝑖 2 -𝜔 𝑗 2 ) 2 + 4𝜉 𝑖 𝜉 𝑗 (𝜔 𝑖 2 + 𝜔 𝑗 2 ) + 4(𝜉 𝑖 2 + 𝜉 𝑗 2 )𝜔 𝑖 2 𝜔 𝑗 2	= 𝜌 ̃𝑘,𝑙 𝜌 𝑖𝑗 .
	According to Subsection 3.2, the mean of the peak response value is:
				𝑛
				𝐹 𝑚𝑎𝑥 = √ ∑ ∑ ∑ 𝜌 𝑖,𝑗 𝜌 ̃𝑘,𝑙 𝐹 𝑖,𝑘 𝐹 𝑗,𝑙	,
				𝑘=𝑥,𝑦,𝑧	𝑙=𝑥,𝑦,𝑧	𝑖,𝑗=1
	For the single component case, we retrieve the classical CQC given in [4,5] (i.e. 𝐹 𝑚𝑎𝑥 =
	√ ∑ 𝑛 𝑖,𝑗=1	𝜌 𝑖,𝑗 𝐹 𝑖,𝑘 𝐹 𝑗,𝑘

  𝑗 𝜌 ̃𝑘,𝑙 𝐹 𝑚,𝑖,𝑘 𝐹 𝑚,𝑗,𝑙

					𝑛
						= √𝑭 𝑚 𝑇 𝑯𝑭 𝑚 ,
		𝑘=𝑥,𝑦,𝑧	𝑙=𝑥,𝑦,𝑧	𝑖,𝑗=1
	Where:
		[𝜌 𝑖,𝑗 ] 𝑖,𝑗=1,…,𝑛	[𝜌 ̃𝑥,𝑦 𝜌 𝑖,𝑗 ] 𝑖,𝑗=1,…,𝑛 [𝜌 ̃𝑥,𝑧 𝜌 𝑖,𝑗 ] 𝑖,𝑗=1,…,𝑛
	𝑯 =	[𝜌 ̃𝑥,𝑦 𝜌 𝑖,𝑗 ] 𝑖,𝑗=1,…,𝑛	[𝜌 𝑖,𝑗 ] 𝑖,𝑗=1,…,𝑛	[𝜌 ̃𝑦,𝑧 𝜌 𝑖,𝑗 ] 𝑖,𝑗=1,…,𝑛	.
	(	[𝜌 ̃𝑥,𝑧 𝜌 𝑖,𝑗 ] 𝑖,𝑗=1,…,𝑛 [𝜌 ̃𝑦,𝑧 𝜌 𝑖,𝑗 ] 𝑖,𝑗=1,…,𝑛	[𝜌 𝑖,𝑗 ] 𝑖,𝑗=1,…,𝑛 )

The correlation between the excitation components is driven by the correlation between the Gaussian processes 𝑌(𝑡) for each direction X and Z. Examples of excitations for instantaneous correlations of 0.8 and -0.6 are presented in Figure 2. Then, the displacements obtained with the response spectrum analysis are compared to the ones obtained with a time-history simulation. To perform the comparison, 200 couples of excitations along the X and Z directions have been generated for each correlation 𝜌 ̃𝑥,𝑧 between -1 and 1. The mean of the maximum displacements with respect to the solicitation correlation is presented in Figure 3. The results in Figure 3 show that the response spectrum method with the GCQC3 rule fits very well the time-history simulation. This emphasizes the efficiency of this rule for taking into account the correlations between the excitation components. Furthermore, Figure 3 confirms that the classical SRSS rule behaves very badly. Indeed, its relative error can be greater than 50%. It is worth to highlight that the SRSS and the GCQC3 rules are equivalent when the correlation coefficient is zero.


G Jouan

F Collin
email: f.collin@ulg.ac.be

P Kotronis

BENDING TESTS ON REINFORCED CONCRETE BEAMS: NUMERICAL MODELLING USING A SECOND GRADIENT THEORY

Keywords: Second Gradient Theory, Crack, Reinforced Concrete, Damage

 and has this far been used mainly to reproduce the behavior of soils. It is here applied to two bending tests of reinforced concrete beams.

INTRODUCTION

Strain localization in quasi-brittle materials, or more generally in materials exhibiting strain softening is a well-known problem [START_REF] Hill | Bifurcation phenomena in the plane tension test[END_REF] [9] [START_REF] Rudnicki | Conditions for the localization of deformation in pressure-sensitive dilatant materials[END_REF]. It is a phenomenon that can be clearly observed in experimental tests, and yet cannot be modeled with classical continuum mechanics models. Analytically, the solution is a crack with zero energy dissipation [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF]; numerically, it leads to a pathological mesh dependency. These shortcomings are due to the lack of an internal length in the continuum model.

Several ways of introducing an internal length have since been proposed. Non local integral models were first formulated by Pijaudier-Cabot et al. [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF] and then extended to damage theory. The gradient models of Aifantis [START_REF] Aifantis | On the microstructural origin of certain inelastic models[END_REF] enrich the conventional plasticity and damage theories with gradient of the internal variables. This type of model can be shown to be equivalent to the integral type models. More recently, strain localization due to damage has been treated using the thick level set approach [START_REF] Moës | A level set based model for damage growth: The thick level set approach[END_REF].

A rather natural way of introducing (indirectly) a length parameter in a continuum model is to somehow account for the microstructure of the material. The general class of so called microstructured models or higher order continuum models allows for the description of the kinematics of the microstructure by using an additional tensor in the displacement field. Higher order continuum theories can be traced back to the works of the Cosserat brothers [START_REF] Cosserat | Théorie des corps déformables[END_REF] and have been generalized and properly formulated by Germain [START_REF] Germain | La méthode des puissances virtuelles en mécanique des milieux continues: Première partie : théorie du second gradient[END_REF][START_REF] Germain | The method of virtual power in continuum mechanics. Part 2[END_REF] using the method of virtual power.

The local second gradient model developed by Chambon et al. can be seen as a particular case of a higher order continuum and has been used to regularize problems involving localization in geomaterials [START_REF] Chambon | Plastic continuum with microstructure, local second gradient theories for geomaterials: localization studies[END_REF][START_REF] Matsushima | Second gradient models as a particular case of microstructured models: a large strain finite elements analysis[END_REF]. It is used here to model two different bending tests on reinforced concrete beams. The first one is a 3 point bending test with a beam of 5m span, the second one is a 4 point bending test with a beam of 6.1m span.

THE SECOND GRADIENT MODEL AND CONSTITUTIVE LAWS

Theoretical framework

As detailed in the seminal work of Germain, using the virtual power method one can chose a field of virtual displacement to describe the proper kinematics of the continuum (including its microstructure). The internal stresses, limit conditions and equilibrium equations then appear naturally as long as the linear form representing the virtual power is correctly defined and that it respects the principle of material independence.

Following this, if one chooses the virtual displacement field as the "field of continuous and continuously differentiable velocities", the principle of virtual work yields:

* + Σ , * Ω - * Ω - * + * = 0 (1)
Where u i is the virtual displacement field, D ij is the symmetric part of its gradient, σ ij is the macro stress (conjugate of first gradient) and Σ ijk is a double stress (conjugate of the second gradient of the velocity field). G i is the classical body force, p i is a contact force and P i a double contact force. Du i refers to the normal derivative: Du i = n k u i / x k . Here and henceforth * denotes a virtual quantity.

The equilibrium equations and boundary conditions are given by:

- Σ + = 0 (2) - Σ - Σ - Σ + Σ - Σ = (3)
Σ = (4) Where Dq/Dx j is the tangential derivative.

Numerical implementation

The second gradient of the displacement in the weak formulation of the problem necessitates the use of C 1 elements in a finite element code. This can be avoided by introducing a new field v ij imposed to be equal to the gradient of u i by lagrange multipliers [START_REF] Matsushima | Large strain finite element analysis of a local second gradient model: application to localization[END_REF]. The new weak formulation of the problem then becomes: * + Σ , * Ω -! " * - * # Ω -$ = 0 (5)

! * " -# Ω = 0 ( 6
)
The problem is then discretized using 9 nodes finite elements, where 8 of the nodes are used for u i , 4 for v ij and the center node for λ ij . This element has been implemented in the finite element code Lagamine and the problem is solved using the classical Newton-Raphson method.

Constitutive laws

Here we suppose that the model uses two constitutive equations, one linking the macro stress to the first gradient of the displacement and the other linking the double stress to the second gradient of the displacement. The two equations are also supposed to be decoupled. The first gradient law can be any classical constitutive law (involving or not damage, plasticity etc…). The constitutive laws can be written in rate form as:

% = & % $ (7)
' % = ( )* % ,)*

If the second gradient law is an isotropic elastic one, the general form of the constitute law involves five material parameters. particular form of this latter constitutive equation has been proposed by Chambon et al. [START_REF] Matsushima | Second gradient models as a particular case of microstructured models: a large strain finite elements analysis[END_REF] with a single modulus B.

+ , , , , , , , , - ' % ... ' % .. ' % . . ' % . ' % .. ' % . ' % . ' % / 0 0 0 0 0 0 0 0 1 + , , , , , , , - 2 0 0 2/2 0 0 2/2 0 0 2/2 0 0 2/2 0 0 2 0 2/2 2/2 0 2/2 0 2/2 0 0 2/2 0 2/ 0 2/2 2/2 0 2/2 0 0 2/2 2/2 0 0 2/2 0 2/2 2/2 0 2 0 0 2/2 0 2/2 0 2/2 0 0 2/2 0
The internal length does not appear clearly from these equations. However, analytical solution of localization in a one problem with a bilinear first gradient law of moduli A 1 and A 2 , exhibits a localized strain band, whose length l s is given by [START_REF] Chambon | One-dimensional localization studied with a second grade model[END_REF]: 56 7 89 : . 2 ;< = > ?@ 9 : . : 56 A9 : 2

It is shown that the internal length mainly a function of the ratio of B

The use of the second gradient model allows for the proper representation of strain localization. It does not however imply uniqueness of the solution. For a one dimensional tension test or a biaxial compression test, the number of bands and their positions can vary [START_REF] Bésuelle | Switching deformation modes in post-localization solutions with a quasibrittle material[END_REF]. This is in accordance with experimental results which tend to be poorly reproducible when strain localization.

The first gradient law used here for t concrete behavior is the classical damage law [START_REF] Mazars | Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture du béton de structure[END_REF]. The equivalent strain is defined as:

B $C D〈B〉 G : 〈B〉 G
Where 〈 〉 G denote the positive part and is the double dot product. The damage is separated in two part D t and D c , one due to the tension, and the other due to compression.

If the second gradient law is an isotropic elastic one, the general form of the constitute five material parameters. A particular form of this latter constitutive equation has been proposed by Chambon et al.

2 0 2/2 2/2 /2 0 0 0 0 2 / 0 0 0 0 0 0 0 1 + , , , , , , , - % .
The internal length does not appear clearly from these equations. However, analytical solution of localization in a one-dimensional problem with a bilinear first gradient law of , exhibits a localized strain y [14]:

= > W (10)
It is shown that the internal length l s is B over A 2 . The use of the second gradient model allows for the proper representation of strain however imply uniqueness of the solution. For a onedimensional tension test or a biaxial compression test, the number of bands and ]. This is in accordance with experimental results which tend to be poorly reproducible when involving

The first gradient law used here for the behavior is the classical Mazars The equivalent strain is [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF] denote the positive part and : he damage is , one due to the tension, and the other due to compression.

NR XY Z (12)
R XY Z [START_REF] Moës | A level set based model for damage growth: The thick level set approach[END_REF] Where ε d0 , A c , A t , B parameters. The total damage is given

[ U \ U ;1 [
[ U is equal to 0 when there is no traction and to 1 when there is no compression. Its value is given by:

[ U ∑^〈B 〉 G B B $C
B U are the principal strains caused by the positive part of the principal stresses

In the tests modeled here, using Mazars damage law in a two dimensional setting, equation ( 10) is not valid, but gives nevertheless an indication of the width of the localization zones at their inception

THREE POINT BENDING

Experimental test

This experimental test was conducted on a reinforced concrete beam according to the specifications of the CEOS.FR benchmark [START_REF]Réponse au benchmark statique monotone du projet national CEOS[END_REF]. The beam is of thickness b=200mm, height h=500mm and 5000mm span (see: figure 2). It was subjected to cyclic controlled load.

Numerical model

The bending test is modeled as a two dimensional problem using the 9 node finite element described above. Two meshes have , B c and A t are material The total damage is given by:

[ U ? \ I (14)
when there is no traction when there is no compression. Its 〉 B U _ [START_REF] Bésuelle | Switching deformation modes in post-localization solutions with a quasibrittle material[END_REF] the principal strains caused by the of the principal stresses. In the tests modeled here, using Mazars damage law in a two dimensional setting,

) is not valid, but gives nevertheless an indication of the width of the ation zones at their inception.

THREE POINT BENDING TEST

experimental test was conducted on a reinforced concrete beam according to the CEOS.FR benchmark The beam is of thickness b=200mm, h=500mm and 5000mm span (see: figure 2). It was subjected to cyclic controlled 3 point bending test: beam dimensions [START_REF]Réponse au benchmark statique monotone du projet national CEOS[END_REF] The bending test is modeled as a two dimensional problem using the 9 node finite element described above. Two meshes have been used for the simulations. The first mesh consists of 5180 elements, 4148 of which are second gradient elements and 1032 truss elements representing the horizontal reinforcement. The average size of the concrete elements for this mesh is of 0.02mx0.035m. The second mesh consists of 13494 elements with an average size of 0.01mx0.017m for the concrete elements. Concrete and steel elements are supposed to be perfectly bonded. The end nodes at each lower extremities of the beam are blocked vertically; the right node is also blocked horizontally. On the upper part where the loading is applied, and near the supports at both ends of the beam, an elastic linear law is introduced to prevent from unwanted damage. An elastoplastic law with isotropic hardening is used for the reinforcement. The elastic modulus B for the second gradient constitutive law is equal to 1.5.MN.

Results

Overall, the numerical results are in good accordance with the experimental data. Figure 5 shows the numerical global force versus displacement at the center of the beam compared to the experimental results (note: for the experimental test the beam was loaded and unloaded cyclically whereas in the simulation the beam was loaded with a monotonic increasing displacement). The force displacement graph exhibits the classical reinforced concrete behavior in three stages: In the first stage, concrete and steel stay both in the elastic regime; then concrete starts to damage and the slope on the force displacement curve changes. Finally, steel enters in its plastic phase and the second change in the slope appears. The final slope of the numerical simulation doesn't match completely the experimental one. This is due to the second gradient constitutive law which considers no damage (elastic).

Figure 5 shows the pattern of the damage variable in concrete at different stages of loading and for the two mesh sizes. The damage pattern develops with sudden "peaks" which experimentally correspond to developing cracks. The crack opening is not modeled directly in this simulation as the displacement field remains continuous, but it can be calculated from the damage model either through energetic equivalence or by simply measuring the jump in displacement between two points located on the opposite sides of a damaged zone. This obviously works only when the damaged bands are clearly separated. The width and separation of the damage bands can be controlled by changing the internal length, which in our case would mean changing the slope of either or both the first gradient and second gradient constitutive laws.

The damage distributions for the two different meshes do not match exactly although the global force displacement curve The damage pattern develops with sudden "peaks" which experimentally correspond to developing cracks. The crack opening is not modeled directly in this simulation as the displacement field remains continuous, but it can be calculated from the damage model either through energetic equivalence or by simply measuring the jump in displacement between two points located on the opposite sides of a damaged zone. This obviously rks only when the damaged bands are clearly separated. The width and separation of the damage bands can be controlled by changing the internal length, which in our case would mean changing the slope of either or both the first gradient and second gradient

The damage distributions for the two different meshes do not match exactly although the global force displacement curves are almost identical. This is due to the fact that, as stated above, the second gradient method doesn't restore solution and small differences different solutions. All the solutions are still physically acceptable.

For similar bending tests, non damage models which define an equivalent strain by averaging over a certain distance have a tendency to develop damage on the upper compressed part of the beam the local strain is not high enough to caus compressive damage. This is due to the averaging over an area. There is no such problem with this model as all the variables are local.

FOUR POINT BENDING T 4.1 Experimental test

In this case, the experimental point bending test on a large reinforced concrete beam. (6.10m x 0.8m x 1.6m) ( figure 7, « Mock up During the experimental test, the beam was cast and then let to shrink freely. was then fixed to a testing bench wit of pre-stressed steel bars with two rows of pistons force of 2250kN and then unloaded completely. . This is due to the fact as stated above, the second gradient tore the unicity of the differences can trigger different solutions. All the solutions are still For similar bending tests, non-local damage models which define an equivalent strain by averaging over a certain distance l c , have a tendency to develop damage on the compressed part of the beam, even when the local strain is not high enough to cause compressive damage. This is due to the averaging over an area. There is no such problem with this model as all the variables

FOUR POINT BENDING TEST

In this case, the experimental test is a 4 point bending test on a large reinforced 10m x 0.8m x 1.6m) (see Mock up », www.ceosfr.org). During the experimental test, the beam was cast and then let to shrink freely. The beam was then fixed to a testing bench with the help stressed steel bars. The beam is loaded with two rows of pistons up to an imposed force of 2250kN and then unloaded The concrete and steel characteristics are given in the following tables:

Numerical model

For this test, it was chosen to model only one half of the beam and to consider a symmetry condition to reduce the computational cost. We used once again the two dimensional 9 nodes second gradient elements to model concrete and one dimensional truss elements for the reinforcement. Two different meshes are used: a coarser one with 5200 elements (4640 second gradient elements and 560 truss elements) and another one with 10000 elements (9209 second gradient elements). The material parameters adopted here don't match exactly the experimental ones. This is because in our calculations only the mechanical loading phase is modeled (the free shrinkage phase that damaged concrete is not simulated). The elastic modulus for the second gradient law B is equal to 1.5.MN.

Results

The force-vertical displacement curve for the numerical model and the experimental data is given in figure 9. The results for the two different meshes are the same and are in good accordance with the experimental ones (note: only the loading phase is reproduced). During the experiment the beam was unloaded; this is the reason why the global force displacement curve shows only two phases compared to the three phases of the previous test. As was the case with the three point bending test, the damage distributions for the two different meshes do not match exactly, although the global force displacement curves are almost identical.

CONCLUSIONS

A second gradient model has been used to model two bending tests of reinforced concrete beams. The results show that the model is able to reproduce the force-displacement curves obtained experimentally. Damage localizes into bands whose width is controlled by the model parameters. The uniqueness of the solution is however not restored. For the two experiments modeled here, the bands are of the same size but their numbers vary (figure 5 and 10). Being a local theory, the second gradient method avoids the limitations caused by the use of a non-local definition of the equivalent strain.

These results are encouraging and represent the first steps toward a wider use of the local second gradient method for concrete structures.

Figure 1 :

 1 Figure 1: 2 nd gradient finite element[START_REF] Matsushima | Second gradient models as a particular case of microstructured models: a large strain finite elements analysis[END_REF]

  ;1 : U ? B $C : U M NO V QR ST NR G. Jouan, F. Collin and P. Kotronis 3

Figure 2 :

 2 Figure 2: 3 point bending test:

Figure 3 :

 3 Figure 3: 3 point bending test: finite element meshes

Figure 4 :

 4 Figure 4: 3 point bending test: force-displacement

Figure 5 Figure 5 :

 55 Figure 5.1. 2mm displacement

Figure 6 :

 6 Figure 6: 4 point bending test:

Figure 7 :

 7 Figure 7: 4 point bending test: reinforcement [19]

Figure 8 :

 8 Figure 8: 4 point bending test: numerical model

Figure 9 :

 9 Figure 9: 4 point bending test: force-displacement

Figure 6

 6 Figure 6 presents the pattern of the damage variable in concrete at different stages of loading and for the two mesh sizes.

Figure 10 :

 10 Figure 10.1. F = 940kN

Table 1 :

 1 3 point bending test: concrete material parameters

	E (GPa)	ε d0	A t	B t	A c	B c	β
	37.2	9.10 -5 0.7	6800 0.42 780	1.1

Table 2 :

 2 3 point bending test: reinforcement material parameters

		Upper trusses
	E (GPa)	σ (MPa)	Area (cm²)
	195	466	16.085 (2HA32)
		Lower trusses
	E (GPa)	σ (MPa)	Area (cm²)
	195	466	1.0053(2HA8)

Table 3 :

 3 4 point bending test: concrete characteristics

	Category E (GPa)	ν	f c (MPa) f t (MPa)
	C50/60	40.2	0.19	63.7	4.65

Table 4 :

 4

		4 point bending test: reinforcement
		characteristics
	Catégorie	E (GPa)	f e (MPa)	ν
	Fe500	200	500	0.3

Table 5 :

 5 4 point bending test: concrete material parameters

	E (GPa)	ε d0	A t	B t	A c	B c	β
	25.2	1.1.10 -4 0.7 6800 0.42 780 1.1

point bending test: experimental test[START_REF] Livrable | Prise en compte des résultats expérimentaux dans les modélisations[END_REF] are indicated in the

ACKNOWLEDGEMENTS

The authors would like to thank the ANR MEFISTO (Maîtrise durable de la fissuration des infrastructures en béton, projet ANR-2008-Sustainable Cities Program, project number: VD08_323065, 2011) and the Programme National CEOS.fr (www.ceosfr.org) for their financial support.

MEFISTO (Maîtrise durable de la fissuration des infrastructures en béton). Projet ANR-2008-Sustainable Cities

Program, project number: VD08_323065, 2011.


Marcela Carvalho

Maxime Ferrera

Alexandre Boulch

Julien Moras

Bertrand Le Saux

Pauline Trouvé-Peloux

Co-learning of geometry and semantics for online 3D mapping

scientifiques 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.

Introduction

Autonomous navigation is conditioned by the ability of sensing and analyzing the environment to take new decisions. In this context, accurate 3D reconstruction and semantic understanding of the scenes are critical. Indeed, building a 3-Dimensional (3D) map of the scene including semantic information allows to plan future trajectories accordingly to the tasks to perform.

Over the past years, improvements on data acquisition techniques and processing made possible reconstructing 3D scenes in multiple ways. On the one hand, active sensors are now mature technology and some variants gain special attention, like LIDARs, which produce dense and reliable point clouds [START_REF] Cole | Using laser range data for 3D SLAM in outdoor environments[END_REF]; and RGB-D sensors that generates corresponding depth maps which can be combined to scene reconstruction [START_REF] Whelan | Elasticfusion: Real-time dense slam and light source estimation[END_REF]. On the other hand, passive approaches like Structurefrom-Motion (SfM) are also commonly adopted to recover 3D relations between points and objects from a set of 2D images. In SfM, we distinguish offline methods from online methods. Offline variants, also denoted as photogrammetry, usually exhaustively process all data before global reconstruction. In opposition, online approaches handle information incrementally to perform reconstruction while data are being acquired, as in Simultaneous Localization and Mapping (SLAM) and in Visual Odometry (VO). However, these latter techniques usually rely on geometric features and not on semantic information, though it is an important feature to perform more specialized and complex navigation tasks.

In this work, we present a new approach to jointly learn geometry and semantics for 3D mapping. The proposed pipeline consists of two steps, corresponding to different levels of data aggregation (Fig. 1). First, at image level, a multi-task network estimates a depth map and a semantic segmentation map. Then, these geometric and semantic features are accumulated into a global representation where the semantic mesh of the scene is extracted from the 3D representation, which allows scene understanding and planning of further actions.

In details, the main contributions of this paper are the following. The first key point is the joint use of geometric and machine learning approaches. As illustrated in Figure [START_REF] Tylecek | The second workshop on 3d reconstruction meets semantics: Challenge results discussion[END_REF], a raw depth map is estimated from a pair of images using stereo and then is refined through a convolutional neural network. A second key point is the co-learning of depth and semantic segmentation from the raw depth map and an RGB images. Hence the proposed network performs multiple tasks at once, with mutual benefit. We show that this approach leads to better performances than independent predictions of depth and semantic segmentation. Finally, on the contrary to global, offline reconstruction methods, our approach is incremental and hence is conceptually compatible with autonomous navigation and robotics.

The paper is organized as follows: section 2 presents works related to the problem, section 3 describes our semantic reconstruction pipeline and finally section 4 evaluates our method with quantitative and qualitative results on the 3D Reconstruction Meets Semantics 2018 (3DRMS) Challenge dataset, which contains series of stereo sequences generated over a simulated garden.

Related work

Perception for autonomous navigation has been a great topic of interest in the last two decades. As cameras became cheap and easy to embed while still offering rich information, vision-based SLAM methods grew more and more popular [4,[START_REF] Forster | SVO: Semidirect Visual Odometry for Monocular and Multicamera Systems[END_REF]. SLAM allows a robot to localize itself with respect to the environment. Either this environment is unknown, and its 3D structure is simultaneously estimated, or the environment is already known, and a previously built map can be used [START_REF] Lynen | Get out of my lab: Large-scale, real-time visual-inertial localization[END_REF][START_REF] Schneider | maplab: An open framework for research in visual-inertial mapping and localization[END_REF]. In the latter case, such maps can be obtained by regular SLAM methods, i.e.. building the map of the environment and then using it for self localization. Maps can also be built offline by SfM algorithms such as Colmap [START_REF] Schönberger | Structure-from-motion revisited[END_REF] or OpenMVG [START_REF] Moulon | Others: Openmvg. an open multiple view geometry library[END_REF] before being used for real-time localization. All these approaches for offline or online map construction, take only the geometric structure of the scene into account. However, a few works proposed to also benefit from semantic information, yielding in semantic SLAM [START_REF] Civera | Towards semantic slam using a monocular camera[END_REF]. Indeed, this allows to get better maps and increase the localization reliability [START_REF] Mccormac | Semanticfusion: Dense 3d semantic mapping with convolutional neural networks[END_REF][START_REF] Schönberger | Semantic visual localization[END_REF]. Using RGB-D data, a pipeline using random forests for creating semantic maps in 2D and 3D was proposed in [START_REF] Hermans | Dense 3d semantic mapping of indoor scenes from rgb-d images[END_REF]. More recently, [START_REF] Ma | Multi-view deep learning for consistent semantic mapping with rgb-d cameras[END_REF] applied joint learning with neural networks over multiple RGB-D views to generate better 2D semantic maps, but did not reconstruct corresponding 3D models. With respect to all these approaches, ours offers a functional pipeline from 2D images to 3D reconstruction with semantics. With respect to the latter ones, semantics and geometry have a better integration directly in the network.

The joint use of geometry and semantics has been investigated in the previous edition of the 3DRMS challenge [START_REF] Strisciuglio | Trimbot2020: an outdoor robot for automatic gardening[END_REF]. The dual objective was 3D geometry reconstruction and semantic classification. The proposed baseline links Colmap [START_REF] Schönberger | Structure-from-motion revisited[END_REF] and SegNet [START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for image segmentation[END_REF] (for 2D classification, then projected in 3D). Both entries [START_REF] Guerry | SnapNet-R: Consistent 3d multi-view semantic labeling for robotics[END_REF][START_REF] Taguchi | Semantic 3D reconstruction using depth and label fusion[END_REF] used semantics during the reconstruction to filter out outliers, but with worst performances than the baseline. With respect to these approaches, we propose to learn the fusion of semantics and geometry through a multi-task network.

Fully Convolutional Networks (FCNs) [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF][START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for image segmentation[END_REF][START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] have been widely used for many tasks in computer vision. In brief, they are dense prediction methods which intend to assign information back onto the original pixels positions. Semantic segmentation is a common domain of application for such dense prediction networks.We focus here on the approaches which benefit from geometric information. FuseNet [START_REF] Hazirbas | Fusenet: Incorporating depth into semantic segmentation via fusion-based cnn architecture[END_REF] uses two interlaced encoders and a single decoder for semantic segmentation from RGB-D data. Alternatively, in [START_REF] Audebert | Semantic Segmentation of Earth Observation Data Using Multimodal and Multi-scale Deep Networks[END_REF], the authors introduce residual fusion using a small network to merge the outputs of two SegNets applied to different sensor modalities. A finer (though more complex) approach, 3D graph neural network [START_REF] Qi | 3d graph neural networks for rgbd semantic segmentation[END_REF], consists in considering information extracted from the local 3D graph of adjacency and using it in the segmentation network. [START_REF] Guerry | SnapNet-R: Consistent 3d multi-view semantic labeling for robotics[END_REF] proposed 3D-consistent data augmentation to incorporate the geometry directly in the training set. Among all these approaches, the one which has most in common with ours is FuseNet [START_REF] Hazirbas | Fusenet: Incorporating depth into semantic segmentation via fusion-based cnn architecture[END_REF], since they share solving the fusion problem by a highly-integrated network. However, our network goes beyond simple fusion, and address a multi-task problem, with semantic segmentation and depth adjustment.

FCNs have also been applied to other tasks such as monocular depth prediction [START_REF] Eigen | Predicting Depth, Surface Normals and Semantic Labels with a Common Multi-Scale Convolutional Architecture[END_REF][START_REF] Xu | Multi-scale continuous crfs as sequential deep networks for monocular depth estimation[END_REF][START_REF] Kendall | What uncertainties do we need in bayesian deep learning for computer vision[END_REF][START_REF] Carvalho | On regression losses for deep depth estimation[END_REF]. These techniques exploit spatial correlation based on structured information (e.g., linear perspective, textures) to produce reliable depth maps from 2D scenes. To cite but a few, we then focus on approaches with open-source code. Based on SegNet, Laina et al. [START_REF] Laina | Deeper depth prediction with fully convolutional residual networks[END_REF] exploit residual connections [START_REF] He | Deep residual learning for image recognition[END_REF] and fast up-projection blocks. In D3-Net [START_REF] Carvalho | On regression losses for deep depth estimation[END_REF], the network consists on a densely connected encoder [START_REF] Huang | Laurens: Densely connected convolutional networks[END_REF] and a U-Net like decoder structure to predict refined depth estimation. With respect to these approaches, our method also uses depth from geometry as an input, and refines the 2D depth map using semantic constraints, which yields in better depths than with stereo or monocular prediction.

Proposed approach

As presented in figure 1, our method is composed of two computation levels: depth and semantic maps generation; 3D data accumulation for surface reconstruction. These tasks are combined sequentially and result in an accurate method for 3D scene reconstruction. From beginning to end, we use a stereo sequence to produce a semantic mesh.

Our main idea is to learn jointly the depth and the semantic segmentation in a multi-task deep neural network framework. Besides, we also benefit from geometric depth estimation methods. Indeed, raw depth map estimated from a pair of stereo images with geometric approach are used as inputs of the multitask network. In the following, we describe in details the four sub-tasks of Fig. 1.

Depth estimation

The first step of the proposed 3D reconstruction pipeline consists in estimating depth maps from stereo views. In brief, the calibration of stereo cameras allows estimating the relative pose of the right camera with respect to the left one, as well as their distortion parameters. Using these informations, the left and right images may be undistorted and rectified in order to be aligned. Once aligned, the depth of corresponding points in both images can be estimated from the known baseline between the cameras, their focal length and the disparity between the two points.

Two different stereo matching algorithms have been tested to compute disparity maps: ELAS [START_REF] Geiger | Efficient large-scale stereo matching[END_REF] and SGBM [START_REF] Hirschmueller | Stereo processing by semiglobal matching and mutual information[END_REF]. ELAS is a probabilistic method based on the triangulation of robust matches to create reliable support points. These support points then serve as priors for the disparity search and allows the computation of a Maximum-A-Posteriori (MAP) estimate over the remaining pixels. ELAS was tested using the LIBELAS implementation, without post-processing. Indeed, some post-processing is commonly applied in order to get better and more dense disparity maps. However, the neural networks might be more prone to handle raw disparity maps than human-interpretable ones. Otherwise, SGBM is a semi-global method which estimates disparity by minimizing an energy function made of the Sum of Absolute Distances (SAD) over a local window and a smoothness term. SGBM was tested using its OpenCV implementation and no post-processing were applied.

Their respective accuracy results on the training sequences 0001 and 0224 of the 3DRMS dataset are displayed in table 1. SGBM gives slightly better results than ELAS, so the depth maps produced by SGBM are the ones we choose to use in the proposed pipeline.

Semantic segmentation and depth enhancement

The task at hand here is the reconstruction of a semantic mesh of the given scene. Hence, the objective is twofold: reconstruct the geometry of the scene (3D localization of the mesh vertices) and identify semantics (attach a label to each mesh element).

Though, the geometric depth estimation previously described is not sufficient for global surface reconstruction. Indeed, the small baseline of the stereo sensor might lead to huge estimation errors. Hence, a refinement step is needed to produce better depth maps. As the geometric errors mostly occur on edges, using the RGB image as an additional information would lead the network to produce sharper edges. Besides, as shown in [START_REF] Hazirbas | Fusenet: Incorporating depth into semantic segmentation via fusion-based cnn architecture[END_REF][START_REF] Qi | 3d graph neural networks for rgbd semantic segmentation[END_REF], semantic segmentation benefits from both RGB images and depth maps.

These considerations motivate the proposed approach of a multi-task fully convolutional neural network for a joint prediction of depth and semantic segmentation. The proposed architecture was inspired in FuseNet [START_REF] Hazirbas | Fusenet: Incorporating depth into semantic segmentation via fusion-based cnn architecture[END_REF] and is presented on figure 2. The Multi-task Network has an encoder-decoder structure, with two branches for the encoder, and 2 independent branches for the decoder (one for semantics, one for depth estimation). Contrarily to the original implementation of FuseNet, we add skip connections between the encoder and decoder parts to improve spatial information flow over the network. Branches in the contractive part take the RGB and raw depth inputs respectively and as feature maps are generated, they are melt from the depth branch to the RGB input branch. Also, depth refinement is performed in a residual manner, adding the correction to the input raw map.

Filtering

Even though depth is enhanced using the multi-task network, a few errors remain when an object occludes another. In this case, the network tends to smooth the transition between objects and overlook small details (such as tree leafs for example). To avoid unwanted outliers in later stages of the 3D reconstruction, we apply the following filtering operations. First, points labeled as sky are removed. Second, points from uncertain object borders are identified and removed. These borders correspond to transitions between objects at different depths, so we compute the gradient of depth over the image and remove all pixels for which the gradient norm is greater than a fixed threshold (empirically set to 0.05). Finally, we also test in the experiments an additional filter: erosion, or removal of the neighbors of a point considered as an outlier. As we will show, this produces more precise but less complete reconstructions.

Iterative 3D map construction

The 3D reconstruction module is based on Truncated Signed Distance Function (TSDF) modeling. This technique estimates a scalar field which represents the approximate distance of every points in the 3D space to the nearest surface. In practice, the field is estimated over a 3D discretization of the world and only close to the surfaces. The distance estimates are signed: positive outside of the object and negative for the inside. Hence, the zero crossing is an implicit representation of the surfaces of the objects present in the scene and a Marching Cube algorithm is used to recovers the mesh. The TSDF implementation used in this paper is based on OpenChisel [START_REF] Klingensmith | Chisel: Real time large scale 3d reconstruction onboard a mobile device[END_REF]. In order to estimate the distance field, the 3D space is discretized into voxels and the filtered depth maps are integrated into the TSDF according to the poses of the camera. The depth maps are first clipped in order to only process 3D points within a clipping range distance from the camera (in practice from 0.5 to 5-10 meters). In addition to distance estimations, we also add semantic classification fusion. Thus the module can take as a new input, either the label image resulting from classification or directly the classification scores (cf. Section 3.2). These semantic inputs are processed in the same way as the depth maps, that is the voxels integrate the semantic scores in addition to the distance-to-nearest-surface values. When all the frames have been integrated, a filter removes the voxels which do not contain accurate enough distance values. For each remaining voxel, the semantic label is selected as the one with highest score. In practice, the voxel grid resolution is set to 3cm. The mesh is finally generated by applying Marching Cube over this voxel grid.

Table 1. Error measurements adopted to evaluate semantic segmentation (left) and depth estimation performances (right). P corresponds to the predictions and GT the ground truth. C is the number of classes. Variables di and di are the ground truth and prediction respectively, and N is the total number of pixels.

Metric Definition

Overall Acc. (OA)

P ∩GT |GT | Average Acc. (A.Acc.) 1 C C i=1 P i ∩GT i |GT i ] Average IoU (A. IoU) 1 C C i=1 P i ∩GT i P i ∪GT i

Metric Definition

Abs. error

1 N N i=0 |d i -di | d i RMS 1 N N i=0 (di -di) 2

Experiments results

In this section, each step of the semantic reconstruction pipeline is evaluated on the 3DRMS dataset. The data consists in four training sequences with ground truth and a test sequence for challenge evaluation (for which the ground truth remains undisclosed to participants). We further divide the training set in train and validation to present evaluation scores and comparable visual results. Precisely, we created two folds from the training data: fold 1 with training scenes from sequences 0128, 0160, 0224 and testing scenes from sequence 0001; and fold 2 with training scenes from sequences 0001, 0128, 0160 and testing scenes from sequence 0224.

In the following, semantic segmentation (Section 4.1), depth estimation (Section 4.2) and global 3D reconstruction (Section 4.3) are evaluated on this dataset and compared to state of the art approaches. Metrics used for comparison are presented in tables 1.

Semantic segmentation

Semantic segmentation from 2D images is one of the tasks of the 3DRMS Challenge. Our architecture is evaluated against the two state-of-the-art approaches: U-Net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] and FuseNet [START_REF] Hazirbas | Fusenet: Incorporating depth into semantic segmentation via fusion-based cnn architecture[END_REF]. Performances are computed according to the metrics presented in Table 1-left. Results are presented in Table 2. It shows a clear improvement of the performance of semantic segmentation when using the proposed multi-task network. The best configuration is multi-task training using the raw depth map from SGBM. In Fig. 3, examples of semantic segmentations of some 2D images from the dataset are displayed. It shows that co-learning enforces consistency with respect to the 3D structure. Indeed, neighbor pixels with the same depth (i.e. also close to each other in 3D) tend to get the same semantic label.

Depth estimation

The quality of 3D reconstruction highly depends on the estimation of an accurate depth map. In this paper, we propose to generate a precise depth map by refining a raw one obtained with a stereo pair. This process is one of the tasks of our multi-task network. In the following, we compare the performances of depth estimation using traditional stereo methods such as ELAS [START_REF] Geiger | Efficient large-scale stereo matching[END_REF] and SGBM [START_REF] Hirschmueller | Stereo processing by semiglobal matching and mutual information[END_REF] with the performances of the refined depth estimates. We also evaluate the performance of a state-of-the-art single-image depth estimation approach, referred to as D3-Net [START_REF] Carvalho | On regression losses for deep depth estimation[END_REF].

The various depth map predictions are first compared using standard error measurements previously used for the same purpose [START_REF] Li | Depth and surface normal estimation from monocular images using regression on deep features and hierarchical crfs[END_REF][START_REF] Eigen | Predicting Depth, Surface Normals and Semantic Labels with a Common Multi-Scale Convolutional Architecture[END_REF] and defined in Table 1(right). We also provide the proportion of points with a deviation less than a given value in Fig. 4(a) and the RMS function with respect to the ground truth distance in Fig. 4(b). Several conclusions can be drawn from this. First, refinement of the geometric depth map using a multi-task neural network highly improves the depth estimation accuracy. Indeed for geometric approaches, only 40 % of the points have a deviation lower than 2m, while it reaches 80 % using the proposed multi-task approaches. One can note that improvement is specifically significant for small depth range, between 0 to 5m, which is crucial for safe autonomous navigation. Second, all multi-task learning approaches show similar results, with slightly better performances when the raw input comes from SGBM. Furthermore, our tests also show that using a state-of-the-art FCN for single-image depth estimation outperforms the purely geometric approaches according to these standard, global metrics. As discussed in the following, this result can be explained by a better depth map segmentation obtained by deep learning approach.

Figure 5 shows examples of depth maps obtained with the various geometric or multi-task approaches. A geometric method such as SGBM results in accurate depth estimates but with a low quality segmentation of the depth map. On the contrary, a deep learning approach such as D3-Net shows an excellent depth segmentation, but produces biased depth values. Finally, the proposed approach which benefits from both geometrical and deep learning techniques shows the best results both in terms of accuracy and quality of depth segmentation.

Reconstruction

The reconstruction is performed with OpenChisel [START_REF] Klingensmith | Chisel: Real time large scale 3d reconstruction onboard a mobile device[END_REF]. In this section, we present the final results of the reconstruction for both test sets 001 and 224 of the 3DRMS dataset. We evaluate the geometric quality of the reconstruction according to the depth map filtering strategy. We also provide quantitative and qualitative results on the semantics in 3D.

Geometric reconstruction As defined in [START_REF] Strisciuglio | Trimbot2020: an outdoor robot for automatic gardening[END_REF], the quality of the reconstruction can be evaluated from two points of view. First, each point of the ground truth must be close to a point of the reconstructed scene, this is the completeness of the reconstruction, i.e.. it express how well the whole scene has been discovered and reconstructed. Second, each point of the reconstruction must be close to a point of the ground truth, this is accuracy. The accuracy aims at evaluating how well the reconstruction fits to the ground truth. In practice a good reconstruction is a compromise between completeness and accuracy; filling the space with points would improve the completeness while selecting only few points, well positioned, would improve the accuracy. We use the following metrics for quantitative results (definition in table 4):

from ground truth to reconstruction: the average distance of GT point to the mesh, and the completeness (the distance d such that 90% of the GT points are at distance less than d to the reconstruction). from reconstruction to GT : the average distance of mesh vertices to GT, and the accuracy (percentage of vertices at distance less than 5cm to the GT)

We compute these metrics using CloudCompare1 . For an easier readability, we restrict the numbers to the test set 001, results on the 224 scene being consistent the previous ones, the prediction method if the method using residual depth output with SGBM data as input.

We compare the effect of the filtering policies. We made reconstruction experiments with the two filtering methods and the baseline (no filtering) for Open-Chisel clipping range r in [START_REF] Forster | SVO: Semidirect Visual Odometry for Monocular and Multicamera Systems[END_REF][START_REF] Civera | Towards semantic slam using a monocular camera[END_REF] meters.

Table 5 presents these results for ranges 5m and 10m. The results are first computed with the full scene ground truth (including complete trees) and then with a cropped ground truth at 1m height (corresponding to the use case of autonomous lawnmower). As expected, giving all the points OpenChisel leads to better completeness but produce a lot of reconstruction artifacts, in particular at transitions between objects or sky. Filtering these points based on gradient produce much better results according to outlier production while ensuring a good completeness. On the opposite, the harder filter given by the gradient and an erosion improves even more the accuracy but drastically reduced the number of 3D points (e.g. top of the bushes and topiaries) which impacts the completeness. Better performances are achieved using a a cropped ground truth. This is mostly due to the small baseline of the stereo images and the ground view, leading to missing or uncertain tree reconstruction.

The compromise between completeness and accuracy is illustrated on figure 6. It presents, on the left, a graph showing the average distance for completeness Ground truth.

No filter, r=6m.

Erosion, r=5m. Gradient, r=10m. Finally, error maps are presented on the left side of Fig. 7. For the GT → Predictions maps, the red points (error greater than 10cm) are the missing parts. For the Predictions → GT maps, red points correspond to hallucinated objects, particularly multiple tree trunks or flowers.

Semantics

Evaluation of 3D semantics is not straightforward: there is no direct correspondence between points of the ground truth and the reconstructed mesh. We use a strategy to create a geometric clone of the ground truth and assigning to each point the label of the nearest vertex of our reconstructed mesh. Thus, we have prediction/GT label pairs usable for metric computation. Table 6 presents the results for the Multi-task Net. with gradient filtering for overall accuracy (OA), average accuracy (Av Acc.) and average intersection over union (Av. IoU). Left side of figure 7 shows snapshots of the surface with semantics labels. Most of the errors are located on the ground, mostly mixed grass and ground. The semantic prediction tend to produce smooth segmentations and fail to create very small connected component, such as pebbles in the grass.

Transfer to real data

The ultimate goal of the reconstruction pipeline is to be applicable to real data. To test the ability of our pipeline to generalize from the synthetic dataset to real outdoor data, we experiment the 3DRMS 2017 test dataset. The results are shown on Fig. 8. For visual comparison, we confront the reconstruction for the synthetic test (first row) and real data (second row), note the high difference between the two sample images. We tested first direct transfer of the neural network to the new dataset. While depth estimation was still efficient (middle image), semantic segmentation was deteriorated. To address this problem, we finetuned the segmentation decoder of the network on the train set of the 3DRMS 2017 dataset for ten epochs. Note that, in order to maintain the depth estimation quality, as the finetuning does not include depth ground truth, we froze the weights of the encoder and the depth decoder. Results are the in the right column: the semantics of the main objects and ground classes are well recovered.

Timings

We give some timing results for each separate block of our pipeline. The experiments were carried with an Intel Xeon CPU E3-1505M and Nvidia GTX1070 GPU.The stereo depth map estimation with SGBM takes 0.03s, the multi-task network depth and semantic inference takes 0.4s and the filtering step has a negligible computation time. For the 3D reconstruction, with a clipping range of 5m and with resolutions of 3cm, 5cm and 10cm, the related run times are respectively 0.4s, 0.15s and 0.1s per depth map. Higher clipping ranges increases the computation load as it significantly augment the number of points to use in the reconstruction. With a range of 10m and a resolution of 3cm, the integration of one depth map takes 2.2s. As our pipeline is designed to process incoming data online, one can expect each stereo pair to be processed in less than 0.85s for a high resolution map (3cm voxels -5m range) and in less than 0.5s for lower resolution maps (≥ 5cm voxels -5m range), which are often sufficient for autonomous navigation.

Conclusion

In this paper, we have presented a 3D reconstruction approach from multiple stereo image pairs. The reconstruction pipeline mixes both the accuracy of geometric approaches and the complex, high-order modeling made possible by the deep neural networks. We show co-learning of depth estimation and pixelwise semantic labeling is possible in robotics scenarios and improves the framework at every stage. Indeed, the multi-task network, while being lighter than separate networks, is also more effective. The proposed approach is compatible with online mapping and does not require global optimization. The method has been evaluated on the 3DRMS 2018 simulated dataset.

A close look at the reconstructed surfaces shows that most of the geometric errors come from the duplication of some objects. Moreover, the main part of semantic errors are due to mis-detected pixels which deteriorate the global score while most of the othe objects have been correctly recognized. To improve these aspects of the method, future works will include performing object detection and tracking during the sequence. First, the object identification between images would reduce the number of instances in the final product and second, labels would be regularized at object level. Even though it would make our method less suitable for online use, we also examine the possibility of a posteriori spatial regularization on the reconstructed surface such as conditional random fields.

Fig. 2 .

 2 Fig. 2. Multi-task Network architecture.

Fig. 3 .

 3 Fig. 3. Comparison of the semantic segmentation maps generated by U-Net trained on RGB images; and the proposed multi-task network architecture trained on RGB and raw depth images on input and extra refined depth on output.

Fig. 4 .

 4 Fig. 4. Comparison of performances of geometric methods (ELAS and SGBM) and the proposed approach, depth estimation through co-learning: (a) Proportion of 3D points with deviation less than a given value; (b) RMS of the depth estimates with respect to the ground truth distances.

Fig. 5 .

 5 Fig. 5. Comparison of various depth estimation approaches: geometric approach (SGBM); depth prediction (D3-Net); and the proposed multi-task network (last column). First row / Last row: depth maps / error maps in mm.

Fig. 6 .

 6 Fig. 6. Influence of the clipping range and the filtering method on 3D reconstruction.

Fig. 7 .

 7 Fig. 7. Reconstructions for two filtering policies: Semantics mesh (left) and error heat maps(right).

Table 2 .

 2 Comparison of semantic segmentation results on the 3DRMS 2018 dataset using state-of-the-art segmentation networks and the proposed Multi-task Network.

	Methods			Test on 0001		Test on 0224
		Input	Output	OA	A. Prec. A. IoU	OA	A. Prec A. IoU
	Baselines
	U-Net	RGB	S	0.9068	0.8286	0.7012	0.9054 0.7496	0.6395
	FuseNet	RGB	S	0.9091	0.8577	0.7371	0.9311 0.8038	0.7169
	Multitask refinement
	Multi-task Net. Multi-task Net.	D+S RGB+D SGBM D+S RGB+D ELAS	0.9363 0.9411	0.8916 0.8965	0.7943 0.7980	0.9277 0.7952 0.9303 0.8017	0.7107 0.7195

Table 3 .

 3 Comparisons of error metrics for depth estimation using geometric, D3-Net depth prediction and multi-task learning on the 3DRMS 2018 dataset.

		Methods		Error↓ Test 0001 Error↓ Test 0224
		Input	Output	rel	rms	rel	rms
	Geometric
	ELAS	RGBx2	D	0.526	2.140	0.444	1.993
	SGBM	RGBx2	D	0.518	1.801	0.439	1.745
	Mono image
	D3-Net Mono RGB	D	0.145	0.755	0.110	0.477
	Refinement
	FuseNet	RGB+DSGBM D	0.057	0.395	0.074	0.454
	Multi-task Net. RGB+DELAS D+S	0.079	0.410	0.066	0.414
	Multi-task Net. RGB+DSGBM D+S	0.082	0.394	0.089	0.436

Table 4 .

 4 Error measurements used for geometric reconstruction quality estimation. min q∈B d(p,q)<0.05 Completeness dist 90% for A → B min d 100 |A| p∈A 1 min q∈B d(p,q)<d = 90 Note: 1 is the indicator function.

	Metric	Definition
	Average distance for A → B	1 |A|	p∈A minq∈B d(p, q)
	Accuracy % < 5cm for A → B	100 |A|	p∈A 1

Table 5 .

 5 Evaluation of the reconstruction on the 3DRMS dataset.

			Full scene			Cropped scene z=1m
	Filtering method GT -→ Recons. Recons. -→ GT GT -→ Recons. Recons. -→ GT
		range Av. Complet.	Av.	Acc.	Av. Complet.	Av.	Acc.
			dist. dist. 90%	dist. % < 5cm	dist. dist. 90%	dist. % < 5cm
	No filtering 5m	0.061 0.145	0.201 32.2%	-	-	-	-
		10m 0.061 0.164	0.311 20.9%	-	-	-	-
	Gradient	5m	0.077 0.208	0.037 77.6%	0.042 0.102	0.027	86.3%
		10m 0.058 0.156	0.047 70.5%	0.043 0.109	0.031	83.9%
	Gradient	5m	0.128 0.427	0.024 87.2%	0.052 0.134	0.020	90.7%
	+ Erosion 10m 0.113 0.356	0.028 85.1%	0.052 0.134	0.022	89.8%

Table 6 .

 6 Evaluation of the 3D semantics. Method range Dataset OA Av. Acc. Av. IoU OA Av. Acc. Av. IoU Gradient 10m 001 0.8950 0.8735 0.7285 0.8640 0.8705 0.7339 Fig. 8. Reconstructions on the tests sets of the 3DRMS datasets 2017 and 2018.

	Full scene	Cropped scene z=1m

CloudCompare: https://www.danielgm.net/cc/