Nonlinear Modal Analysis of a Full-Scale Aircraft G. Kerschen, * M. Peeters, † and J. C. Golinval * University of Liège, 4000 Liège, Belgium and C. Stéphan ‡ ONERA The French Aerospace Lab, 92320 Châtillon, France Nonlinear normal modes, which are defined as a nonlinear extension of the concept of linear normal modes, are a rigorous tool for nonlinear modal analysis. The objective of this paper is to demonstrate that the computation of nonlinear normal modes and of their oscillation frequencies can now be achieved for complex, real-world aerospace structures. The application considered in this study is the airframe of the Morane-Saulnier Paris aircraft. Ground vibration tests of this aircraft exhibited softening nonlinearities in the connection between the external fuel tanks and the wing tips. The nonlinear normal modes of this aircraft are computed from a reduced-order nonlinear finite element model using a numerical algorithm combining shooting and pseudo-arclength continuation. Several nonlinear normal modes, involving, e.g., wing bending, wing torsion, and tail bending, are presented, which highlights that the aircraft can exhibit very interesting nonlinear phenomena. Specifically, it is shown that modes with distinct linear frequencies can interact and generate additional nonlinear modes with no linear counterpart. = acceleration at wing side, m∕s 2 x rel = relative displacement across the nonlinear connection, m y o = linear mode shape of the full size finite element model y r = linear mode shape of the reduced finite element model z = state space vector, m and m∕s z = prediction of the next state space vector, m and m∕s z 0 = initial state space vector, m and m∕s I. Introduction D URING the development of a new aircraft, testing plays a key role for flutter qualification. Before flight flutter tests, ground vibration testing (GVT) is performed on the full scale aircraft. For small aircraft for which no finite element (FE) model is developed, GVT is thus the only means of providing a modal basis. This basis is used to form a mathematical model of the structural dynamic behavior that is used for flutter predictions. For larger aircraft, the main GVTobjective is to update and validate the finite element model of the aircraft, which, in turn, will be used for making reliable flutter predictions.

Aircraft GVT has now become standard industrial practice [START_REF] Goge | Ground Vibration Testing of Large Aircraft State of the Art and Future Perspectives[END_REF][START_REF] Peeters | Modern Solutions for Ground Vibration Testing of Large Aircraft[END_REF]. However, nonlinearity is frequently encountered during these tests and represents a significant challenge for aerospace engineers [START_REF] Rodriguez Ahlquist | Assessment of Nonlinear Structural Response in A400M GVT[END_REF][START_REF] Kehoe | Aircraft Ground Vibration Testing at the NASA Dryden Flight Research Facility[END_REF]. For instance, besides the nonlinear fluid structure interaction, typical nonlinearities include backlash and friction in control surfaces and joints, hardening nonlinearities in engine to pylon connections, and saturation effects in hydraulic actuators. Although there has been some recent progress in nonlinearity detection and characterization in aerospace structures [START_REF] Gloth | Influence and Characterisation of Weak Non Linearities in Swept Sine Modal Testing[END_REF][START_REF] Gloth | Handling of Non Linear Structural Characteristics in Ground Vibration Testing[END_REF], further developments for accounting for nonlinear dynamic phenomena during aircraft certification are still needed [START_REF] Luber | Application of Analytical Methods for Advanced Flight Flutter Test[END_REF]. They would clearly improve the confidence practicing engineers have in both aircraft testing and modeling.

In this context, the nonlinear normal mode (NNM) theory offers a solid theoretical tool for interpreting a wide class of nonlinear dynamical phenomena, yet NNMs have a clear conceptual relation to the linear normal modes (LNMs) [START_REF] Vakakis | Normal Modes and Localization in Nonlinear Systems[END_REF][START_REF] Kerschen | Nonlinear Normal Modes, Part I: A Useful Framework for the Structural Dynamicist[END_REF]. However, most engineers still view NNMs as a concept that is foreign to them. One reason suppor ting this statement is that most existing constructive techniques for computing NNMs are based on asymptotic approaches [START_REF] King | An Energy Based Formulation for Computing Nonlinear Normal Modes in Undamped Continuous Systems[END_REF][START_REF] Touzé | Nonlinear Normal Modes for Damped Geometrically Nonlinear Systems: Application to Reduced Order Modelling of Harmonically Forced Structures[END_REF]. Because these techniques rely on fairly involved mathematical developments, they are usually limited to the analysis of systems with low dimensionality.

There have been few attempts to compute NNMs using numerical methods [12 15]. The objective of this paper is to demonstrate that the computation of NNMs of complex, real world structures is now within reach and that it can be of great help for understanding and interpreting experimental observations. Specifically, the algorithm proposed in [START_REF] Peeters | Nonlinear Normal Modes, Part II: Toward a Practical Computation Using Numerical Continuation Techniques[END_REF] is used in the present study to compute the NNMs of the airframe of the Morane Saulnier Paris aircraft, for which the ground vibration tests have exhibited some nonlinear structural behaviors. Very interesting nonlinear dynamics, including nonlinear modal interactions, will be revealed.

The paper is organized as follows. Section II briefly reviews the two main definitions of an NNM. Section III describes the algorithm used for NNM computation. Section IV introduces the mathematical modeling of the considered aircraft and discusses carefully the results of nonlinear modal analysis. The conclusions of the present study are summarized in Sec. V.

II. Brief Review of Nonlinear Normal Modes

A detailed description of NNMs and of their fundamental properties (e.g., frequency energy dependence, bifurcations, and stability) is given in [START_REF] Vakakis | Normal Modes and Localization in Nonlinear Systems[END_REF][START_REF] Kerschen | Nonlinear Normal Modes, Part I: A Useful Framework for the Structural Dynamicist[END_REF]. For completeness, the two main definitions of an NNM are briefly reviewed in this section.

The free response of discrete conservative mechanical systems with n degrees of freedom (DOFs) is considered, assuming that continuous systems (e.g., beams, shells, or plates) have been spatially discretized using the finite element method. The governing equations of motion are

M xt Kxt f nl fxtg 0 (1)
where M and K are the mass and stiffness matrices, respectively, xt and xt are the displacement and acceleration vectors, respectively, and f nl is the nonlinear restoring force vector.

There exist two main definitions of an NNM in the literature due to Rosenberg [START_REF] Rosenberg | Normal Modes of Nonlinear Dual Mode Systems[END_REF] and Shaw and Pierre [START_REF] Shaw | Non Linear Normal Modes and Invariant Manifolds[END_REF]:

1) Targeting a straightforward nonlinear extension of the LNM concept, Rosenberg defined an NNM motion as a vibration in unison of the system (i.e., a synchronous periodic oscillation).

2) To provide an extension of the NNM concept to damped systems, Shaw and Pierre defined an NNM as a two dimensional (2 D) invariant manifold in phase space. Such a manifold is invariant under the flow (i.e., orbits that start out in the manifold remain in it for all time), which generalizes the invariance property of LNMs to nonlinear systems.

At first glance, Rosenberg's definition may appear restrictive in two cases. First, it cannot be easily extended to nonconservative systems. However, the damped dynamics can often be interpreted based on the topological structure of the NNMs of the underlying conservative system [START_REF] Kerschen | Nonlinear Normal Modes, Part I: A Useful Framework for the Structural Dynamicist[END_REF]. Second, in the presence of internal resonances, the NNM motion is no longer synchronous, but it is still periodic.

In the present study, an NNM motion is therefore defined as a (nonnecessarily synchronous) periodic motion of the conservative mechanical system (1). As we will show, this extended definition is particularly attractive when targeting a numerical computation of the NNMs. It enables the nonlinear modes to be effectively computed using algorithms for the continuation of periodic solutions.

For illustration, the NNMs of the two DOF system

x 1 2x 1 -x 2 0.5x 3 1 0 x 2 2x 2 -x 1 0 (2)
are depicted in Fig. 1. Because of the frequency energy dependence of nonlinear oscillations, an appropriate graphical depiction is to represent NNMs in a frequency energy plot (FEP). An NNM is represented by a point, which is drawn at a frequency corresponding to the minimal period of the periodic motion and at an energy equal to the conserved total energy during the motion (i.e., the sum of kinetic and potential energies). A branch, represented by a solid line, is a family of NNM motions possessing the same qualitative features. The backbone of the plot is therefore formed by two branches, which represent in phase (S11) and out of phase (S11-) NNMs. The FEP clearly shows that the nonlinear modal parameters have a strong dependence on the total energy in the system based on the following:

1) The frequency of both the in phase and out of phase NNMs increases with the energy level, which reveals the hardening characteristic of the system.

2) The modal curves change for increasing energies. The in phase NNM tends to localize to the second DOF (i.e., it resembles a vertical curve), whereas the out of phase NNM localizes to the first DOF (i.e., it resembles a horizontal curve). Such a plot gives a very clear indication of the effects of nonlinearity on modal parameters.

III. Numerical Computation of NNMs

The numerical method proposed here for the NNM computation relies on two main techniques, namely, shooting and pseudo arclength continuation. A detailed description of the algorithm is given in [START_REF] Peeters | Nonlinear Normal Modes, Part II: Toward a Practical Computation Using Numerical Continuation Techniques[END_REF].

A. Shooting Method

The equations of motion of system (1) can be recast into state space form:

_ z gz (3)
where z x _ x is the 2n dimensional state vector and star denotes the transpose operation, and

gz _ x -M -1 Kx f nl x; _ x ( 4
)
is the vector field. The solution of this dynamical system for initial conditions z0 z 0 x 0 _ x 0 is written as zt zt; z 0 in order to exhibit the dependence on the initial conditions, z0; z 0 z 0 . A solution z p t; z p0 is a periodic solution of the autonomous system (3) if z p t; z p0 z p t T; z p0 , where T is the minimal period.

The NNM computation is carried out by finding the periodic solutions of the governing nonlinear equations of motion [START_REF] Rodriguez Ahlquist | Assessment of Nonlinear Structural Response in A400M GVT[END_REF]. In this context, the shooting method is probably the most popular numerical technique. It solves numerically the two point boundary value problem defined by the periodicity condition Hz p0 ; T ≡ z p T; z p0 -z p0 0

(5)

Hz 0 ; T zT; z 0 -z 0 is called the shooting function and represents the difference between the initial conditions and the system response at time T.

The shooting method consists of finding, in an iterative way, the initial conditions z p0 and the period T that realize a periodic motion.

Starting from some assumed initial conditions z 0 p0 , the motion z 0 p t; z 0 p0 at the assumed period T 0 can be obtained by numerical time integration methods (e.g., Runge Kutta or Newmark schemes). A Newton Raphson iteration scheme is therefore used to correct the initial guess and to converge to the actual solution.

The phase of the periodic solutions is not fixed. If zt is a solution of the autonomous system (3), then zt Δt is geometrically the same solution in state space for any Δt. Hence, an additional condition, termed the phase condition, has to be specified in order to remove the arbitrariness of the initial conditions. An isolated NNM is therefore computed by solving the augmented two point boundary value problem defined by

Fz p0 ; T ≡ Hz p0 ; T 0 hz p0 0 (6)
where hz p0 0 is the phase condition.

B. Continuation of Periodic Solutions

Because of the frequency energy dependence, the modal parameters of an NNM vary with the total energy. An NNM family, governed by Eqs. ( 6), therefore traces a curve, termed an NNM branch, in the 2n 1 dimensional space of initial conditions and period z p0 ; T. Starting from the corresponding LNM at low energy, the computation is carried out by finding successive points z p0 ; T of the NNM branch using methods for the numerical continuation of periodic motions [START_REF] Seydel | Practical Bifurcation and Stability Analysis, from Equilibirum to Chaos[END_REF][START_REF] Nayfeh | Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods[END_REF]. The so called pseudo arclength continu ation method is used herein.

Starting from a known solution z p0 j ; T j , the next periodic solution z p0;j1 ; T j1 on the branch is computed using a predictor step and a corrector step.

Predictor Step

At step j, a prediction zp0;j1 ; Tj1 of the next solution z p0;j1 ; T j1 is generated along the tangent vector to the branch at the current point z p0;j :

zp0 j1

Tj1

z p0;j T j s j p z j p T;j [START_REF] Luber | Application of Analytical Methods for Advanced Flight Flutter Test[END_REF] where s j is the predictor stepsize. The tangent vector p j p z j p T;j to the branch defined by ( 6) is solution of the system with the condition kp j k 1.

Corrector Step

The prediction is corrected by a shooting procedure in order to solve Eq. ( 6) in which the variations of the initial conditions and the period are forced to be orthogonal to the predictor step. At iteration k, the corrections z k1 p0;j1 z k p0;j1 Δz k p0;j1 and T k1 j1 T k j1 ΔT k j1 are computed by solving the overdetermined linear system using the Moore Penrose matrix inverse: 2 where the prediction is used as initial guess. This iterative process is carried out until convergence is achieved. The convergence test is based on the relative error of the periodicity condition, i.e., kHz p0 ; Tk∕kz p0 k < ϵ where ϵ is usually set to 10 -6 . An important technical remark is that the method requires the evaluation of the 2n × 2n Jacobian matrix:

∂H ∂z 0 z 0 ; T ∂zt; z 0 ∂z 0 tT -I ( 10
)
where I is the 2n × 2n identity matrix. The classical finite difference approach requires us to perturb successively each of the 2n initial conditions and integrate the nonlinear governing equations of motion. This is extremely computational intensive for large scale finite element models such as the one considered in this study. Targeting a substantial reduction of the computational cost, sensitivity analysis is exploited for determining ∂zt; z 0 ∕∂z 0 . It amounts to differentiating the equations of motion (3) with respect to the initial conditions z 0 , which leads to

d dt ∂zt; z 0 ∂z 0 ∂gz ∂z zt;z 0 ∂zt; z 0 ∂z 0 ( 11
)
with ∂z0; z 0 ∂z 0 I ( 12
)
because z0; z 0 z 0 . Hence, the matrix ∂zt; z 0 ∕∂z 0 at t T can be obtained by numerically integrating over T the initial value problem defined by linear ordinary differential equations [START_REF] Touzé | Nonlinear Normal Modes for Damped Geometrically Nonlinear Systems: Application to Reduced Order Modelling of Harmonically Forced Structures[END_REF] with initial conditions [START_REF] Slater | A Numerical Method for Determining Nonlinear Normal Modes[END_REF]. The complete algorithm is shown in Fig. 2.

IV. Nonlinear Modal Analysis of the Morane-Saulnier Paris Aircraft

The numerical computation of the NNMs of a complex, real world structure is addressed in this section. The structure is the airframe of the Morane Saulnier Paris aircraft, represented in Fig. 3. This French jet aircraft was built during the 1950s and was used as a trainer and liaison aircraft. The structural configuration under consideration corresponds to the aircraft without its jet engines and standing on the ground through its three landing gears with deflated tires. General characteristics of the aircraft are listed in Table 1. A specimen of this airplane is present in ONERA's laboratory. It is used for training engineers and technicians for GVT and for research purposes. Ground vibration tests exhibited nonlinear behavior in the connection between the wings and external fuel tanks located at the wing tip. Figure 4 shows that this connection consists of bolted attachments. The FE model of the linear part of the full scale aircraft, illustrated in Fig. 5, was elaborated from drawings at ONERA. The wings, T tail, and fuselage are modeled by means of 2 D elements such as beams and shells. Three dimensional (3 D) spring elements, which take into account the structural flexibility of the tires and landing gears, are used as boundary conditions of the aircraft. At each wing tip, the front and rear connections between the wings and the fuel tanks are modeled using beam elements. The FE model, originally created in NASTRAN, was converted into the SAMCEF FE software. The complete FE model has more than 80,000 DOFs.

The natural frequencies of the underlying linear system in the [0 50 Hz] frequency range are given in Table 2. The first nine modes correspond to aircraft rigid body modes, i.e., six modes are landing gear suspension modes, whereas three modes involve rigid body motions of the control surfaces. The frequency range of the rigid body modes is comprised between 0.09 and 3.57 Hz, i.e., noticeably lower than the first elastic mode located at 8.19 Hz. The modal shapes of different elastic normal modes of vibrations are depicted in Fig. 6. Figure 6a represents the first wing bending mode. The first and second wing torsional modes are depicted in Figs. 6b and 6c and correspond to symmetric and antisymmetric wing motions, respectively. These latter modes are of particular interest because there is a significant deformation of the nonlinear connec tions between the wings and fuel tanks. The other modes concern the aircraft tail and are consequently almost unaffected by the nonlinear connections.

Reduced Order Modeling

Because the nonlinearities are spatially localized, condensation of the linear components of the model can be achieved using the Craig Bampton reduction technique [START_REF] Craig | Coupling of Substructures for Dynamic Analysis[END_REF]. This will lead to a substantial decrease in the computational burden without degrading the computational accuracy, at least in the frequency range of interest. We stress that many practical applications possess spatially localized nonlinearities, and so the procedure developed herein is of wide applicability.

The Craig Bampton method expresses the complete set of initial DOFs in terms of retained DOFs and internal vibration modes of the primary structure clamped on the retained nodes. To introduce the nonlinear behavior of the connections the wings and the tanks, the reduced order model of the aircraft is constructed by keeping one node on both sides of the attachments. For each wing, four nodes are retained, namely, two nodes for the front attachment and two nodes for the rear attachment. In total, eight nodes of the initial FE model possessing, each 6 DOFs and 500 internal modes of vibrations, are kept in the reduced model; the FE model is thus reduced to 548 DOFs.

To assess the accuracy of the reduced model, its LNMs are compared to those predicted by the initial FE model. The comparison is performed in the space of the initial model after projecting the reduced modes back into the original space. The deviation between the mode shapes of the original y o and reduced y r models is determined using the modal assurance criterion (MAC):

MAC jy o y r j 2 jy o y o jjy r y r j (13)
MAC values range from 0 in the absence of correlation to 1 for a complete correspondence. In the [0 100 Hz] range, MAC values between modes shapes are all greater than 0.999, and the maximum relative error on the natural frequencies is 0.2%. The accuracy of this linear reduced model is therefore excellent. We note that much less than 500 internal modes are sufficient to build a good reduced model in the [0 100 Hz] frequency range. However, a larger number of modes was deliberately chosen for two main reasons. On the one hand, it serves to illustrate that our NNM algorithm can deal with systems of relatively high dimensionality. On the other hand, due to the presence of harmonics, nonlinear modal interactions may occur between a mode in the frequency range of interest and another mode outside this range. More internal modes are therefore necessary to guarantee the accuracy of the reduced model in nonlinear regimes of motion.

Aircraft Nonlinearities

The existence of a softening nonlinear behavior was evidenced during ground vibration tests conducted at ONERA. In particular, frequency response function measurements revealed a decrease in resonant frequencies when the excitation level was increased. The connections between the wings and fuel tanks were assumed to cause this observed nonlinear effect. To confirm this hypothesis, the front and rear connections of each wing were instrumented, and measure ments dedicated to nonlinearity characterization were carried out. Specifically, accelerometers were positioned on both the wing and tank sides of the connections, and two shakers were located at the tanks. The instrumentation at one of the connection is shown in Fig. 7.

The dynamic behavior of these connections in the vertical direction was investigated using the restoring force surface method [START_REF] Masri | A Nonparametric Identification Technique for Nonlinear Dynamic Problems[END_REF]. Newton's second law applied at the wing side of one connection writes

m c x c t f r;c fx c t; _ x c tg 0 ( 14
)
where f r;c is the restoring force at the considered DOF. The index c is related to the connection under consideration (i.e., either the rear or front attachment of the left or right wing). From Eq. ( 14), the restoring force is obtained by

f r;c -m c x c ( 15
)
Because the restoring force f r;c includes both the force related to the connection of interest and the force generated by the elastic deformation of the wing, a quantitative assessment of the nonlinear effects cannot be achieved precisely. However, plotting the measured acceleration signal x c t against the relative displacement and velocity across the connection can provide meaningful qualitative information about the type of nonlinearity present in the bolted connections. To this end, the aircraft was excited close to the second torsional mode, which is known to activate nonlinear behavior, using a band limited swept sine excitation. Figure 8 presents the resulting 3 D plots for the rear connections of the left and right wings. 2 D sections corresponding to zero relative velocities are also shown to highlight elastic nonlinearities. A softening elastic nonlinearity with a piecewise linear characteristic can clearly be observed. This result seems to be compatible with what was previously reported for bolted connections in the literature [START_REF] Gaul | Nonlinear Dynamics of Structures Assembled by Bolted Joints[END_REF][START_REF] Hartwigsen | Experimental Study of Nonlinear Effects in a Typical Shear Lap Joint Configuration[END_REF]. Similar nonlinear effects were also observed for the front connections. They are not shown here, because these connections participated in the aircraft response to a lesser extent.

Because nonsmooth nonlinearities require dedicated and very specific time integration methods, which are not yet available in our NNM computation algorithm, the nonlinearities in the connections were regularized using polynomial expansions. This is why the A value of -10 13 N∕m 3 was found to be adequate for the cubic spring. Damping in the connections was not considered because the focus of this study is on the NNMs of the underlying Hamiltonian system. The reduced order model, generated in the SAMCEF software, was exported to MATLAB® where nonlinearities were imple mented. Four nonlinearities were added, each one being located between the two nodes defining a nonlinear connection.

B. Nonlinear Normal Modes and Corresponding Oscillation Frequencies

The NNM computation was carried out in the MATLAB® environment using the nonlinear reduced order model built in the preceding sections. For conciseness, only the modes that are the most relevant ones for the purpose of this study are described in what follows.

The first mode examined herein is the nonlinear counterpart of the first T tail symmetric bending mode (i.e., mode 13 in Table 2). As in Sec. II, the computed backbone branch and related NNM motions are depicted in a FEP in Fig. 9. Because there is no visible wing deformation, the modal shape and the corresponding oscillation frequency remain practically unchanged when the total energy is increased in the system. Despite the fact that a nonlinear aircraft model is considered, this mode is not affected by nonlinearity and can be considered as a purely linear mode.

The nonlinear extension of the two node wing bending mode (i.e., mode 10 in Table 2) is illustrated in Fig. 10. The FEP reveals that this mode is weakly affected by the nonlinearities, at least until energies compatible with physical constraints. The frequency of the NNM motions on the backbone slightly decreases with increasing energy levels, which results from the softening characteristic of the nonlinearity. A quantitatively similar decrease in frequency between low level and high level swept sine excitations was observed during the experimental tests. Even though the frequency decrease is small, we remark that an accuracy of 0.001 Hz is usually sought during GVTs; it is therefore important to account for it. The MAC value between the NNMs at low and high energy levels (see Figs. 10b and10c) is 0.99. Accordingly, the modal shapes do not change much over the considered energy range and resemble the corresponding LNM.

Figure 11 represents the FEP of the first symmetric wing torsional mode (i.e., mode 19 in Table 2). For this mode, the relative motion of the fuel tanks is more important, which enhances the nonlinear effect of the connections. As a result, the oscillation frequency has a marked energy dependence along the backbone branch. This decay of the natural frequency is in agreement, at least qualitatively, with what observed experimentally and confirms the relevance of this nonlinear modal analysis. Conversely, the modal shapes on the backbone branch are only weakly altered by the nonlinearities; the MAC value between the NNMs at low and high energy levels is equal to 0.98.

In addition to the main backbone branch, three other NNM branches that are localized to a specific region of the FEP can be noticed. These tongues bifurcate from the backbone branch of the considered mode and bifurcate into the backbone branch of another mode, thereby realizing an internal resonance between two modes. Because of harmonics, modal interactions with distinct linear frequencies can be generated in nonlinear regimes of motion. Specifically, the 3∶1, 5∶1, and 9∶1 internal resonances in Fig. 11 involve third, fifth, and ninth harmonics of the fundamental frequency of the first symmetric wing torsional mode. A close up of the 3∶1 internal resonance is depicted in Fig. 12. Modal shapes at three different locations on the tongue are also represented. When the energy gradually increases along the tongue, a smooth transition from the first wing torsional mode to a higher frequency tail torsional mode occurs. Both modes are present in the linear modal basis, but their mixing realized by the NNM in Fig. 12c is not. It is therefore a mode with no linear counterpart. These tongues of internal resonance represent an inherently nonlinear dynamical phenomenon. Focusing now on the frequency content along the branch, a third harmonic progressively appears for increasing energies. In the vicinity of the NNM in Fig. 12c, the aircraft undergoes a subharmonic motion with two dominant frequencies. The wing vibrates at a frequency around 31 Hz, whereas the tail vibrates at a frequency close to 93 Hz. Beyond this point, the relative importance of the third harmonic component grows, until a complete coalescence with the tail torsional mode at 93 Hz is realized.

A completely similar dynamical mechanism exists for the 5∶1 and 9∶1 internal resonances, but the wing torsional mode interacts with other modes of the linear basis at 153 and 275 Hz, respectively. These modes correspond to local modes of some aircraft parts and are not shown here. We note that such nonlinear modal interactions were also observed in [START_REF] Kerschen | Nonlinear Normal Modes, Part I: A Useful Framework for the Structural Dynamicist[END_REF] for the two DOF system (2). However, they appeared at very high energies, and this is why they were not present in Fig. 1.

Finally, the second (antisymmetric) wing torsional mode (i.e., mode 20 in Table 2) is plotted in Fig. 13. The oscillation frequency undergoes a decrease around 5% in the considered energy range. The MAC value between the modal shapes at low and high energy levels is 0.97. The modal shapes are therefore affected by nonlinearity but not to a large extent. As for the first wing torsion mode, one tongue of internal resonance is present in the FEP and indicates that nonlinear modal interactions are in fact generic in such a large scale structure.

The computation of the backbone branch up to the tongue in Fig. 13 requires 20 min of CPU time with a single 2.67 GHz processor. Because of the presence of bifurcations, the computation of the branch of internal resonance demands about 1 h. This shows that the computation of NNMs of real world structures is now within reach using a conventional laptop computer.

V. Conclusions

The objective of this paper was to compute the nonlinear normal modes (NNMs) of a full scale aircraft together with the corres ponding oscillation frequencies. NNMs were represented in a frequency energy plot, in which the impact that nonlinearities have on modal parameters can readily be seen.

Parts of aircraft, such as landing gears, connections between the horizontal tail plane and vertical tail plane, and pylons, are known to be nonlinear. When these components are excited during ground vibration testing, complex nonlinear phenomena that cannot be described using linear tools may be observed. Therefore, the first practical implementation of this work is to provide a rigorous framework for the understanding and interpretation of these phenomena. For instance, the present study showed that only a few aircraft modes were influenced by the nonlinear bolted connections. Their frequencies were found to decrease by a few percent in the considered energy range, whereas their mode shapes were practically unaffected. All this information is relevant for the practicing aerospace engineer. Because of mode bifurcations, supernumerary nonlinear modes with no linear counterpart (i.e., that are not the direct extension of linear normal modes) were also reported in this study. To the authors' knowledge, this is the first time that such modes are evidenced for a real life structure. The practical realization of these nonlinear resonances is still an open question in view of the small frequency region in which they live. This will be addressed in subsequent research by applying, numerically and experimentally, a slow sweep sine excitation and nonlinear force appropriation to this aircraft.

Finally, this research paves the way for the constructive use of nonlinearity for design. Some nonlinear phenomena with no linear counterpart such as targeted energy transfer could be exploited for protecting more sensitive parts of aircraft by transferring or redistributing energy from one substructure to another substructure.

NomenclatureF

  = augmented two point boundary value problem, m and m∕s f nl = restoring force vector, N f r;c = restoring force at wing side, N g = vector field, m∕s and N∕kg H = shooting function, m and m∕s h = phase condition I = identity matrix, 1∕s K = stiffness matrix, N∕m k = linear cubic spring, N∕m k - nl = nonlinear cubic spring, N∕m 3 M = mass matrix, kg m c = mass at wing side, kg n = number of degrees of freedom p = tangent vector, m; m∕s; and s T = motion period, s T = prediction of the motion period, s x = displacement vector, m _ x = velocity vector, m∕s x = acceleration vector, m∕s 2 x c

Fig. 1

 1 Fig.1Frequency-energy plot of the two-DOF system.

Fig. 2

 2 Fig. 2 Algorithm for NNM computation.

Fig. 4 a

 4 Fig. 4 a) Connection between external fuel tank and wing tip (top view). Close-up of b) front-and c) rear-bolted attachments.

Fig. 6 aFig. 7

 67 Fig. 6 a) First wing bending mode, b) first wing torsional mode, and c) second wing torsional mode.

Fig. 5

 5 Fig. 5 Finite element model of the Morane-Saulnier Paris aircraft.

Fig. 9 FEP

 9 of an NNM involving T-tail symmetric bending.

Fig. 8

 8 Fig. 8 Left: restoring force surface, right: stiffness curves (zero relative velocity). Top: rear tank of the left wing, bottom: rear tank of the right wing.

A. Aircraft Structural Model 1 .

 1 Finite Element Model of the Underlying Linear Structure

Fig. 10 FEP

 10 Fig.10FEP of an NNM involving wing bending (linear frequency 8.19 Hz). NNMs at energy levels marked in the FEP are inset; they are given in terms of the initial displacements (m) that realize the periodic motion (with zero initial velocities assumed).

Fig. 11 FEP

 11 Fig.11FEP of an NNM involving symmetric wing torsion (linear frequency 31.08 Hz). NNMs at energy levels marked in the FEP are inset; they are given in terms of the initial displacements (m) that realize the periodic motion (with zero initial velocities assumed).

Fig. 12 3∶1

 12 Fig.123∶1 internal resonance between the first wing torsional mode and a higher tail mode. a) Close up in the FEP of the 3∶1 tongue of Fig.11. Bottom plots: NNM motions at b) beginning of the tongue (in the vicinity of the backbone of the first wing torsional mode), c) middle of the tongue, and d) extremity of the tongue. From top to bottom: NNM shapes, time series of the vertical displacements at the rear tip of the right tank ( --) and at the right side of the horizontal tail (---), and Fourier coefficients of both displacements (in gray and black, respectively).

  stiffness curve was reconstructed by considering the mathematical model f r;c kx rel k -

1 Fig. 13 FEP

 113 Fig.13FEP of an NNM involving antisymmetric wing torsion. NNMs at energy levels marked in the FEP are inset; they are given in terms of the initial displacements (m) that realize the periodic motion (with zero initial velocities assumed).

Table 1

 1 Aircraft properties

	Length, m Wingspan, m Height, m Wing area, m 2 Weight, kg
	10.4	10.1	2.6	18	1945

Table 2

 2 Natural frequencies of the linear finite element model

	Mode	Frequency,	Modal shape
	number	Hz
	1	0.0936	Rigid body mode
	2	0.7260	Rigid body mode
	3	0.9606	Rigid body mode
	4	1.2118	Rigid body mode
	5	1.2153	Rigid body mode
	6	1.7951	Rigid body mode
	7	2.1072	Rigid body mode
	8	2.5157	Rigid body mode
	9	3.5736	Rigid body mode
	10	8.1913	Two node wing bending
	11	9.8644	Fin bending
	12	16.1790	Fin torsion
	13	21.2193	First T tail symmetric bending
	14	22.7619	Front fuselage torsion
	15	23.6525	T tail torsion
	16	25.8667	Antielevator torsion
	17	28.2679	Two node vertical fuselage bending
	18	29.3309	Three mode wing bending
	19	31.0847	Symmetric wing torsion
	20	34.9151	Antisymmetric wing torsion
	21	39.5169

x (fore and aft) symmetric wing bending 22 40.8516 Antisymmetric wing torsion front fuselage 23 47.3547 Three node vertical fuselage bending 24 52.1404 Three node T tail bending


DISAT Amelia Carolina Sparavigna

Peruvian Transverse Dunes in the Google Earth Images

Keywords: Sand dunes, Transverse dunes, Peru, Barchans, Google Earth

A sand sea where winds create an endlessly moving collection of long dunes, like waves on the surface of water: this is the landscape of transverse dunes that we can see in Peru by means of satellite images of Google Earth. Here we will show some examples of these dune fields.

Introduction

Dunes are hills or ridges made of sand blown by winds, which appear in some arid and coastal regions of the world. Dunes begin to grow when some sand is deposited, because an obstacle is impeding the flow of winds. Over time, sand dunes can grow further or shrink or move, when winds are blowing in a prevailing direction, shifting of several meters in a year.

How dunes form and move is the subject of constant experimental and theoretical researches, which started several decades ago with the pioneering works of Ralph Bagnold [START_REF] Bagnold | Physics of blown sand and desert dunes[END_REF][START_REF] Pye | Aeolian sand and sand dunes[END_REF]. In fact, the current studies are increasingly addressing the problem of dune stability in relation to climate changes [START_REF] Thomas | Remobilization of southern African desert dune systems by twenty-first century global warming[END_REF][START_REF] Hiza Redsteer | Monitoring and analysis of sand dune movement and growth on the Navajo Nation, Southwestern United States[END_REF]. In nature, dunes cannot move when are stabilized by vegetation, because it reduces the impact of wind. For this reason, it is important to preserve this natural situation, or, alternatively, to create a proper environment for an artificial stabilization by planting vegetation or by using some other materials and measures [START_REF] Zhang | A comparison of different measures for stabilizing moving sand dunes in the Horqin sandy land of Inner Mongolia, China[END_REF].

When dunes are moving, their motion can be a challenge for human activities. It is therefore important a survey able of giving an estimate of their migration rate. In some recent papers, we have shown how to determine the drift of dunes using the satellite images of Google Earth, in particular its time series [START_REF] Sparavigna | A study of moving sand dunes by means of satellite images[END_REF][START_REF] Sparavigna | Moving dunes on the Google Earth[END_REF][START_REF] Sparavigna | The GNU Image Manipulation Program applied to study the sand dunes[END_REF][START_REF] Sparavigna | A case study of moving sand dunes: The barchans of the Kharga Oasis[END_REF]. In the given references, we have investigated the motion of barchans. These dunes have a crescent shape which is aligned with the wind direction. Barchans are observed in areas where there is a hard ground surface, a moderate supply of sand and a constant wind direction. Besides the crescent one, the dunes can have different shapes, depending on the amount of sand and the blowing of winds [START_REF] Vv | Types of dunes[END_REF][START_REF] Vv | Deserts and winds[END_REF]. For instance, when there is abundant supply of sand and constant wind direction, the transverse dunes are created. These are large dunes that look like sand ripples on a large scale. We will show here some beautiful fields of transverse dunes that we can find in Peru. These dunes move and, as in the case of barchans, we can study their motion using Google Earth.

Moving dunes and time series of Google Earth

A survey able of monitoring the migration of dunes can be obtained by using time series of satellite images. Since the evaluation of migration rates of dunes requires large length and time scales, these time series turn out to be the best resource. Moreover, some of them are freely available in the World Wide Web, such as those collected by Google Earth [START_REF] Sparavigna | A study of moving sand dunes by means of satellite images[END_REF][START_REF] Sparavigna | Moving dunes on the Google Earth[END_REF].

To obtain information from the satellite images we can use GIMP, the GNU Image Manipulation Program, a program which is freely available on the Web [START_REF] Sparavigna | The GNU Image Manipulation Program applied to study the sand dunes[END_REF]. In this reference we have discussed how some of GIMP features can be applied to study the dunes. To illustrate the method, here we use the Figure 1: in it we can see two examples from images of Google Earth, showing some barchans in Peru. As discussed in [START_REF] Sparavigna | The GNU Image Manipulation Program applied to study the sand dunes[END_REF], GIMP allows measuring the drift of dunes. In the Figure 2, another example of moving Peruvian barchans is proposed.

Figure 1 -Moving dunes in the time series of Google Earth images (coordinates: 15°22'55.58''S, 74°54'48.81''W, and 15°28'29.68''S, 74°36'21.66''W; Courtesy Google Earth). Using GIMP, we can measure the speed of dunes. Two images, recorded at different time, are combined in the same image. In fact, GIMP allows laying an image on the other [START_REF] Sparavigna | The GNU Image Manipulation Program applied to study the sand dunes[END_REF]. Then, using a tool of GIMP, the pair of compasses, we can measure the displacement of dunes in pixels. Since from Google Earth, we know distances on the ground, we can convert pixels in meters. In the image on the right, we can see three dunes, the smaller is moving faster than the other two [START_REF] Bagnold | Physics of blown sand and desert dunes[END_REF].

Peruvian transverse dunes

The dunes can have several shapes [START_REF] Vv | Peruvian Transverse Dunes in the Google Earth Images[END_REF]. Among the others, we have the transverse dunes, which are dunes looking like stripes of sand perpendicular to the direction of the blowing wind. The characteristic form, with gentle windward and steep leeward slopes, is maintained when they migrate. Transverse dunes usually form in the regions where wind is blowing almost constantly from a given direction. Of the other types of dunes, let us just talk about the longitudinal dunes, which are ridges with about the same slope on both sides, elongated in the direction of prevailing winds.

In the previous two figures we have shown examples of barchans. In fact, Google Earth is offering images of beautiful transverse dunes, such as those in the Province of Caraveli, Department of Arequipa, Peru. In the Figure 3, we can see the dunes near Acari: these dunes move as shown by the time series of images. Since the series has several images, we can prepare also a movie (it is available at http://dunes-on-earth.blogspot.it/2014/12/transverse-dunes-of-peru.html). In the Figure 4, it is shown another field of transverse dunes in the same province of Peru; such straight parallel ridges of moving dunes are rare to observe, being a result of the peculiar environment in which the sand moves. The huge sand sea of Caraveli is also producing amazing landscapes, such as that given in Figure 5.

Besides barchans, the proposed images are giving another example of the use of Google Earth in studying the motion of sand dunes. In a study which requires large length and time scales, Google Earth is providing a large amount of data on which we can work.

Figure 2 -

 2 Figure 2 -Another example of moving dunes in the time series of Google Earth images (coordinates: -15.038653,-74.681023, Courtesy Google Earth).

Figure 3 -

 3 Figure 3 -Transverse dunes look like waves on the surface of water. These dunes appear where there is abundant supply of sand and the direction of winds is constant. In the lower part of the figure, the images are showing the region of Acari and the dune field. In the upper part of it, we can see two images of the Google Earth time series (October 2, 2011, on the left and December 29, 2012, on the right).

Figure 4 -

 4 Figure 4 -Transverse dunes in the Caraveli Province of Peru. Note the almost perfect alignment

Figure 5 -

 5 Figure 5 -Dunes near Colca, Caraveli, Arequipa.


M Cloirec

P A Boucard

L Champaney

S Guinard

J Sen Gupta

STUDY FOR DESIGN AND IDENTIFICATION OF A BOLTED JOINT MODEL

Keywords: bolted joint, macromodel, contact, friction

The study deals with the behavior of bolted joints between the wings and the fuselage of an airplane. To study this complex structure, we use a classical finite element code with a simplified model (called macromodel) to avoid a large number of degrees of freedom implied by a fine discretization of the local geometry and non-linearities (contact, friction). During previous studies, the macromodel was constructed by simple springs. The use of this type of macromodel in assembly simulations enables to define the most loaded bolt. However, these simple springs do not allow to take into account non-linearities such as contact with friction. To face this problem, a new macromodel is considered.

Introduction

The aim of this work is to solve mechanical problems including frictional contact in bolted joint using a macromodel. This macromodel should be implemented in an usual finite element code and used for each bolted joint in the complete structure simulation. The 2-nodes macromodel element can be characterized via parameters (a gap, a tangential stiffness, an axial stiffness, a friction ratio, and a pretension for example). Once the nature of the macromodel is established, an identification of the parameters is proceed. To this purpose, simulations on a 3D structure defined as a mesomodel can be used. This structure is composed of two aluminium plates, one composite plate and a bolt. To carry out the simulations, we use a dedicated software [CHA 95] based on the LATIN method [LAD 99], which enables us to handle frictional contact between the different parts of the assembly (see figure 1). The mesomodel parameters seizure is lead through a specific multi-resolution process allowing an important wage of time [BOU 03]. In the next sections, two types of macromodel are developed: a contact macromodel and a frictional contact one. Finally, some results are presented. The contact macromodel in one dimension can be seen like a rheological model (represented by figure 2.b), and is composed of two different stiffnesses (k, k ' ) and a gap (j). The contact macromodel behavior can be defined by: F = α e kδu + α c k ' δu, where F is the force, δu the jump in displacement, k, k ' the element stiffness before and after contact and α e , α c the ratio of the jump in displacement before and after contact. Considering this behavior, the ratio of the jump in displacement after contact can be written as: α c = <|δu|-j> + |δu| and its complementary term (before contact) as: α e = 1α c . Using a standard linear discretization, the specific element behavior becomes:

{F } = α e [K]{U } + α c [K ' ]{U }; K = k -k -k k and K ′ = k ′ -k ′ -k ′ k ′ (1)
To solve the non-linear problem of contact, a standard newton algorithm is used in a quasistatic study. Consequently, the residual {R} is computed such that:

{R} = α e [K]{U + ∆U } + α c [K ' ]{U + ∆U } and the tangent matrix [K T ] is chosen such that [K T ] = [K] if α c = 0 or [K T ] = [K ' ] if α c > 0.
In a three-dimensional configuration, a simplified beam model is used to connect two shells and consequently defines a patch representing a contact joint between two parts (cf. fig. 2.a). To link the shell elements, the simplified beam model contains three translation dof and three rotation dof for each node (similar to a shell junction element integrated into a standard finite element code such as those exposed in [VAD 06], [MAY 07]). Considering the contact from a translational point of view, relations between translations and rotations implied by the element stiffness matrix are decoupled. As a result, the macromodel element stiffness matrix (K el ) is defined by:

K el = M -M -M M M ij = 0 ∀i = j M 11 = k n ; M 22 = M 33 = k t ; M 44 = M 55 = M 66 = k r (2)
where k n is the axial stiffness, k t the radial stiffness and k r the rotational stiffness. Furthermore, the terms of the joint element matrix coupling the rotations are chosen high to impose rigid body motion for rotation. To know when the contact occurs, a normed jump in displacement is computed in the local tangential plane of the specific element and compared with the gap. In this contact model, the parameters are two radial stiffnesses (k t ,k ′ t ), an axial stiffness (k n ) and a gap (j). Thus, considering the shell assembly, the joint is considered as a rigid body taken into account the translational contact.

Macromodel for frictional contact

For frictional contact (cf. fig. 2.d), three states are considered: sticking, sliding and contact. First, the plates are stuck. When a critical load (G) is attained, one plate slide over the second plate and finally the bolt contacts the plates. To model the friction, a coulomb law is used with a threshold force G = f (P c , µ), where P c is the pretension in the bolt and µ the friction ratio. The states can be expressed as:

states conditions displacement ratio sticking F t • sign( Ut ) < G α e = 1; α g = 0; α c = 0 sliding |δu t | -G k e < j and F t • sign( Ut ) ≥ G α e = G k e •|δut| ; α g = 1 -α e ; α c = 0 contact |δu t | -G k e ≥ j and F t • sign( Ut ) ≥ G α e = 1; α g = j |δut| ; α c = 1 -(α e + α g )
Consequently, the residual is computed by the relation:

{R} = α e [K e ]{U + ∆U } + α g [K g ]{U + ∆U } + α c [K c ]{U + ∆U }, and the tangent matrix ([K T ]) is chosen such that: [K T ] = [K e ] for sticking, [K T ] = [K g ] for sliding, [K T ] = [K c
] for contact, where F t is the radial force, δu t the radial jump in displacement, K e the element stiffness before contact, K g the regularization term during sliding, K c the element stiffness after contact, Ut the radial displacement speed, G the threshold friction force and α e , α g , α c the ratio of the jump in displacement related to the three states.

Conclusion

The key point of using the macromodel and the mesomodel presented is that they rest on a physical meaning. On one hand, the geometry and the initial conditions allow an identification of the macromodel's gap and pretension. On the other hand, real or virtual (mesomodel simulations presented in figure 1) elementary tests could enable to identify the macromodel's stiffnesses and friction ratio. Once the macromodel's parameters will be defined, the final aim will be to carry out a simulation on a complex structure joining the wings to the fuselage composed of plates.

  Figure 1: Bolted joint -mesomodel

Figure 2 :

 2 Figure 2: Bolted joint -macromodel definitions

Figure 3 :

 3 Figure 3: Bolted joint -macromodel results


M Cloirec

P A Boucard

L Champaney

S Guinard

J Sen Gupta

Study for design and identification of a bolted joint model

Keywords: bolted joint, macromodel, contact, friction

The study deals with the behavior of bolted joints between the wings and the fuselage of an airplane. To study this complex structure, we use a classical finite element code with a simplified model (called macromodel) to avoid a large number of degrees of freedom implied by a fine discretization of the local geometry and non-linearities (contact, friction). During previous studies, the macromodel was constructed by simple springs. The use of this type of macromodel in assembly simulations enables to define the most loaded bolt. However, these simple springs do not allow to take into account non-linearities such as contact with friction. To face this problem, a new macromodel is considered.

Introduction

The aim of this work is to solve mechanical problems including frictional contact in bolted joint using a macromodel. This macromodel should be implemented in an usual finite element code and used for each bolted joint in a complete structure simulation. The 2-nodes macromodel element can be characterized via parameters (a gap, a tangential stiffness, an axial stiffness, a friction ratio, and a pretension for example). Once the nature of the macromodel is established, an identification of the parameters is proceed. To this purpose, simulations on a 3D structure defined as a mesomodel can be used. This structure is composed of two aluminium plates, one composite plate and a bolt. To carry out the simulations, we use a dedicated software [CHA 95] based on the LATIN method [LAD 99], which enables us to handle frictional contact between the different parts of the assembly (see figure 2). In the next sections, mesomodel and macromodel are described. Finally, some results are presented.

Used strategy

The study deals with links between the wings and the fuselage of an airplane.

To study this complex structure, a classical finite element code is used with a simplified model, called macromodel, to avoid a large number of degrees of freedom implied by a fine discretization of the local geometry and nonlinearities (contact, friction).

During previous studies, the macromodel was constructed by simple springs. The use of this type of macromodel in assembly simulations enables to define the most loaded bolt. These simple springs do not allow to take into account non-linearities such as contact with friction. To face this problem, a new macromodel is considered. an example of this type of macromodel is shown figure (1). The parameters of this macromodel could be identified by the way of 3D numerical simulations of a bolted shell, called mesomodel (cf. figure 2). The 16 sub-structures 52 interfaces Fig. 2. mesomodel sample next figure (3) sum up, in a globally manner, our approach. At first, the mesomodel is set up using geometry, material properties already defined and a set of parameters (like a gap between the bolt and the holes of the plates, a friction ratio between the different parts, a pretension of the bolt, ...etc). Then, from the mesomodel's results, an evolution of the jump in displacement, between the two plates, function of the load is obtained. Those evolutions are used to realize the identification of the macromodel's parameters (stiffnesses in particular). Once the macromodel's parameters are defined, we can use it in a complete simulation of a shell assembly.

Mesomodel

Mesomodel employed

The aim of the mesomodel is to enable us to define the value of the macromodel's parameters using an identification process. In fact, the mesomodel's results define the reference solution for the identification of the macromodel's parameters.

A one-bolt mesomodel

The mesomodel is based on a test assembly which is composed of two aluminium plates, a composite plate, a screw (cf. figure 4) and a nut (cf. figure 5) made of titane. The symmetries implied by the geometry are used to define the model shown figure 7. The outer diameter of the nut and the shaft as well as the thicknesses are equal deducing medium values from initial dimensions. Besides, the length of the screw is equal to the total thickness of the joint.

The mesh contains 47 892 6-nodes prismatic elements and 29 930 nodes. The mesh of the frictional contact zone, between the plates, is fined for a better description of fricitonal contact phenomena. The material law used are linear elastic law and an orthotrope behavior is associated to the composite plate. The

Conditions limites

As it is shown figure 8, the geometry used is the quarter of the joint, and symmetry conditions (in white) are applied following the planes (x, z) and (x, y) on the different parts of the joint. On one of the tips of the composite plate (in black), a zero-displacement condition is imposed for the x -direction. The quarter of the total load F is imposed on a face of the aluminium plate (in purple).

Finally, a pretension is imposed to the screw to ensure the clamping of the joint. This pretension is defined by the way of a negative gap between the shaft and the screw (in blue). The definition of the pretension is described, more in details, in the next section 3.3.

Conditions aux interfaces

During the use of Cofast, parameters concerning the gaps, the pretension and the frictional contact zone are defined creating interfaces between the differents parts of the joint. As shown figure 8, frictional contact interfaces are identified in red, gap interface in blue and symmetries in white.

the pretension value is monitored getting back the effort at the interface between the shaft and the screw. First of all, the shaft is subjected to an elongation by the way of a negative gap between the shaft face and its respective screw face. This gap is linearly applied function of a time interval until a maximum gap. Then an evolution of the force function of the gap is obtained (cf. figure 9) and, for instance, the necessary gap to obtain a pretension of 21 kN is about -0.05 mm.

Numerical experiments

The Les résultats obtenus avec le méso-modèle sont présentés dans cette section. Ils ont été obtenus avec une charge de 40 kN, un jeu au diamètre de 0,04 mm entre le boulon et l'alésage des plaques, une précharge imposée par un jeu relatif entre la vis et le boulon de 0,05 mm (soit une valeur de précharge d'environ 21 kN et pour une valeur de coefficient de frottement de 0,3.

Les figures (10) et (11) présentent les déformées des deux plaques où on peut observer l'influence de la précharge sur les plaques engendrant une légère flex- La figure (18) présente l'évolution du saut de déplacement en fonction de la charge. Le saut de déplacement δU est défini par la valeur absolue de la différence entre un déplacement moyenné sur une zone (proche du boulon) dans l'épaisseur de la plaque composite et un déplacement moyenné sur une zone (proche du boulon) dans l'épaisseur de la plaque aluminium (cf. 19) tel Avant de décrire le comportement propre de l'élément définissant le macromodèle, la partie suivante traite de l'intégration de l'élément spécifique au sein de l'assemblage. Pour connecter les éléments plaque, l'élément spécifique est composé de 3 degrés de liberté en translation et 3 degrés de liberté en rotation (élément du même type que des éléments de jonction entre éléments plaque tels que ceux présentés dans les études développées dans [VAD 06] et [MAY 07]). On considère, dans ce macro-modèle, que le contact est pris en compte uniquement pour les mouvements de translation. Dans un premier temps, les relations entre les translations et rotations, au sein du système matriciel impliquant les raideurs, sont découplées. Par conséquent, la matrice de raideur du macro-modèle (K el ) est définie par :

0 4 0 6 0 8 0 0 1 0 2 1 6 0 . 0 - 5 0 . 0 - 4 0 . 0 - 3 0 . 0 - 2 0 . 0 - 1 0 . 0 - 0 F ∆L
δU = |U plc -U pla | (1) (a) F x (b) F y (c) F z
{F } = [K el ] {U} (2)
K el =    M -M -M M    with M =                  k n k t 0 k t k r 0 k r k r                  (3)
{∆δU} = {∆δU e } + {∆δU g } (5)
où ∆δU, ∆δU e , ∆δU g sont respectivement, l'incrément en saut de déplacement total, l'incrément en saut de déplacement adhérent,, l'incrément en saut de déplacement glissant.

Par conséquent, les différentes relations entre incrément en effort et incrément en saut de déplacement peuvent être d'éfinies : pour l'état de glissement par la relation :

{∆F } = [K e ] {∆δU e } (6)
pour l'état d'adhérence par la relation :

{∆F } = [K e ] {∆δU e } = [K g ] {∆δU g } = ([K g ] [K e ])([K g ] + [K e ]) -1 {∆δU} (7)
pour l'état de contact par la relation :

{∆F } = [K c ] {∆δU c } (8)
où K e , K g , K c et ∆δU e , ∆δU g , ∆δU c sont respectivement les raideurs et les incrément en saut de déplacementet pour les trois états, où ∆δU est l'incrément en saut de déplacement total et ∆F et l'incrément en effort. K g joue ici le rôle d'un terme de régularisation.

Concernant l'implémentation, les données accessibles de la sous-routine UEL d'Abaqus, utilisée pour définir le macromodèle, sont la matrice tangente (K T ) et le résidu (R). Dans notre cas, la matrice tangente est équivalente à la matrice de raideur de l'état d'adhérence K e et le résidu égal à la force induite par le comportement lors des trois états (F ).

L'algorithme suivant décrit le comportement de l'élément au niveau de son plan tangent et par conséquent a pour résultat de fournir l'effort dans le plan tangent. L'effort normal et les différents moments sont déduits de manière classique grâce du syst me matriciel 2 en utilisant la la matrice de raideur K el (cf. 3). Le point de départ de l'algorithme est d'évaluer le saut de déplacement tangentiel du macromodèle. La différence entre les déplacements des deux noeuds du macromodèle est tout d'abord évaluée puis normée tel que :

{δU(t)} = {U noeud2 el (t)} -{U noeud1 el (t)} (9) et δU t (t) = (δU y (t)) 2 + (δU z (t)) 2 (10)
où δU t (t) est le saut de déplacement dans le plan tangent du macromodèle au temps t. Ensuite une force prédictive F p (t) est calculée telle que :

F p (t) = K e • (δU t (t) -δU g t (t -1)) (11)
avec δU g t (t -1) la partie "glissante" du saut de déplacement dans le plan tangent au pas de temps précédent. Ensuite on calcule une force de glissement F g au pas de temps précédent telle que :

F g (t -1) = K g • δU g t (t -1) (12)
Puis on en déduit la fonction seuil (Φ(t)) telle que : Compute yield function:

Φ(t) = |F p (t) -F g (t -1)| -G (13)
Alors si Φ(t) < 0, c'est l'état d'adhérence et la force prédicitive correspond à la force au sein du macromodèle : F (t) = F p (t)), δU g t (t) = δU g t (t -1) et δU c t (t) = 0. Si Φ(t) ≥ 0, alors un multiplicateur (λ(t), équivalent à un multiplicateur plastique) est déduit de la valeur de la fonction seuil tel que :

λ(t) = Φ(t) (K e + K g ) (14)
Ensuite, le signe de la fonction seuil est évaluée telle que :

sign(Φ(t)) = |F p (t) -F g (t -1)| F p (t) -F g (t -1) (15)
et la partie glissement du saut de déplacement tel que :

δU g t (t) = δU g t (t -1) + λ(t) (16)
Alors, si celui-ci est supérieur à la moitié du jeu, soit :

|δU g t (t)| > j 2 (17)
il y a contact. Par conséquent le signe de la partie glissement du saut de déplacement est évaluée telle que :

sign(δU g t (t)) = |δU g t (t)| δU g t (t) (18)
pour définir la partie contact du saut de déplacement tel que :

δU c t (t) = (δU g t (t) - j 2 ) • sign(δU g t (t)) (19)
Enfin, la valeur du multiplicateur est réactualisée telle que : L'éprouvette sur laquelle s'appuie la définition de ce deuxième méso-modèle est constituée de deux plaques en aluminium, d'une plaque en composite, de quatre vis et de quatre 'un boulons en titane. Les symétries impliquées par la géométrie de l'assemblage sont utilisées afin de définir le modèle représenté sur la figure 29. Le diamètre extérieur de chaque tëte de vis et de chaque boulon ainsi que leur épaisseur respective sont identiques en prenant une valeur moyenne sur les diamètres initiaux. De plus, la longueur des tiges de vis est équivalente à l'épaisseur de l'assemblage. Le maillage contient 5 312 éléments prismatiques à 6 noeuds et 4 616 noeuds. Les lois matériaux utilisées sont les mêmes que pour le mésomodèle à une fixation.

λ(t) = λ(t) • 1 - |δU c t (t
Les paramètres du calcul sont décrits par le tableau 3 suivant : 3 Parameters of the 4-bolts mesomodel où µ représente le coefficient de frottement pour toutes les surfaces de contact, P n c la précharge imposée au boulon n , et j n le jeu au diamètre entre le boulon n et l'alésage lui correspondant. La charge appliquée est de 40 kN, et les conditions aux limites (symétrie et encastrement) sont les mêmes que pour le mésomodèle à une fixation.

  Fig. 1. macromodel sample

  Fig. 3. used method

Fig. 6 .

 6 Fig. 4. screw definition

  Fig. 7. used mesh

Fig. 8 .

 8 Fig. 8. Interface and boundary conditions

Fig. 9 .

 9 Fig. 9. Pretension evolution

Fig. 10 .

 10 Fig. 10. Warp (x20) aluminium plate

Fig. 16 .

 16 Fig. 16. σ yy -Composite plate

Fig. 17 .

 17 Fig. 17. Interface forces between the plates 4 Definition of the macromodel

Fig. 18 .ε

 18 Fig. 18. Curve F-δu -tension/compression loading/unloading

  Fig. 22. macromodel sample

Fig. 23 .

 23 Fig. 23. Macromodel connecting two plates

KFig. 28 .

 28 Fig. 28. Résultats liaison à une fixation

Fig. 29

 29 Fig. 29. mesh used P 1 c (kN) P 2 c (kN) P 3 c (kN) P 4 c (kN) 21.3 26 17.5 13 µ j 1 (mm) j 2 (mm) j 3 (mm) j 4 (mm) 0.3 0.03 0.035 0.04 0.045 Table3Parameters of the 4-bolts mesomodel

Fig. 30 .

 30 Fig. 30. Warp (x20) of the plates Les figures (31), (32), (33), (34) suivantes présentent les contraintes au niveau des deux plaques dans deux directions (x et y).

Fig. 31 .

 31 Fig. 31. σ xx -plaque Aluminium Les figures (35) présentent les efforts à l'interface entre les plaques. Celles-ci permettent de visualiser l'effet local de la précharge en fonction de l'effort normal (F z ) et de mesurer les efforts qui peuvent être transmis par frottement (pour F x et F y : efforts dans le plan tangent).

Fig. 32 .

 32 Fig. 32. σ yy -Aluminium plate

Fig. 36 .Fig. 37 .

 3637 Fig. 35. Interface forces between the plates l'épaisseur de la plaque aluminium (cf. 37) tel que :

  table 1 sum up the material properties employed.

		E (MPA)		Young modulus E poisson ratio ν
		aluminium	70 000		0.3
		screw and nut	110 000	0.3
	E, G (MPA)	E1	E2	E3	G12	G13	G23	ν12 ν13 ν23
	composite	90 740 37 700 37 700 18 190 18 190 13 860 0.35 0.35 0.35
	Table 1
	Material properties

  où k n est la raideur axiale, k t la raideur radiale et k r la raideur pour tous les mouvements de rotations. De plus, les valeurs des termes de la matrice de raideur du macro-modèle couplant les rotations sont choisies assez grandes pour imposer des rotations équivalentes aux deux noeuds du macro-modèle. Enfin, chacun des deux éléments plaque (connectés au macro-modèle), par l'intermédiaire de ces 4 noeuds, sont liés par un mouvement de corps rigide au noeud du macro-modèle leurs correspondant. En conclusion, en considérant un assemblage de plaques, l'ensemble constitué du macro-modèle et des deux éléments plaque qu'il relie peut être vu comme un corps rigide prenant en compte le contact en translation.

			Ft
		G = Pc . μ
	1 st cycle					Ut .
	next cycles		-G
							A-A
							initial position
	1: sticking state	2 Fig. 25. Coulomb law used	3
	2: sliding state					radial stiffnesses:
	3: contact state	f	k t 1		k t 1	= Ke -((Ke Kg)/(Ke+Kg))
	P c friction ratio: f pretension:	k t 2	1	k k t 3	k n P	= K c -((Ke Kg)/(Ke+Kg)) = Ke k k t 3 = Kg k t 2 k t 4
	axial stiffness: k n gap: j	k t 4	j			with: K e sticking state stiffness
							K g	sliding state stiffness
							K c	contact state stiffness
		Fig. 24. Joint behavior
							les plaques en
	aluminium glissent par rapport à la plaque composite. Enfin, les plaques en-
	trent en contact avec le boulon. Durant le déchargement, l'état de contact se Du point de vue du macromodèle (représenté figure 26), nous considérons
	maintient. Ensuite, lorsque l'effort au sein de la liaison devient inférieur á la trois états : adhérence, glissement et contact. Ce macromodèle comprend les
	force critique, l'état devient adhérent. La même séquence se répète lors de la paramètres suivants : trois raideurs dans le plan tangent pour les trois états
	compression. considérés (K e pour l'état d'adhérence, K g pour l'état de glissement et K c pour
	l'état de contact), une raideur normale (K n ), un jeu (j), et enfin, un coefficient Afin de représenter le comportement illustré par la figure 24, un algorithme de frottement (µ) et une précharge affectée au boulon (P c ) permettants de spécifique est utilisé. Celui-ci s'appuie sur l'algorithme développé dans [VAL 02], qui traite de problème de plasticité. Le comportement de contact avec frotte-déterminer le seui de glissement (G).
	ment est considéré d'un point de vue tangentiel, c'est à dire au niveau de la

4.3 Frictional contact behavior model

Pour les problèmes impliquant du contact avec frottement, nous considérons trois états : adhérence, glissement, contact. La figure (24) illustre le comportement de la liaison, durant ces différents états, lorsque celle-ci est soumise à une "traction-compression" durant des cycles de charge et de décharge successifs. Les cercles définis sur la figure 24 représentent une vue en coupe de la liaison au niveau de l'alésage. Tout d'abord, une traction est imposée. La précharge appliquée à la liaison permet l'adhérence entre les plaques. Celles-ci commencent par se déformer. Une fois la force critique (G) atteinte, surface de contact entre les plaques. Par conséquent, celui-ci s'applique aux déplacements tangentiels du macromodèle. Afin de représenter le frottement, une loi de Coulomb est utilisée telle que :

G = P c • µ

(

4) où G est le seuil de glissement (considéré constant dans cette étude), P c la précharge imposée au boulon et µ le coefficient de frottement entre les plaques. La loi est illustrée par la figure (25) où Ut est la vitesse de déplacement radial et F t la force radiale. Le comportement de la liaison lors des états d'adhérence et de glissement peut c Fig. 26. Macromodel employed être vu comme un comportement élasto-plastique cinématique. C'est pourquoi, durant la phase de glissement, l'incrément en saut de déplacement est divisé en une partie adhérente (élastique) et en une partie glissante (plastique) telle que :

  La figure 27 présente le cas traité. Il s'agit de soumettre une représentation de la liaison á un cycle de chargement (décrit sur la figure27) en traction compression. La liaison contient une plaque en aluminium, une plaque composite et le macromodèle. Les caratéristiques matériaux des plaques sont les même que pour le mésomodèle (cf. Tab 1).Le maillage contient 302 noeuds, 200 éléments plaque et le macromodèle. La charge appliquée est de XX kN. Les valeurs des paramètres du macromodèles sont présentés dans le tableau 2. Le nombre de pas de temps est : 128.

	Si la partie glissement du saut de déplacement est inférieure à la moitié du
	jeu, alors :
	δU c t (t) = 0		(22)
	Enfin, la force au sein du macromodèle est calculée
	F (t) = F p (t) -(λ(t) • (Ke • sign(Φ(t))) + ((Kc -Ke -Kg) • δU c t (t))(23)
	Remarque : ce comportement est basé sur un algorithme décrivant une plas-
	ticité cinématique dans le cas de la traction et de la compression. Par conséquent,
	lors d'une simulation avec ce macromodèle, la direction du chargement doit
	être constante mais le sens de chargement peut varier.
	4.4 Macromodel parameters Identification
	F
				F
		loading	3
				2
			T
				1
			Fig. 27. Problem data
			)| t (t) -δU g |δU g t (t -1)|	(20)
	et on calcule la nouvelle valeur, au signe près, de la partie glissement du saut
	de déplacement.
	δU g t (t) =	j 2	• sign(δU g t (t))	(21)

Results from macromodel

Analysis 6 Conclusion


J John
email: john@gcetly.ac.in

M S Malchijah Raj

THE RESTRAINED STEINER NUMBER OF A GRAPH

Keywords: geodetic number, Steiner number, restrained geodetic number, restrained Steiner number AMS Subject Classification: 05C12

For a connected graph G = (V, E) of order p, a set

has no isolated vertices. The minimum cardinality of a restrained Steiner set of G is the restrained Steiner number of G, and is denoted by s r (G). The restrained Steiner number of certain classes of graphs are determined. Connected graphs of order p with restrained Steiner number 2 are characterized. Various necessary conditions for the restrained Steiner number of a graph to be p are given. It is shown that, for integers a, b and p with 4 ≤ a ≤ b ≤ p, there exists a connected graph G of order p such that s(G) = a and s r (G) = b. It is also proved that for every pair of integers a, b with a ≥ 3 and b ≥ 3, there exists a connected graph G with s r (G) = a and g r (G) = b.

Introduction

By a graph G = (V, E), we mean a finite undirected graph without loops or multiple edges. The order and size of G are denoted by p and q respectively. 1 For every vertex v ∈ V , the open neighborhood N (v) is the set {u ∈ V /uv ∈ E} and the closed neighborhood of v is the set

N [v] = N (v) ∪ {v}. The degree of a vertex v ∈ V is deg G (v) = deg(v) = |N (v)|.
The minimum and maximum degrees of a graph G are denoted by δ = δ(G) and ∆ = ∆(G), respectively. The distance d(u, v) between two vertices u and v in a connected graph G is the length of a shortest u -v path in G. An u -v path of length d(u, v) is called an u -v geodesic. It is known that this distance is a metric on the vertex set V (G). For basic graph theoretic terminology we refer to [START_REF] Harary | Graph Theory[END_REF][START_REF] Buckley | Distance in Graphs, Addition-Wesley[END_REF]. For a vertex v of G, the eccentricity e(v) is the distance between v and a vertex farthest from v. The minimum eccentricity among the vertices of G is the radius, radG and the maximum eccentricity is its diameter, diamG of G. Two vertices x and y are antipodal if d(x, y) = diamG. A vertex v is an extreme vertex of a graph G if the subgraph induced by its neighbors is complete. If e = {u, v} is an edge of a graph G with d(u) = 1 and d(v) > 1, then we call e a pendant edge, u a leaf or end vertex and v a support. A geodetic set of G is a set S ⊆ V (G) such that every vertex of G is contained in a geodesic joining some pair of vertices in S. The geodetic number g(G) of G is the minimum cardinality of its geodetic sets and any geodetic set of cardinality g(G) is a minimum geodetic set. The geodetic number was introduced in [START_REF] Harary | The geodetic number of a graph[END_REF] and further studied in [START_REF] Chartrand | On the geodetic number of a graph[END_REF][START_REF] Hernando | On the Steiner, geodetic and hull number of graphs[END_REF]. A set of vertices S in a graph G is a restrained geodetic set if S is a geodetic set and the subgraph G[V -S] induced by V -S has no isolated vertices. The minimum cardinality of a restrained geodetic set, denoted by g r (G), is called the restrained geodetic number of G. A g r (G) -set is a restrained geodetic set of cardinality g r (G). The restrained geodetic number of a graph was introduced in [START_REF] Abdollahzadeh Ahangar | The restrained geodetic number of a graph[END_REF].

For a nonempty set W of vertices in a connected graph G, the Steiner distance d(W ) of W is the minimum size of a connected subgraph of G containing W . Necessarily, each subgraph is a tree and is called a Steiner tree with respect to W or a Steiner W -tree. S(W ) denotes the set of all vertices that lie on Steiner W -trees. A set W ⊆ V (G) is called a Steiner set of G if every vertex of G lies on some Steiner W -tree or if S(W ) = V (G). A Steiner set of minimum cardinality is a minimum Steiner set or simply a sset and this cardinality is the Steiner number s(G) of G. The Steiner number of a graph was introduced in [START_REF] Chartrand | The Steiner Number of a Graph[END_REF] and further studied in [START_REF] Hernando | On the Steiner, geodetic and hull number of graphs[END_REF][START_REF] Pelayo | Comment on "The Steiner number of a graph[END_REF][START_REF] Santhakumaran | The Edge Steiner Number of a Graph[END_REF][START_REF] Santhakumaran | The Upper Steiner Number of a Graph, Graph Theory Notes of New York LIX[END_REF][START_REF] Santhakumaran | The forcing Steiner number of a graph[END_REF][START_REF] Santhakumaran | The Forcing Geodetic and the Forcing Steiner Numbers of a Graph[END_REF].

The following theorems are used in the sequel.

s(G) =      2 if m = n = 1 n if n ≥ 2, m = 1 min{m, n} if m, n ≥ 2 Theorem 1.3. ([1]) Let G be a connected graph of order n ≥ 3. Then g r (G) = 2 if and only if g(G) = 2 and diamG ≥ 3 or G = K 2 ∨ H, where δ(H) ≥ 1.
1 v 2 v 3 v 4 v 5 v 6 G Figure 2.1 Example 2.2.
For the graph G given in Figure 2.1, it is clear that W = {v 1 , v 3 , v 5 } is a s -set of G and so s(G) = 3, and

W 1 = {v 1 , v 3 , v 4 , v 5 } is a s r -set so that s r (G) = 4.
(G) ≤ s r (G). Also, since V (G) is a restrained Steiner set, it is clear that s r (G) ≤ p.
Remark 2.7. The bounds in Theorem 2.6 are sharp. For any non trivial path P p (p ≥ 2), the set of two end vertices is the unique Steiner set so that s(P p ) = 2. For the cycle C 3 , s(C 3 ) = s r (C 3 ). For the complete graph K p , s r (K p ) = p. Also the inequalities are strict. For the graph G given in Theorem 2.9. There is no graph of order p with s r (G) = p -1.

Proof. Since every restrained Steiner set of G is a Steiner set of G and the complement of each restrained Steiner set has cardinality different from one, we have s r (G) = p -1.

Theorem 2.10. Let G be a non trivial tree which is not a star. Then s r (G) is equal to the set of all end vertices of G.

Proof. Let W be the set of all end vertices of G.

Then W ia a restrained Steiner set of G. By Theorem 2.4, s r (G) ≤ |W |. Since G = K 1 , p-1 , the subgraph G[V -W ] has no isolated vertices. Therefore W is a restrained Steiner set of G so that s r (G) = |W |.
Corollary 2.11. For any tree T with p ≥ 3 vertices, s r (T ) = p if and only if T is a star. Corollary 2.12. For any tree T with p ≥ 4 vertices, s r (T ) = p -2 if and only if T is a double star.

Theorem 2.13. For the cycle G = C p (p ≥ 3), s r (G) =     
p for p ∈ {3, 4, 5} 2 for p ≥ 6 and p is even 3 for p ≥ 7 and p is odd Proof. It is clear that s r (C p ) = p for p ∈ {3, 4, 5}. Let p ≥ 6. Let p be even and W = {u, v} be a set of antipodal vertices of C p . Then W is a Steiner set of G. Since the subgraph G[V -W ] has no isolated vertices, W is a restrained Steiner set of G so that s r (C p ) = 2. Let p be odd. It is clearly verified that no two element subset of C p is a Steiner set of G and so s r (G) ≥ 3. For any vertex u, let v and w be the two antipodal vertices. Let

W 1 = {u, v, w}. Then W 1 is a Steiner set of G. Since p ≥ 7, the subgraph G[V -W 1 ] has no isolated vertices. Then W 1 is a restrained Steiner set of G so that s r (C p ) = 3. Theorem 2.14. For the complete bipartite graph K m,n (m, n ≥ 2), s r (K m,n ) = m + n. Proof. Let G = K m,n (2 ≤ m ≤ n). Let X and Y be the bipartite sets of G, where X = {x 1 , x 2 , ...., x m } and Y = {y 1 , y 2 , ...., y n }. Let 2 ≤ m < n. By Theorem 1.2, W = X is a s -set of G. Since the subgraph G[V -W ] has isolated vertices, W is not a restrained Steiner set of G and so s r (G) ≥ m+1. Let W be a restrained Steiner set of G with |W | ≥ m+1. Then either W ⊆ Y or W X ∪Y . If W ⊆ Y , then W = Y is the only Steiner set of G. Since the subgraph G[V -W ] has isolated vertices, W is not a restrained Steiner set of G, which is a contradiction. If W X ∪Y , then < W > is connected. Then the Steiner W -tree contains elements of W only. Therefore W is not a restrained Steiner set of G. Hence W = X ∪Y . This implies s r (G) = m+n. Similarly if m = n, we can prove that s r (G) = m + n.
Theorem 2.15. For the wheel

W p = K 1 + C p-1 , (p ≥ 5), s r (W p ) = p -3. Proof. Let v be the vertex of K 1 and let v 1 , v 2 , ...., v p-1 , v 1 be the cycle of C p-1 . For p = 5, let W = {v 1 , v 3 } . Then W is a Steiner set of G. Since the subgraph G[V -W ] has no isolated vertices, W is a restrained Steiner set of G so that s r (W p ) = 2 = p -3. Let p ≥ 6. Let W be a Steiner set of G. If v ∈ W , then < W > is connected. Then the Steiner W -tree contains elements of W only. Therefore v / ∈ W . Hence W ⊆ V (C P -1 ). Let W = {v 1 , v 3 , v 4 , ...., v p-2 }. Then W is a restrained Steiner set of G and so s r (G) ≤ p -3. We prove that s r (G) = p -3. If not let W be a restrained Steiner set of G with |W | ≤ p -4. Then v / ∈ W . Therefore there exists at least one v i / ∈ W (1 ≤ i ≤ p -1). Then v i (1 ≤ i ≤ p -1)
does not lie on any Steiner W -tree of G and so W is not a restrained Steiner set of G, which is a contradiction. Hence s r (G) = p -3.

Theorem 2.16. Let G be a connected graph of order p ≥ 2. Then s r (G) = 2 if and only if g r (G) = 2. Proof. Let g r (G) = 2. By Theorem 1.3, g(G) = 2 and diamG ≥ 3 or G = K 2 ∨ H, where δ(H) ≥ 1. Then it follows that s r (G) = 2. Conversly let s r (G) = 2. Let W = {u, v} be a s r -set of G. Since every Steiner W -tree is a u -v geodesic and the subgraph G[V -W ] has no isolated vertices, W is a restrained geodetic set of G so that g r (G) = 2.
Corollary 2.17. For the hyper cube

Q n (n ≥ 3), s r (Q n ) = 2.
Proof. Q n has 2 n vertices, which may be labeled (a 1 a 2 a 3 ...a n ), where each

a i (1 ≤ i ≤ n) is either 0 or 1. It is easily seen that {(0, 0, 0, ... ,0), (1, 1, 1, ...,1)} is a s r (Q n ) set for n ≥ 3. Hence s r (Q n ) = 2.
Theorem 2.18. Let G be a connected graph of order p with exactly one vertex of degree p -1 which is not a cut vertex. Then s r (G) ≤ p -2.

Proof. Let x be the vertex of degree p-1 which is not a cut vertex of G. Since < N (x) > is non complete, there exist at least two non adjacent vertices say y and z such that y and z belong to < N (x) > . Let x 1 be a vertex in y -z geodesic. Then

W = V -{x, z} is a Steiner set of G. Since the subgraph G[V -W ] has no isolated vertices, W is a restrained Steiner set of G so that s r (G) ≤ p -2. v 1 v 2 v 3 v 4 x G Figure 2.2
W = {v 1 , v 3 , v 4 } is a restrained Steiner set of G so that s r (G) = 3 = p -2.
Theorem 2.20. Let G be a connected graph which contains a cut vertex of degree p -1. Then s r (G) = p. Proof. Let X = {x 1 , x 2 } and Y = {y 1 , y 2 , ...., y p-3 } be the bipartite sets of K 2,p-3 . Without loss of generality we subdivide the edge x 1 y 1 by the vertex x. If W = {x 1 , x 2 }, then the vertices x and y 1 do not lie on any Steiner W -tree of G and so W is not a restrained Steiner set of G. If W = {y 1 , y 2 , ...., y p-3 }, then the vertices x and x 1 do not lie on any Steiner W -tree of G.

Proof. Let x be a cut vertex of degree p -1. Let W be a Steiner set of G. If x ∈ W , then W = V is the only restrained Steiner set of G. If x / ∈ W , then W = V -{x} is the only restrained Steiner set of G and so s r (G) ≥ p -1. Since s r (G) = p -1, s r (G) = p. v 1 v 2 v 3 v 4 v 5 v 6 v 7 G Figure 2.3
If W ⊂ Y or W ⊂ X ∪ Y , then it is easily verified that W is not a restrained Steiner set of G. If W ⊂ Y ∪ {x}, then W is not a restrained Steiner set of G. Similarly if W ⊂ X ∪ {x}, then W is not a restrained Steiner set of G. Therefore W ⊂ X ∪ Y ∪ {x}. If W X ∪ Y ∪ {x}
x l-2 x l-3 x 2 x 1 u 1 u 2 u 3 u 4 u 2r u 2r-1 . . . . . . G Figure 2.4
It is easily verified that the eccentricity of each vertex of G is r so that radG = diamG = r. Let X = {x 1 , x 2 , ...., x l-2 } be the set of all extreme vertices of G with |X| = l -2. By Theorem 2.4, X is a subset of every restrained Steiner set of G. It is clear that X is not a restrained Steiner set of G so that s r (G) ≥ l -1. It is easily observed that X ∪ {v}, where v / ∈ X is not a restrained Steiner set of G so that s r (G) ≥ l. Let y and z be the two antipodal vertices of

x i (1 ≤ i ≤ l -2). Then W = X ∪ {y, z} is a restrained Steiner set of G and so s r (G) = l -2 + 2 = l.
Case 2. r < d. Let C 2r : v 1 , v 2 , ...., v 2r , v 1 be a cycle of order 2r and let P d-r+1 : u 0 , u 1 , u 2 , ...., u d-r be a path of order d -r + 1 and length d -r. Let H be the graph obtained from C 2r and P d-r+1 by identifying v 1 in C 2r and u 0 in P d-r+1 . Then add l -2 new vertices w 1 , w 2 , ...., w l-2 to H and join each vertex w i , (1 ≤ i ≤ l -2) to the vertex u d-r-1 and obtain the graph G of Figure 2.5. Proof. Case 1. a = b. Let P p-a+1 : u 1 , u 2 , ...., u p-a+1 be a path of order p-a+1. Let G be the graph obtained from the path P p-a+1 and the complete graph K a-1 by joining each vertex of K a-1 to the vertex u p-a+1 . The graph G is shown in Figure 2.6 and its order is p.

v r+1 v 2r v 1 u 0 v 2 u 1 u d-r-1 u d-r
Now W = V (K a-1 ) ∪ {u 1 }. By Theorem 1.1, W is a subset of every Steiner set of G and so s(G) ≥ a. Since S(W ) = V , W is a Steiner set of G and so s(G) = a. Since the subgraph G[V -W ] has no isolated vertices, W is a restrained Steiner set of G. Thus s r (G) = a. u 1 u 2 u p-a+1 K a-1 . . . G Figure 2.6 Case 2. 4 < a < b = p. Let G = K a,
p-a be the complete bipartite graph with partite sets X = {x 1 , x 2 , ...., x a } and Y = {y 1 , y 2 , ...., y p-a }. By Theorems 1.2 and 2.14, we have s(G) = a and s r (G) = b.

Case 3. 4 < a < b < p. Let G 1 = K 2,
b-a be the complete bipartite graph with partite sets X = {x 1 , x 2 } and Y = {y 1 , y 2 , ...., y b-a } and let P : v 1 , v 2 , ...., v p-b+1 be a path of order p -b + 1. Let H be the graph obtained from G 1 and P by identifying the vertex x 2 in G 1 with the vertex v 1 in P and subdividing the edge x 1 y 1 of G 1 by the vertex x. Add a -3 new vertices {z 1 , z 2 , ...., z a-3 } to H and join each vertex z i (1 ≤ i ≤ a -3) with the vertex v 2 and obtain the graph G of order p as shown in Figure 2.7.

y 1 y 2 y b-a x 1 x 2 v 1 x v 2 v 3 v p-b+1 z 1 z 2 z a-3 . . . . . . . . . G Figure 2.7
Let Z = {z 1 , z 2 , ...., z a-3 , v p-b+1 } be the set of all end vertices of G. First we show that s ( G) = a. By Theorem 1.1, Z is a subset of every Steiner set of G and so s(G) ≥ a -2. Since S(Z) = V , Z is not a Steiner set of G and so s(G) ≥ a -1. It is easily verified that Z ∪ {u} where u / ∈ Z is not a Steiner set of G and so s(G) ≥ a. Case 1. 3 ≤ a < b. Let P : z 1 , z 2 , z 3 , z 4 , z 5 , z 6 , z 7 be a path of order 7. Let P i : u i , v i , w i (1 ≤ i ≤ b -a) be a copy of path of order 3. Let H be the graph obtained from the paths P and P i (1 ≤ i ≤ b -a) by joining the vertex z 4 with each u i (1 ≤ i ≤ b -a) and joining the vertices z 1 and z 7 with each w i (1 ≤ i ≤ b -a). Let G be the graph obtained from H by adding a -2 new vertices x 1 , x 2 , ...., x a-2 and joining z 4 with each

However W = Z ∪ {x 1 , y 1 } is a Steiner set of G so that s(G) = a. Next we show that s r (G) = b. Since the subgraph G[V -W ] contains an isolated vertex, W is not a restrained Steiner set of G. By Theorem 2.4, Z is a subset of every restrained Steiner set of G. Let W 1 = Z ∪ Y . It is easily seen that W 1 is not a restrained Steiner set of G and so s r (G) ≥ b -1. It is easily observed that W 1 ∪ {v} where v / ∈ W 1 is not a restrained Steiner set of G and so s r (G) ≥ b. However W 2 = W 1 ∪ {x, x 1 } is a restrained
x i (1 ≤ i ≤ a -2). The graph G is shown in Figure 3.3. z 4 z 3 z 2 z 1 w 1 z 7 z 6 z 5 x 1 x 2 x a-2 v 1 u 1 u 2 v 2 w 2 u b-a v b-a w b-a . . . . . . . . . . . . G Figure 3.3
First we show that s r (G) = a. Let X = {x 1 , x 2 , ...., x a-2 } be the set of all end vertices of G. By Theorem 2.4, X is a subset of every restrained Steiner set of G. It is clear that X is not a restrained Steiner set of G so that s r (G) ≥ a -1. It is easily observed that X ∪ {u}, where u / ∈ X is not a restrained Steiner set of G and so s r (G) ≥ a.

However W = X ∪ {z 1 , z 7 } is a restrained Steiner set of G so that s r (G) = a. Next we show that g r (G) = b. Since the vertices u i (1 ≤ i ≤ b -a) and v i (1 ≤ i ≤ b -a)
do not lie on any geodesic joining a pair of vertices in W , W is not a restrained geodetic set of G. Let S be any restrained geodetic set of G. By Theorem 1.1, X ⊆ S. Since the vertices u Let X = {x 1 , x 2 , ...., x b-2 } be the set of all end vertices of G. Since X is a subset of every restrained geodetic set, g r (G) ≥ b-2. It is clear that X is not a restrained geodetic set of G and so g r (G) ≥ b -1. It is easily verified that X ∪ {v}, where v / ∈ X is not a restrained geodetic set of G and so g r (G) ≥ b. However S = X ∪ {u 1 , u 4 } is a restrained geodetic set of G so that g r (G) = b. Next we show that s r (G) = b + 1. By Theorem 2.4, each x i (1 ≤ i ≤ b -2) belongs to every restrained Steiner set of G and so s r (G) ≥ b -2. It is clear that X is not a restrained Steiner set of G and so s r (G) ≥ b -1. Also it is easily seen that neither X ∪{v} nor X ∪{u, v}, where u, v / ∈ X is a restrained Steiner set of G. However W = X ∪ {u 1 , u 3 , u 4 } is a restrained Steiner set of G so that s r (G) = b -2 + 3 = b + 1. Hence the proof.

i (1 ≤ i ≤ b -a) and v i (1 ≤ i ≤ b -a) lie on the w i -x j (1 ≤ i ≤ b -a, 1 ≤ j ≤ a -2) geodesics, each w i (1 ≤ i ≤ b -a) belongs to S and so g r (G) ≥ a -2 + b -a = b -2. Since the vertices z 2 , z 3 lie on z 1 -x i (1 ≤ i ≤ a -2)

Theorem 1 . 1 .

 11 ([5, 6]) Each extreme vertex of a graph G belongs to every geodetic set of G as well as every Steiner set of G . Theorem 1.2. ([6]) For the complete bipartite graph G = K m,n ,

Theorem 1 . 4 .

 14 ([1]) For any tree T different from star, g r (T ) = |L(T )|, where L(T ) is the set of all leaves of T . Moreover, g r (K 1,r ) = r + 1.Throughout the following G denotes a connected graph with at least two vertices.2 The Restrained Steiner Number of a GraphDefinition 2.1. A set W of vertices of a graph G is a restrained Steiner set if W is a Steiner set, and if either W = V or the subgraph G[V -W ] induced by V -W has no isolated vertices. The minimum cardinality of a restrained Steiner set of G is the restrained Steiner number of G, and is denoted by s r (G). A restrained Steiner set of minimum cardinality is called the s r -set of G.

  v

Figure 2 . 1 ,

 21 s(G) = 3, s r (G) = 4 and p = 6. Remark 2.8. If s(G) = p or p -1, then s r (G) = p. The converse is not true. For the cycle C 4 , s(G) = 2 = p -2 and s r (G) = 4 = p.

Remark 2 . 19 .

 219 The bound in Theorem 2.18 is sharp. For the graph G given in Figure2.2,

Remark 2 . 21 .

 221 The converse of Theorem 2.20 need not be true. For the graph G given in Figure2.3 s r (G) = 7 = p. But G has no cut vertex of degree p -1.Theorem 2.22. Let G be the graph obtained from K 2,p-3 , p ≥ 4 by subdividing an edge once. Then s r (G) = p.

  such that W contains at least one vertex in X and W contains at least one vertex in Y , then the Steiner W -tree contains the vertices adjacent to x and the elements of W only and henceS(W ) = V . Then it follows that W is not a restrained Steiner set of G. Therefore W = X ∪ Y ∪ {x}. Hence s r (G) = p. If G is C 4 or C 5 or K m,n or a connected graph which contains a cut vertex of degree p -1 or the graph obtained from K 2,p-3 , p ≥ 4 by subdividing an edge once or the graph given in Figure 2.3, then s r (G) = p. So we have the following open problem. Problem 2.23. Characterize graphs for which s r (G) = p. Theorem 2.24. For positive integers r, d and l ≥ 2 with r ≤ d ≤ 2r, there exists a connected graph G with radG = r, diamG = d and s r (G) = l. Proof. When r = 1, then d = 1 or d = 2. If d = 1, then G = K l has the desired properties. If d = 2, then G = K 1,l-1 has the desired properties. If r = d = 2, then G = K 2,l-2 has the desired properties. Now, let r ≥ 3. Construct a graph G with the desired properties as follows: Case 1. r = d. Let C 2r : u 1 , u 2 , ...., u 2r , u 1 be a cycle of order 2r. Let G be the graph obtained from C 2r by adding l -2 new vertices x 1 , x 2 , ...., x l-2 and joining each x i (1 ≤ i ≤ l -2) with the vertices u 1 and u 2 of C 2r . The graph G is shown in Figure 2.4.

w 1 w 2 w l- 2 .Then

 2 radG = r and diamG = d. Let W = {w 1 , w 2 , ...., w l-2 , u d-r } be the set of all end vertices of G. By Theorem 2.4, W is a subset of every restrained Steiner set of G. Since no vertex on the cycle C 2r lies on the Steiner W -tree, we see that W is not a restrained Steiner set of G.Let W = W ∪ {v r+1 }, where v r+1 is the antipodal vertex of u d-r in G. Then it is clear that W is a restrained Steiner set of G. Hence s r (G) = l.In view of Theorem 2.6, we have the following realization result. Theorem 2.25. Let a, b and p be integers such that 4 ≤ a ≤ b ≤ p. Then there exists a connected graph G of order p with s(G) = a and s r (G) = b.

3 . 1 . 10 G 3 . 1 Example 3 . 2 .

 31103132 Steiner set of G and so s r (G) = b. 3 The restrained Steiner and the restrained geodetic numbers of a graph Example For the graph G given in Figure 3.1, W = {v 1 , v 4 , v 6 } is a restrained Steiner set of G so that s r (G) = 3. Since the vertices v 9 and v 10 do not lie on a geodesic joining any pair of vertices of W , W is not a restrained geodetic set of G and so g r (G) ≥ 4. However W = {v 1 , v 4 , v 5 , v 6 } is a restrained geodetic set of G so that g r (G) = 4. Figure For the graph G given in Figure 3.2, W = {v 1 , v 3 , v 5 } is a restrained geodetic set of G so that g r (G) = 3. Since v 4 is a cut vertex of degree 4 = p -1, by Theorem 2.20 we have s r (G) = 5.

5 G 3 . 2 FromTheorem 3 . 3 .

 53233 Figure Examples 3.1 and 3.2, we see that there is no relation between s r (G) and g r (G). So we have the following realization result. For every pair of integers a, b with a ≥ 3 and b ≥ 3, there exists a connected graph G with s r (G) = a and g r (G) = b. Proof. For a = b, let G be K 1,a-1 . Then by Corollary 2.11 and Theorem 1.4, s r (G) = a and g r (G) = a.

  geodesic and, the vertices z 5 , z 6 lie on z 7 -x i (1 ≤ i ≤ a -2) geodesic, z 1 and z 7 belong to S and so g r(G) ≥ b -2 + 2 = b. Let S = W ∪ {w 1 , w 2 , ...., w b-a }. Then S is a restrained geodetic set of G so that g r (G) = a + b -a = b.Case 2. b < a and a = b + 1. Let P : x, y, z be a path of order 3. Add b -2 new vertices y 1 , y 2 , ...., y b-2 to the path P and join each vertex y i (1 ≤ i ≤ b -2) with the vertex y and obtain the graph H. Let G be the graph obtained from H by adding a -b -1 new vertices h 1 , h 2 , ...., h a-b-1 and joining each vertex h i (1 ≤ i ≤ a -b -1) with the vertices x, y and z. The graph G is shown in Figure 3.4.

First

  we show that g r (G) = b. Let Y = {y 1 , y 2 , ...., y b-2 } be the set of all end vertices of G. Since Y is a subset of every restrained geodetic set,g r (G) ≥ b -2. It is easily seen that Y is not a restrained geodetic set of G and so g r (G) ≥ b -1. It is easily verified that Y ∪ {v}, where v / ∈ Y is not a restrained geodetic set of G and so g r (G) ≥ b. However W = Y ∪ {x, z} is a restrained geodetic set of G so that g r (G) = b. Next we show that s r (G) = a.Since the vertices h i (1 ≤ i ≤ a -b -1) do not lie on any Steiner W -tree, W is not a restrained Steiner set of G. Since y is a full degree cut vertex of G, it follows from Therorem 2.20 that s r (G) = a. Case 3. a = b + 1. Let C : u 1 , u 2 , u 3 , u 4 , u 5 , u 6 be a cycle of order 6. Let H be the graph obtained by joining the vertices u 2 and u 6 of the cycle C. Add b -2 new vertices x 1 , x 2 , ...., x b-2 to H and join each vertex x i (1 ≤ i ≤ b -2 with the vertex u 6 and obtain the graph G as shown in Figure 3.5.

  Thus the Steiner number and the restrained Steiner number of a graph are different. Every restrained Steiner set is a Steiner set, and the converse is not true. For the graph G given in Figure 2.1, W is a Steiner set of G, however it is not a restrained Steiner set of G. Also every restrained Steiner set of G is a Steiner set of G. For the complete graph K p (p ≥ 3), s r (K p ) = p. Theorem 2.6. For any connected graph G, 2 ≤ s(G) ≤ s r (G) ≤ p. Proof. Any Steiner set needs at least two vertices and so s(G) ≥ 2. Since every restrained Steiner set is a Steiner set, it follows that s

	Remark 2.3. Since every restrained Steiner set of G is a Steiner set of G, the next
	result follows from Theorem 1.1.
	Theorem 2.4. Each extreme vertex of a connected graph G belongs to every
	restrained Steiner set of G.
	Corollary 2.5.