F Georgiades
email: fotios.georgiades@ulg.ac.be.

M Peeters
email: m.peeters@ulg.ac.be.

G Kerschen
email: g.kerschen@ulg.ac.be.

J C Golinval
email: jc.golinval@ulg.ac.be.

M Ruzzene

Modal Analysis of a Nonlinear Periodic Structure with Cyclic Symmetry

The objective of this paper is to examine nonlinear normal modes and their bifurcations in cyclic periodic structures. The nonlinear normal modes are computed using a numerical technique that combines shooting and pseudoarclength continuation. Unlike perturbation techniques, the resulting algorithm can investigate strongly nonlinear regimes of motion. This study reveals that modal interactions may occur without necessarily having commensurate natural frequencies in the underlying linear system. In addition, a countable infinity of such modal interactions are shown to exist in the system. M ODAL analysis refers to the analysis of the dynamics of a vibrating structure in terms of its modal parameters. For linear systems, these parameters, termed the mode shapes, natural frequencies, and damping ratios, can be computed from the governing equations of motion merely by solving an eigenvalue problem. A number of methods that can extract them directly from the system response have also been developed [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF]. Modal analysis and testing of linear mechanical structures has been developed over the past 40 50 years, and the techniques available today are really quite sophisticated and advanced.

Modal analysis proved useful for analyzing the dynamics of linear periodic or nearly periodic structures (i.e., structures composed of a number of similar or identical elements coupled together, practical examples of which are bladed disk assemblies and space antennas). In particular, energy localization phenomena occurring in slightly mistuned periodic structures were shown to exist due the presence of strongly localized modes [2 4]. Realizing that nonlinearity is a frequent occurrence in real life applications, the dynamics of nonlinear periodic structures were also investigated in several studies. Wei and Pierre examined the effects of dry friction on a nearly cyclic structure using the harmonic balance method [START_REF] Wei | Effects of Dry Friction Damping on the Occurrence of Localized Forced Vibrations in Nearly Cyclic Structures[END_REF]. In a series of papers, Vakakis et al. demonstrated that, in contrast to the findings of linear theory, nonlinear mode localization may occur in perfectly cyclic nonlinear systems [6 9]. Other studies dealing with mode localization in nonlinear cyclic systems are [10 12].

Because any attempt to apply traditional linear analysis to nonlinear systems is bound to fail, the nonlinear normal mode (NNM) theory is used herein. NNMs and asymptotic analysis were used in [6 9] to examine the free and forced vibrations of nonlinear periodic systems with cyclic symmetry. In the present study, we revisit the free dynamics of such systems using the numerical algorithm for NNM computation proposed in [START_REF] Kerschen | Nonlinear Normal Modes, Part I: A Useful Framework for the Structural Dynamicist[END_REF][START_REF] Peeters | Nonlinear Normal Modes, Part II: Practical Computation Using Numerical Continuation Techniques[END_REF]. This computational framework allows the relaxation of the assumptions of weak nonlinearity and/or weak coupling used in [6 9]. In addition, it is shown that modal interactions can take place in nonlinear systems without necessarily having commensurate natural frequencies in the underlying linear system. This dynamical phenomenon is rarely discussed in the literature.

This paper is organized as follows. In the next section, the two main definitions of an NNM are provided. The fundamental properties of NNMs are also described, and a numerical algorithm for their computation is introduced. In Sec. III, the NNM theory is applied to a nonlinear cyclic system, and the resulting dynamics are analyzed in detail. The conclusions are summarized in Sec. IV.

II. Nonlinear Normal Modes

The free response of discrete mechanical systems is considered, assuming that continuous systems have been spatially discretized using, for example, the finite element method. The equations of motion are

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

A. Definitions

NNMs were pioneered in the 1960s thanks to the seminal work of Rosenberg [START_REF] Rosenberg | On Nonlinear Vibrations of Systems with Many Degrees Of Freedom[END_REF]. They were regarded as a theoretical curiosity until the beginning of the 1990s, when they were given a new impetus through the efforts of Vakakis et al. [16 18] and Shaw and Pierre [START_REF] Shaw | Non Linear Normal Modes and Invariant Manifolds[END_REF][START_REF] Shaw | Normal Modes for Non Linear Vibratory Systems[END_REF].

Two main definitions of an NNM exist in the literature. During the normal mode motion of a linear conservative system, each system component moves with the same frequency and with a fixed ratio amongst the displacements of the components. Targeting a straightforward nonlinear extension of the linear normal mode (LNM) concept, Rosenberg defined an NNM as a vibration in unison of the system (i.e., a synchronous oscillation). This definition requires that all material points of the system reach their extreme values and pass through zero simultaneously and allows all displacements to be expressed in terms of a single reference displacement. One fundamental difference with linear systems is that the modal lines of nonlinear systems are generally curves, resulting from the nonlinear relationship between the coordinates during the periodic motion.

Shaw and Pierre proposed a generalization of Rosenberg's definition that provides a direct and elegant extension of the NNM concept to damped systems [START_REF] Shaw | Normal Modes for Non Linear Vibratory Systems[END_REF]. Based on geometric arguments and inspired by the center manifold technique, they defined an NNM as a two dimensional 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 extends the invariance property of LNMs to nonlinear systems.

At first glance, Rosenberg's definition may appear restrictive in two cases:

1) It cannot be easily extended to nonconservative systems.

2) In the presence of internal resonances (i.e., when two or more NNMs interact), some coordinates may have a dominant frequency component that is different than that of the other coordinates (e.g., some coordinates may vibrate faster than others). In this case, the system no longer vibrates in unison.

However, these two limitations can be circumvented. First, as shown in [START_REF] Kerschen | Nonlinear Normal Modes, Part I: A Useful Framework for the Structural Dynamicist[END_REF][START_REF] Lee | Complicated Dynamics of a Linear Oscillator with a Light, Essentially Nonlinear Attachment[END_REF], the damped dynamics can often be interpreted based on the topological structure and bifurcations of the NNMs of the underlying undamped system. Secondly, realizing that the motion is still periodic in the presence of internal resonances, Rosenberg's definition of an NNM can be extended to a (nonnecessarily synchronous) periodic motion of the system. This extended definition is particularly attractive when targeting a numerical computation of the NNMs. As will be shown in Sec. II.C, the definition also enables the nonlinear modes to be effectively computed using advanced algorithms for the continuation of periodic solutions. This NNM definition is considered throughout the present study.

B. Fundamental Properties

A detailed description of NNMs and their fundamental properties is beyond the scope of this paper. However, for completeness, the frequency energy dependence, as well as bifurcations and stability of NNMs, are briefly reviewed in this section. For further details, the reader can refer to [START_REF] Kerschen | Nonlinear Normal Modes, Part I: A Useful Framework for the Structural Dynamicist[END_REF][START_REF] Vakakis | Normal Modes and Localization in Nonlinear Systems[END_REF][START_REF] Vakakis | Non Linear Normal Modes and Their Applications in Vibration Theory: An Overview[END_REF].

A typical dynamical feature of nonlinear systems is the frequency energy dependence of their oscillations. One important consequence is that the frequency response functions (FRFs) of nonlinear systems are no longer invariant. The modal curves and frequencies of oscillation of NNMs also depend on the total energy in the system. In contrast to linear theory, this energy dependence prevents the direct separation of space and time in the governing equations of motion, which complicates the analytical calculation of the NNMs.

Because of frequency energy dependence, an appropriate graphical depiction of the NNMs is key to their exploitation. In this study, the representation of the NNMs of the underlying Hamiltonian system in a frequency energy plot (FEP) is adopted. An NNM is represented by a point in the FEP, 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. A branch, represented by a solid line, corresponds to a family of NNM motions possessing the same qualitative features (e.g., the in phase NNM motions of a 2 degree of freedom (DOF) system). For illustration, a 2 DOF system, the motion of which is governed by the equations

x 1 2x 1 x 2 0:5x 3 1 0 x 2 2x 2 x 1 0 (2)
is considered, and its FEP is represented in Fig. 1. The backbone of the plot is formed by two branches, which represent in phase (S11) and out of phase (S11 ) synchronous NNMs. The FEP clearly shows that the nonlinear modal parameters have a strong dependence on the total energy in the system:

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) 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 an horizontal curve).

A second salient feature of nonlinear systems is that NNMs may interact during a general motion of the system. An energy exchange between the different modes involved may therefore be observed during the internal resonance. For instance, exciting a high frequency mode may produce a large amplitude response in a low frequency mode. Internally resonant NNMs have no counterpart in linear systems and are generated through bifurcations. The FEP in Fig. 1 does not seem to feature internally resonant NNMs. At higher energies, however, another branch of periodic solutions, termed a tongue, emanates from the backbone branch S11 (Fig. 2). On this tongue, denoted S31, there is a 3:1 internal resonance between the in phase and out of phase NNMs. We note that other resonance scenarios (e.g., 2:1, 4:1, 5:1, etc.) exist in this seemingly simple system and are described in [START_REF] Kerschen | Nonlinear Normal Modes, Part I: A Useful Framework for the Structural Dynamicist[END_REF].

A third fundamental property of NNMs is that their number may exceed the number of DOFs of the system. Because of mode bifurcations, not all NNMs can be regarded as a nonlinear continuation of normal modes of linear systems. Internally resonant NNMs are one possible example. Moreover, NNMs can be stable or Fig. 1 Frequency-energy plot of system [START_REF] Pierre | Localization of Vibrations by Structural Irregularity[END_REF]. NNM motions depicted in the configuration space are inset. The horizontal and vertical axes in these plots are the displacements of the first and second DOFs, respectively; the aspect ratio is set so that increments on the horizontal and vertical axes are equal in size.

unstable, which is in contrast to linear theory, in which all modes are neutrally stable. In this context, instability means that small perturbations of the initial conditions that generate the NNM motion lead to the elimination of the mode oscillation. Therefore, unstable NNMs are not physically realizable. In the present study, the stability analysis is performed numerically through the eigenvalues of the monodromy matrix. Bifurcations and stability are interrelated concepts, because a change in stability necessarily occurs through a bifurcation.

NNMs have a fundamental limitation compared with their linear counterpart, that is, the general motion of a nonlinear system cannot be expressed as a superposition of individual NNM motions. The practical utility of NNMs might appear, at first, questionable. A first motivation to compute and exploit the NNMs is that forced resonances in nonlinear systems occur in their neighborhoods. The knowledge of the NNMs can therefore provide valuable insight into the structure of the resonances, a feature of considerable engineering importance [START_REF] Vakakis | Normal Modes and Localization in Nonlinear Systems[END_REF]. A second motivation is that NNMs are attractors for the damped dynamics. For example, during the targeted energy transfer phenomenon [START_REF] Vakakis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF], the NNMs attract the free dynamics, even though the motion is initiated using arbitrary initial conditions. A third motivation is that an experimental structure can be excited in one of its NNMs provided an appropriate exciter force pattern is applied. Nonlinear normal modes can therefore be isolated using force appropriation techniques, similar to what is achieved for linear structures (e.g., ground vibration testing of aircraft).

C. Numerical Computation

There have been very few attempts to compute NNMs using numerical methods (see, for example, [21,23 25]). Interestingly, algorithms for the continuation of periodic solutions are really quite sophisticated and advanced [START_REF] Doedel | Software for Continuation and Bifurcation Problems in Ordinary Differential Equations[END_REF][START_REF] Seydel | Practical Bifurcation and Stability Analysis, from Equilibrium to Chaos[END_REF], yet they have not been fully exploited for the computation of nonlinear modes.

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

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, with denoting the transpose operation, and

g z _ x M 1 Kx f nl 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 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] Bendiksen | Localization Phenomena in Structural Dynamics[END_REF]. In this context, the shooting method, probably one of the most popular numerical techniques, solves numerically the two point boundary value problem defined by the periodicity condition:

H z p0 ; T z p T; z p0 z p0 0 (5)
where 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. Unlike forced motion, the period T of the free response is not known a priori.

The shooting method consists of finding, in an iterative way, the initial conditions z p0 and the period T that realize a periodic motion. To this end, the method relies on direct numerical time integration and on the Newton Raphson algorithm.

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). In general, the initial guess z 0 p0 ; T 0 does not satisfy the periodicity condition [START_REF] Wei | Effects of Dry Friction Damping on the Occurrence of Localized Forced Vibrations in Nearly Cyclic Structures[END_REF]. A Newton Raphson iteration scheme is therefore used to correct the initial guess and to converge to the actual solution. The corrections z 0 p0 and T 0 are found by expanding the nonlinear function

H z 0 p0 z 0 p0 ; T 0 T 0 0 (6)
in a Taylor series and neglecting higher order terms. 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 phase space for any t. Hence, an additional condition, termed the phase condition, has to be specified to remove the arbitrariness of the initial conditions. This is discussed in detail in [START_REF] Peeters | Nonlinear Normal Modes, Part II: Practical Computation Using Numerical Continuation Techniques[END_REF].

In summary, an isolated NNM is computed by solving the augmented two point boundary value problem defined by

F z p0 ; T Hz p0 ; T 0 hz p0 0 (7)
where hz p0 0 is the phase condition.

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. [START_REF] Vakakis | Mode Localization in a Class of Multidegree of Freedom Nonlinear Systems with Cyclic Symmetry[END_REF], 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 (also called path following methods) [START_REF] Seydel | Practical Bifurcation and Stability Analysis, from Equilibrium to Chaos[END_REF]. The space z p0 ; T is termed the continuation space.

Among the different methods for numerical continuation found in the literature, the so called pseudoarclength continuation method is used herein.

Starting from a known solution z p0;j ; T j , the method consists of computing the next periodic solution z p0 j1 ; T j1 on the branch using a predictor step and a corrector step.

1) Predictor step: At step j, a prediction z p0;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 : where s j is the predictor step size. The tangent vector p j p T z;j p T j T to the branch defined by Eqs. [START_REF] Vakakis | Mode Localization in a Class of Multidegree of Freedom Nonlinear Systems with Cyclic Symmetry[END_REF] is the solution of the system: @H @z p0 j z p0;j ;T j @H @T j z p0;j ;T j @h @z p0 j z p0;j 0

zp0;j1 Tj1 z p0 j T j s j p z;j p T;j (8)
" # p z;j p T j 0 0 (9)
with the condition kp j k 1. This normalization can be taken into account by fixing one component of the tangent vector and solving the resulting overdetermined system using the Moore Penrose matrix inverse; the tangent vector is then normalized to 1.

2) Corrector step: The prediction is corrected by a shooting procedure to solve Eqs. [START_REF] Vakakis | Mode Localization in a Class of Multidegree of Freedom Nonlinear Systems with Cyclic Symmetry[END_REF], 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 T k1 j1 T k j1 T k j1 ( 10
)
are computed by solving the overdetermined linear system using the Moore Penrose matrix inverse:

@H @z p0 j z k p0;j1 ;T k j1 @H @T j z k p0;j1 ;T k j1 @h @z p0 j z k p0;j1 0 p T z;j p T j 2
T k j1 " # H z k p0;j1 ; T k j1 h z k p0;j1 0 2 6 6 4 3 7 7 5 (11)
where the prediction is used as initial guess, that is, z 0 p0;j1 zp0 j1 and T 0 j1 Tj1 . The last equation in Eq. ( 11) corresponds to the orthogonality condition for the corrector step.

This iterative process is carried out until convergence is achieved. The convergence test is based on the relative error of the periodicity condition:

kHz p0 ; Tk kz p0 k kz p T; z p0 z p0 k kz p0 k < ( 12
)
where is the prescribed relative precision.

Algorithm for NNM computation

The algorithm proposed for the computation of NNM motions is a combination of the shooting and pseudoarclength continuation methods, as shown in Fig. 3. It has been implemented in the MATLAB® environment. Other features of the algorithm, such as the step control, the reduction of the computational burden, and the method used for the numerical integration of the equations of motion, are discussed in [START_REF] Peeters | Nonlinear Normal Modes, Part II: Practical Computation Using Numerical Continuation Techniques[END_REF].

Thus far, the NNMs have been considered as branches in the continuation space z p0 ; T. As explained in Sec. II.B, the NNMs are represented in an FEP in this study. This plot can be computed in a straightforward manner: 1) the conserved total energy is computed from the initial conditions realizing the NNM motion, and 2) the frequency of the NNM motion is calculated directly from the period.

III. Modal Analysis of a Nonlinear Cyclic Structure

The objective of this section is to determine the essential features of the modal parameters of nonlinear periodic structures with cyclic symmetry. To this end, the modal shapes and frequencies of oscillation of the different NNM branches are computed from the free response of the undamped system.

The system considered throughout this paper is a simplified mathematical model of a bladed disk assembly. This model, though simplified, allows the investigation of interesting and complex dynamic phenomena related to the presence of nonlinearity. The lumped parameter model admits a single degree of freedom for each blade and includes a similarly simplified representation of the flexible disk. The bladed disk is composed of 30 sectors assembled with cyclic periodicity; a single sector is represented in Fig. 4. Each sector is modeled using disk (M) and blade (m) lumped masses coupled by linear (k) and cubic (k nl ) springs. The nonlinear springs can, for instance, be representative of geometrically nonlinear effects in slender blades. The disk masses are connected together by linear springs, K. The equations of motion of this 60 DOF system are LNM Shooting NNM motion at low energy z p0, [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF] , T (1) (j = 1)

Predictor step

Computation of the tangent and control of the stepsize

z (0) p0,(j +1) , T (0) (j +1) (k = 0)
Corrector step (Shooting)

Residual evaluation

Time integration:

H = H (z (k ) p0,(j +1) )
Cvg ?

|| H || < ε|| z (k ) p0,(j +1) ||
Computation of the corrections:

z (k +1) p0,(j +1) = z (k ) p0,(j +1) + ∆ z (k ) p0 , (j +1) T (k +1) (j +1) = T (k ) (j +1) + ∆ T (k ) (j +1) k = k + 1 NO YES T (j +1) = T (k ) (j +1) z p0,(j +1) = z (k ) p0,(j +1) j = j + 1
Next NNM motion on the branch:

x i kx i X i k nl x i X i 3 0 M X i KX i X i1 KX i X i 1 kX i x i k nl X i x i 3 0 (13)
for i 1; . . . ; 30, where X 31 X 1 and X 0 X 30 (conditions of cyclic periodicity). X i and x i are the displacements of the disk and blade masses of the ith sector, respectively. The values M 1, m 0:3, K 1, k 1, and k nl 0:1 are used in this study.

A. Modal Analysis of the Underlying Linear System

Before studying the nonlinear bladed disk assembly, the natural frequencies and mode shapes of the underlying linear system are first discussed. All bladed assemblies with circumferential symmetry exhibit single and double modes [START_REF] Ewins | Structural Dynamics Characteristics of Bladed Assemblies[END_REF].

Double modes represent the majority. They have the same natural frequency and similar mode shapes. In fact, no unique mode shapes can be specified for these modes. Rather, it is sufficient to specify two suitably orthogonal shapes and to note that, when vibrating freely at that natural frequency, the structure can assume any form given by a linear combination of the two specified shapes. Further, at the corresponding natural frequency, the assembly can vibrate in any combination of cos n and sin n circumferential distributions of displacement around the assembly, that is, in a shape of the form cos n . The mode shape is characterized by n nodal diameters because the displacement is constrained to be zero along n equally spaced diametral lines. The mode shapes of a mode pair have mutually orthogonal nodal diameters.

Single modes correspond to motion with all the blades undergoing either in phase (zero nodal diameter) or out of phase (N=2 nodal diameters) motion.

The natural frequencies of the underlying linear bladed assembly are listed in Table 1, in which the modes are denoted by the integer pair n; m, which corresponds to the number of nodal circles n and nodal diameters m for the considered mode. In the model [START_REF] Kerschen | Nonlinear Normal Modes, Part I: A Useful Framework for the Structural Dynamicist[END_REF], the nodal circle parameter n can only take the values of n 0 or n 1 according to whether the blade and disk masses undergo in phase or out of phase motion, respectively. One observes the existence of 28 pairs of double modes and 4 single modes. Figure 5 representative LNMs of the bladed assembly, namely, mode (0,15) and one mode of the mode pairs (0,1), [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Wei | Effects of Dry Friction Damping on the Occurrence of Localized Forced Vibrations in Nearly Cyclic Structures[END_REF] and [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Peeters | Nonlinear Normal Modes, Part II: Practical Computation Using Numerical Continuation Techniques[END_REF].

depicts four

B. Nonlinear Normal Modes

Modal analysis of the nonlinear bladed assembly is carried out in this section using the algorithm described earlier. NNM branches are computed by starting from the corresponding LNMs at low energy and gradually increasing the total energy in the system. These branches, termed backbone branches, are represented in Fig. 6 and form the skeleton of the FEP. As we shall see, other NNM branches bifurcate from and coalesce into these backbone branches.

The first noticeable feature in Fig. 6 is the frequency energy dependence of the NNMs. The oscillation frequency of the modes with one nodal circle is strongly affected by the nonlinearities in the system. For these modes, the blade and disk masses vibrate in an out of phase fashion, which enhances nonlinear effects. On the other hand, the oscillation frequency of the modes with zero nodal circles is much less affected. This is because the blade and disk masses vibrate in an in phase fashion for these modes.

Similar and Nonsimilar Nonlinear Normal Modes

In addition to the dependence of their oscillation frequency, the NNMs may also have their modal shapes vary with the total energy in the system. According to Rosenberg's terminology [START_REF] Rosenberg | On Nonlinear Vibrations of Systems with Many Degrees Of Freedom[END_REF], a similar NNM corresponds to an (energy independent) straight modal line in the configuration space and occurs only in systems presenting certain spatial symmetries. A nonsimilar NNM corresponds to a curve in the configuration space, the shape of which varies with the total energy. Because of its symmetry properties, the system possesses both similar and nonsimilar NNMs. Two examples of similar NNMs in the bladed disk are the nonlinear extension of the LNMs with zero nodal diameters, namely, modes (0,0) and (1,0). Mode (0,0) is a rigid body mode, which is obviously unaffected by nonlinearity. The FEP of mode (1,0) in Fig. 7 clearly depicts that, although the NNM frequency is altered by the nonlinearities in the system, the modal shape remains unchanged.

Nonsimilar NNMs resemble the corresponding LNMs at low energy. The structure (i.e., the number of nodal circles and diameters) is preserved and, as for the modes of the linear system, they mostly appear in pairs. Nonsimilar NNMs in this system are either weakly, moderately, or strongly affected by nonlinearity for increasing energy levels. Figure 8 represents a mode of the mode pair (0,2), whose shape is almost energy independent. Figure 9 shows that the NNM motions of mode pair (0,14) have a marked energy dependence.

A remarkable property of the NNM motions of mode [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Peeters | Nonlinear Normal Modes, Part II: Practical Computation Using Numerical Continuation Techniques[END_REF] is that the vibrational energy localizes to a limited number of sectors (four in this case), the rest of the system being virtually motionless (see Fig. 10). The resulting spatial confinement of the energy causes the responses of some blades to become dangerously high and might lead to premature failure of the blades. For illustration, the time series corresponding to such an NNM motion is displayed in Fig. 11. This localization phenomenon was also observed in linear mistuned bladed assemblies [START_REF] Castanier | Modeling and Analysis of Mistuned Bladed Disk Vibration: Status and Emerging Directions[END_REF], but, here, it occurs even in the absence of structural disorder and direct interblade coupling. Localization is in fact the result of the frequency energy dependence inherent to nonlinear oscillations, as discussed in [START_REF] Vakakis | A Multiple Scales Analysis of Nonlinear, Localized Modes in a Cyclic Periodic System[END_REF].

All these NNM motions correspond to standing wave motions in the sense that the system coordinates vibrate in a synchronous manner. They are represented by lines or curves in the configuration space. The phase condition used for their computation assumes that all initial velocities are zero. One therefore starts the motion from a maximum of the potential energy (see, for example, Fig. 11).

Modal Interaction: Internally Resonant Nonlinear Normal Modes

A first example of modal interaction is the 1:1 internal resonance between the two modes of a mode pair. This resonance scenario results in NNM motions, which take the form of traveling waves and are represented by ellipses in the configuration space. A detailed

Frequency (rad/s) Fig. 6 Evolution of the NNM frequencies with the total energy in the system. The solid lines represent an NNM with zero nodal circles and the dashed lines represent an NNM with 1 nodal circle. Fig. 8 FEP of one mode of the mode pair (0,2). NNMs represented by bar graphs are inset; they are given in terms of the initial displacements that realize the periodic motion (with zero initial velocities assumed). The blade and disk masses are shown in black and gray, respectively. Frequency (rad/s) Fig. 14 A 3:1 internal resonance between modes of the mode pairs (0,6) and [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Samaranayake | Resonant Vibrations in Harmonically Excited Weakly Coupled Mechanical Systems with Cyclic Symmetry[END_REF]. The solid line corresponds to the backbone of one mode of the mode pair (0,6), which is continued by a tongue of internally resonant NNMs. The dashed line corresponds to the backbone of one mode of the mode pair [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Samaranayake | Resonant Vibrations in Harmonically Excited Weakly Coupled Mechanical Systems with Cyclic Symmetry[END_REF] represented at the third of its dominant frequency. analytical study of these modes is given in [START_REF] Vakakis | Dynamics of a Nonlinear Periodic Structure with Cyclic Symmetry[END_REF]. Because of the existence of a phase difference between the coordinates, a different phase condition is considered for the NNM computation: only one initial velocity is set to zero, which is compatible with a traveling wave motion. For instance, Fig. 12 depicts the NNMs corresponding to a 1:1 internal resonance between the modes of mode pair [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Bendiksen | Localization Phenomena in Structural Dynamics[END_REF]. The mode structure is preserved in the sense that this traveling wave motion also features one nodal circle and three nodal diameters. Representative time series are shown in Fig. 13 and clearly highlight that the motion is no longer synchronous. In this particular case, the traveling wave is propagating in the counterclockwise direction, but its companion propagating in the clockwise direction also exists. Fig. 16 Intricate succession of modal interactions between mode (0,6) and other modes of the system.

x 10 Fig. 17 Splitting of one mode of the mode pair [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Wei | Effects of Dry Friction Damping on the Occurrence of Localized Forced Vibrations in Nearly Cyclic Structures[END_REF]. NNMs represented by bar graphs are inset; they are given in terms of the initial displacements that realize the periodic motion (with zero initial velocities assumed). The blade and disk masses are shown in black and gray, respectively. These modes have an important practical significance because they can be excited with an appropriate engine order excitation. Other resonance scenarios can be observed in this system with the occurrence of tongues of internally resonant NNMs (see Sec. II.B). Unlike backbone branches, tongues are localized to a specific region of the FEP. They bifurcate from the backbone branch of a specific mode and coalesce into the backbone branch of another mode, thereby realizing an internal resonance between the two modes. For instance, Fig. 14 depicts a 3:1 internal resonance between modes (0,6) and [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Samaranayake | Resonant Vibrations in Harmonically Excited Weakly Coupled Mechanical Systems with Cyclic Symmetry[END_REF] in the FEP. To better understand the resonance mechanism, the backbone of mode [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Samaranayake | Resonant Vibrations in Harmonically Excited Weakly Coupled Mechanical Systems with Cyclic Symmetry[END_REF] is represented at the third of its characteristic frequency (this is relevant because a periodic solution of period T is also periodic with period 3T). This demonstrates that a smooth transition from mode (0,6) to mode [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Samaranayake | Resonant Vibrations in Harmonically Excited Weakly Coupled Mechanical Systems with Cyclic Symmetry[END_REF] occurs on the tongue. A further illustration is that modes M1 and M2, which are the modes right after and before the coalescence of the two NNM branches, are almost identical.

During this 3:1 internal resonance, the system vibrates along a subharmonic NNM, that is, an NNM motion characterized by more than one dominant frequency component. On the branch of mode (0,6), the motion is characterized by one dominant frequency com ponent, say, !. As we move along the tongue from this branch, a third harmonic progressively appears, and the system vibrates with two dominant frequency components, ! and 3!. As we progress further on the tongue, the third harmonic tends to dominate the component at the fundamental frequency until this latter completely disappears. At this precise moment, a transition to mode [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Samaranayake | Resonant Vibrations in Harmonically Excited Weakly Coupled Mechanical Systems with Cyclic Symmetry[END_REF] is realized. This transition is illustrated in Fig. 15 using a time series representative of the NNM motion at three different locations on the tongue.

Surprisingly, the ratio of the linear natural frequencies of modes (0,6) and (1,12) is far from 3; it is equal to 2.47. A 3:1 internal resonance between the two modes can still be realized, because the frequency of mode (0,6) increases much less rapidly than that of mode [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Samaranayake | Resonant Vibrations in Harmonically Excited Weakly Coupled Mechanical Systems with Cyclic Symmetry[END_REF], as shown in Fig. 6. From this discus sion, it turns out that a 3:1 internal resonance is not the only pos sible interaction between modes (0,6) and [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Samaranayake | Resonant Vibrations in Harmonically Excited Weakly Coupled Mechanical Systems with Cyclic Symmetry[END_REF]. Depending on the relative evolution of the frequencies on the backbones of these modes, other n:m resonances with n and m as relatively prime integers can exist. This highlights that NNMs can be internally resonant without necessarily having commensurate linear natural frequencies, a feature that is rarely discussed in the literature.

Another interesting finding is that there is a countable infinity of branches of internally resonant NNMs in this system, similar to what was reported for a 2 DOF system in [START_REF] Kerschen | Nonlinear Normal Modes, Part I: A Useful Framework for the Structural Dynamicist[END_REF]. Figure 16 depicts the same FEP as in Fig. 14, but the algorithm is not stopped after the tongue. Clearly, there is an intricate succession of modal interactions, each one a different realization of an internal resonance between mode (0,6) and another mode of the system.

Additional Mode Bifurcations

According to the degree of detuning for linear mistuned bladed assemblies, double modes with identical frequencies may split into two different modes with distinct natural frequencies [START_REF] Ewins | Structural Dynamics Characteristics of Bladed Assemblies[END_REF]. Because of mode bifurcations, mode splitting may still occur in nonlinear systems with cyclic symmetry. A direct consequence is that a mode pair can bifurcate into two mode pairs, and the number of NNMs exceeds the number of DOFs of the system, as mentioned in Sec. II.B. For illustration, the splitting of one mode of the mode pair [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Wei | Effects of Dry Friction Damping on the Occurrence of Localized Forced Vibrations in Nearly Cyclic Structures[END_REF] is depicted in Fig. 17 (1,14) undergoes a bifurcation during which a new NNM branch is generated. Interestingly, this branch is characterized by modes that are localized to two sectors only. Even more interesting is the bifurcation of mode [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Rosenberg | On Nonlinear Vibrations of Systems with Many Degrees Of Freedom[END_REF], which generates mode shapes that are localized to only one sector. Depending on their stability properties, these NNMs may be excited in practice, giving rise to a potentially harmful motion that must be accounted for.

Mode Stability

In the present study, a stability analysis is performed numerically using the eigenvalues of the monodromy matrix (i.e., the Floquet multipliers), which are a by product of the proposed algorithm. NNM stability is an important concept because it dictates whether or not a mode is physically realizable. For instance, Fig. 20 shows the stability properties of mode [START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Rosenberg | On Nonlinear Vibrations of Systems with Many Degrees Of Freedom[END_REF]. From very low energies to ener gies slightly above 1J, the Floquet multipliers lie on the unit circle. The NNM motions are stable and, hence, physically realizable. From this latter energy, the Floquet multipliers leave the unit circle, and the NNM motions become unstable. This stability change occurs through a bifurcation, which coincides exactly with the generation of the branch of NNMs localized to one sector in Fig. 19.

IV. Conclusions

In this study, the free vibrations of a nonlinear periodic structure with cyclic symmetry were examined using the NNM theory. The NNMs were computed numerically by combining a shooting technique with pseudoarclength continuation. One advantage of this approach is that it provides the capability for analysis of strongly nonlinear regimes of motion.

A very complicated structure of NNMs, including similar and nonsimilar NNMs, nonlocalized and localized NNMs, and bifurcating and internally resonant NNMs, was observed. One important finding of this study is that modal interactions can occur without necessarily having commensurate natural frequencies in the underlying linear system. Furthermore, a countable infinity of such modal interactions were shown to exist in this system.

Further research will investigate how forced resonances of nonlinear periodic structures relate to the topological structure of their NNMs. This will provide useful insight into the forced dynamics, which is of particular practical significance.

  ith disk mass displacement x i = ith blade mass displacement xt = displacement _ xt = velocity xt = acceleration zt = state space vector T = correction to the period during shooting z p0 = correction to the initial conditions during shooting = prescribed relative precision I. Introduction

Fig. 2

 2 Fig.2Frequency-energy plot of system (2) featuring a 3:1 internal resonance between the in-phase and out-of-phase NNMs.

Fig. 3 Fig. 4

 34 Fig. 3 Algorithm for NNM computation.

  m

Fig. 5 Representative

 5 LNMs of the bladed assembly in which the blade and disk masses are shown in black and gray, respectively: a) mode (0,15) and one mode of the mode pair, b) (0,1), c) (1,5), and (d)[START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Peeters | Nonlinear Normal Modes, Part II: Practical Computation Using Numerical Continuation Techniques[END_REF].

Fig. 7 FEP

 7 Fig.7FEP of mode[START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF]0). NNMs represented by bar graphs are inset; they are given in terms of the initial displacements that realize the periodic motion (with zero initial velocities assumed). The blade and disk masses are shown in black and gray, respectively.

Fig. 9 FEP

 9 Fig. 9 FEP of one mode of the mode pair (0,14). NNMs represented by bar graphs are inset; they are given in terms of the initial displacements that realize the periodic motion (with zero initial velocities assumed). The blade and disk masses are shown in black and gray, respectively.

Fig. 10 FEPFig. 11

 1011 Fig.10FEP of one mode of the mode pair[START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Peeters | Nonlinear Normal Modes, Part II: Practical Computation Using Numerical Continuation Techniques[END_REF]. NNMs represented by bar graphs are inset; they are given in terms of the initial displacements that realize the periodic motion (with zero initial velocities assumed). The blade and disk masses are shown in black and gray, respectively.

Fig. 12 FEPFig. 13

 1213 Fig. 12 FEP of the traveling-wave NNM corresponding to a 1:1 internal resonance between the modes of mode pair (1,3) (solid line). For comparison, the dashed line represents the backbone of one standingwave NNM of the mode pair (1,3). NNMs represented by bar graphs are inset.

  Time series corresponding to NNM motions on the tongue of a 3:1 internal resonance (solid line: blade 1; dashed line: disk 10; dotted line: disk 14): a) beginning of the tongue (in the vicinity of the branch of mode (0,6)), b) middle of the tongue, and c) extremity of the tongue (in the vicinity of the branch of mode (1,12)).

Fig. 18

 18 Fig.18Bifurcation of the mode pair[START_REF] Ewins | Modal Testing: Theory, Practice and Application[END_REF][START_REF] Peeters | Nonlinear Normal Modes, Part II: Practical Computation Using Numerical Continuation Techniques[END_REF].

  . Clearly, after the bifurcation, two NNM branches exist and are characterized by different oscillation frequencies and modal curves. Other examples of mode bifurcations are shown in Figs. 18 and 19 for the mode pair (1,14) and mode (1,15), respectively. Mode pair

Fig. 19

 19 Fig. 19 Bifurcation of mode (1,15).

Fig. 20

 20 Fig. 20 Stability analysis of mode (1,15): a) FEP, in which a solid (dashed) line indicates stability (instability); and b) unit circle and Floquet multipliers, which are represented by squares.

Table 1

 1 Natural frequencies of the underlying linear bladed assembly

	Mode	Nodal circles	Nodal diameters Freq., rad=s	Mode	Nodal circles	Nodal diameters	Freq., rad=s
	1	0	0	0.000	31	1	0	2.082
	2,3	0	1	0.183	32,33	1	1	2.084
	4,5	0	2	0.363	34,35	1	2	2.092
	6,7	0	3	0.536	36,37	1	3	2.104
	8,9	0	4	0.700	38,39	1	4	2.123
	10,11	0	5	0.850	40,41	1	5	2.147
	12,13	0	6	0.985	42,43	1	6	2.178
	14,15	0	7	1.103	44,45	1	7	2.215
	16,17	0	8	1.202	46,47	1	8	2.258
	18,19	0	9	1.282	48,49	1	9	2.304
	20,21	0	10	1.346	50,51	1	10	2.350
	22,23	0	11	1.394	52,53	1	11	2.394
	24,25	0	12	1.428	54,55	1	12	2.431
	26,27	0	13	1.452	56,57	1	13	2.460
	28,29	0	14	1.465	58,59	1	14	2.478
	30	0	15	1.470	60	1	15	2.485


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

APPLICATIONS OF GAUSS'S PRINCIPLE OF LEAST CONSTRAINT TO THE NONLINEAR HEAT-TRANSFER PROBLEM

B. VUJANOVIC and B. BACLIC:

Faculty of Technical Sciences, University of Novi Sad, 21000 Novi Sad, Yugoslavia

Abstract-An approximate direct method for solving linear and nonlinear heat conduction problems, based on the Gauss's principle of least constraint is presented. In every particular case, the problem is reduced to the algebraic minimization of a quadratic form with respect to some complex of physical parameters. By principles in use at the present time in the study of heat transfer, are built up in formal analogy with known principles of classical mechanics. This fact confirms the claim of Heisenberg that : " Mechanics could be a model and a basis for all other fields of science ... " ( [START_REF] Heisenberg | Philosophic Problems of Nuclear Sciences[END_REF], p. 89). The most popular principles of mechanics which are advantageously used in many branches of engineering as a starting point for the direct calculation are the D'Alambert's principle of virtual work and Hamilton's principle of stationary action. In addition, when a physical problem is stated in the form of Hamilton's principle it becomes a variational problem in the sense of variational calculus. The direct methods applied on the D'Alambert's principle are generally known as the Galerkin's method, while the numerous direct computations based on the Hamilton's principle are known as Ritz's method, method of partial integration, Rayleigh's method etc.

In this paper we shall study the possibility of applications of Gauss's principle of least constraint to the nonlinear heat transfer. This principle is a true minimum principle in contrast with the two aforementioned principles of D'Alambert and Hamilton, which are generally not the minimum principles. • Contrary to *The variational principle of Hamilton may be occasionally a true minimum principle.

the general impression that the Gauss's principle is of purely theoretical interest in classical (point) mechanics, the purpose of this, and the previous article [START_REF] Vujanovic | The practical use of Gauss's principle of least constraint[END_REF] is, to point out that this principle can lead to considerable practical advantages in nonlinear heat transfer.

GAUSS'S PRINCIPLE

For the sake of clarity we will briefly describe the Gauss principle of ordinary mechanics in a form suitable for applications in heat transfer. where p is the density, k(T) and c(T) are the thermophysical coefficients which are supposed to be the functions of temperature T.

Consider the heat-transfer analog of equation (2) :

z = L [X-YY d v (9)
where Vis the volume which is engaged in the process of heat transfer and

X = div(k •grad T)
y = pc oT [START_REF] Jaeger | Conduction of heat in a solid with a power law of heat transfer at its surface[END_REF] ar are spatial and temporal parts respectively. As in the previous case the following two variational rules for Z =!I (F;-m;a;) 2

(2)

minimizing (9) are possible i=l bX=!-0, bY=O (11)
and suppose that the configuration (r;), velocity (v;) and and forces (F;) of the system are prescribed at timet, i.e. and the remaining "complex"-m;a; is then the only one to be varied, i.e.

(4)

The Gauss's principle of least constraint states that under the conditions (3) and (4) the quantity Z can assume its absolute minimum, which is zero. The proof is simple. If ( -m;a;) represents the actual inertial force, and -m;a;+b(m;a;) represents any other possible inertial force, we have:

As a second possibilityt of achieving the absolute minimum of "constraint" Z defined by (2), we can suppose that the configuration, velocity and acceleration of the system are prescribed at time t and the impressed force F; is then only one to be varied, i.e.

.5r, = 0, bv; = 0, b(m;a;) = 0, bF, =!-0.

(7)

If the system in question is subject to holonomic or nonholonomic constraints the constraint Z is in minimum again but introduction of Lagrange's undetermined multipliers is necessary. We will show now how one can treat the equations of nonlinear heat transfer in the same way as the dynamical systems discussed above. *In the original version of Gauss's principle the expression in the parenthesis of ( 2) is multiplied by the factor 1/m;. tFor more details see [START_REF] Vujanovic | The practical use of Gauss's principle of least constraint[END_REF]. bX = 0, bY=!-0.

(

) 12
This variational principle may be formulated by means of the generalized coordinates instead of the components of the field itself. In numerous problems it is possible to guess that a solution belongs to a family of functions with one or more unknown parameters, which is more or less characteristic feature for all approximate methods. In each particular case we must be able to identify the characteristic complex of parameters which represent X or Y and minimization of Z has to be performed with respect to one of these complexes. The method proposed is very simple because the whole technique for obtaining approximate solution is of the algebraic nature-minimization of a quadratic form with respect to a complex. If some of generalized coordinates are coupled by one or more algebraic •relations the minimization of Z should be performed using Lagrange's undetermined multipliers.

The efficiency of the method can be illustrated by obtaining approximate solution of a simple linear heat-transfer problem.

Consider a one-dimensional thermally insulated semi-infinite body with constant thermo-physical coefficients p = p 0 , c = c0 and k = k0 . The body is initially at the temperature T = 0. At t = 0 the face of the body, located at x = 0, is suddenly brought to a constant temperature T = T0 where rx = k0 /p 0 c0 , together with

T(O, t) = To; T(x, 0) = 0. ( 14
)
Following the concept ofthe penetration depth assume the solution in the form

T = To (1-q~tJ (15)
where q(t) is penetration depth.

Let us solve the problem minimizing the constraint:

Z = lq [aT -rt. azT]z dx Jo at ax 2 (16)
with respect to the temporal part aTjat. From (15)

we have

Y= aT =2To(1-~)~•W at q q a 2 T 2To (17)
axz qz where W=~ (18) q is the "temporal complex" with respect to which we will minimize the expression (16). Substituting (17) into (16) integrating and retaining only the terms with W we have

W 2 q Wrt. Z(W)= - -- 10 q (19)
where a constant multiplicative factor has been omitted. Minimizing (19) with respect to W, i.e.

az(W) =O aw and using ( 18) we find the following differential equation

qq = 5rt.
the solution of which, with respect to the initial condition q(O) = 0 is

q = J(lOrt.t). ( 20
)
Finally let us use the second possibility, minimizing ( 16) with respect to the spatial complex a2T X = rt. ,; yz = 2KT0

aT =2T 0 (1-~)~q, at q q 2
and ( 21) into ( 16), retaining the terms with K, we have after integration

(23)
where a constant multiplicative factor has been omitted.

The equation az(K)jaK = 0, together with (22) will yield the equation qq = 6rt.

i.e.

q = J(12rt.t). ( 24
)
These two results (20) and ( 24) are in good agreement with the exact solution. Besides, they are identical with the approximate solutions obtained by the Galerkin and integral methods respectively.

It should be noted that the Gauss principle of least constraint has been applied on the linear heatconduction problems by Samoilovich in [START_REF] Yu | Gauss's principle in heat conduction theory[END_REF]. However, the author deals with a transformed form of the basic heat conduction equation similar with the Biot's quasivariational method.

It is the purpose of this note to show that the Gauss's principle can be applied directly on the governing heat-conduction equations, and in addition, the nonlinear problems are generally treated in the same way as the linear ones.

As was shown above, the applications of direct methods, using Gauss's principle are simple and straightforward. In the next section several more complex examples will be presented .

EXAMPLES (A) Unsteady two-dimensional nonlinear heat conduction through the prism-like infinite bodies with a given cross section

As the first example we shall study the temperature distribution through the prism-like infinite bodies. The thermal conductivity is supposed to be a linear function of temperature hence the differential equation is of the form aT -r~.!...[o+aT) aTJ-r~. i [o+aT) 0 TJ=o (25) ot ax ox oy ay

where rt. and a are given constants.

Initially, the body is at the uniform temperature, i.e.

T(O , x, y)

= T0 = const. ( 26
)
and the surfaces of the body are maintained at the zero temperature

T(t, X, y)ls = 0 ( 27
)
where Is is the symbol for the external surfaces of the prism. Consequently the problem is to find the approximate solution of (25) in the presence of initial and boundary conditions (26) and (27).

The suitable form of the trial solution will be :

T = tf>(t)B(x, y) (28)
where B(x, y) is a specified function of the position which has to be chosen in accordance with the shape of the cross section of the body. The function B(x, y)

has the following properties: B(x , y)ls = 0, and B(x, y) > 0 for every x andy inside the region bounded by S. Hence, the problem is reduced on finding the unknown function of time tf>(t), using the Gauss's principle.

Let us consider the constraint in the same sense of Gauss:

Z = --rt.-(1+uT)- f x,iy'{oT a l 8TJ xo Yo 0( OX ax -rt.!_ l(l +aT) oT]} 2 dxdy (29)
ay ay

where(x0 ,y0 ) and (x 11 yJ) are to be selected in accordance with the contour in question.

The minimization of Z will be performed with respect to the temporal "complex" ¢. Introducing (28) into (29), integrating, and omitting all terms not containing ¢,we find: (

Z = 11 rP 2 -2r:xl2 4>¢-2rxa !3 4> 2 rP
) 32
The initial condition </J(O) will be determined by minimizing the initial square residual of the form:

J = [ ' [ ' [To-</J(O)B(x,y)J2dxdy
with respect to the arbitrary constant of the general solution of differential equation (32).* As the particular examples we will consider two characteristic shapes of cross section.

(a) The triangle. Consider the triangular cross-section bounded by x = 0, y = 0, x+ y-1 = 0. For this case the trial solution is supposed to be of the form:

( x y)x y T= ¢>(t) 1------ 1 I I I (33) hence ( x y)x y B(x y) = 1------ , I I II (34)
and the limits of integration in (31) are

Xo=Yo=O, Xt=l, Yt=l-x. hence the solution is of the form For this case we will suppose (39)

42To exp(-56r:xt fl 2 ) T(x, y, t) = _ _ Tt _ _ _ _ _ _ _ _ 1 +a 2 ° [1-exp(-56r:xt j/ 2 )] X ( 1 _ T _ n y r (
(40)
and proceeding as in the previous case we find that the temperature distribution is 25T0 exp[ -10 ~:~~2 at]

T(x ,y, t) = [ ( 12 b 2 ) ] 1 +-HaT0 1-exp -10 1 ~2 r:xt X (7-~:)(~-;:). ( 41
)
It is interesting to note that for the linear case a= 0, the corresponding results obtained from (38) and (41) are identical with those obtained by Tsoi in [START_REF] Tsoi | Transient heat conduction in the cylindrical bodies of arbitrary cross-section[END_REF], who used the approximate method based on Laplace transforms. Tsoi reports that for the case of a rectangle his results are in good agreement with the exact solution. Unfortunately for more complicated geometry as triangular for example, the comparison is not possible because the exact analytical solution is not available.

It is reasonable to suppose that the solutions (38) and (41) are of some validity for the moderate range of parameter a.

The corresponding results for both cases (a) and (b) are presented graphically on Figs. 1 and2, where the influence of nonlinearity is presented also.

(B) A melting problem

Consider a semi-infinite solid initially at the melting temperature eP whose surface x = 0 is raised suddenly to the temperature 80 and held there for t ~ 0. We will assume the temperature distribution only in the liquid phase. Such a simplification was proposed by

Goodman [START_REF] Goodman | Application of integral methods to transient nonlinear heat transfer[END_REF] and greatly enhances the use of penetration depth concept in trial solution as far as the latter becomes identical with the location of the melting line ~(t).

Introducing dimensionless temperature T = (8 -8p)/

(80 -Bp), and assuming a liquid of constant thermophysical properties, the governing equation is given by (13), where a modification is made introducing dimensionless time Fo = r:xtj/ 2 where only penetration depth ~(Fo) and the constant A remain for evaluation in accordance with the condition at the interface (45) and the governing heat-conduction equation (13).

Upon substituting (47) into equation (45), one is Jed to the following differential equation for ~ (48) where .u = 2(80 -8p)c0 /L , which when solved for the initial condition ~(0) = 0 gives

~ = 2/J.J(Fo) (49)
where p == .J( -A.u/4).

(50)

It can easily be seen that in this manner all aforementioned conditions are fulfilled by ( 47) and (49). But as trial function (47) fails to satisfy the governing differential equation ( 13), we may proceed, by forming the constraint (16). The temperature profile (47), when substituted into (16), and the integration with respect to X , this time from X= 0 to X= ~(Fo) performed, yields

wz Z = ~(2A 2 -A+2) ~ +!(A 2 -A-2)KW+4(A+1) 2 K 2 ~ (51)
where W = ~ and K = l/~2 are physical components of oT joFo and o 2 T joX 2 respectively. We have two possibilities to obtain the approximate solution of the problem in consideration, i.e. by minimization of the constraint (51): (i) with respect to W , and (ii) with respect to K. In the exact solution J1. is related to the melting constant p by [START_REF] Carslaw | Conduction of Heat in Solids[END_REF] :

Jl. = 2(n:)tp exp(/J 2 )erf(pj. ( 57
)
Since the relation p = P(JJ.) is of most interest in the solution of heat-conduction problems with moving boundaries, in Fig. 4 the percentage error E = lOO(Papprox-Pexac1)/Pexact is plotted against the dimensionless quantity Jl. defined above. J1. = 3 was the largest value of this parameter considered so that values required in practice are covered by this range. The improvement of the accuracy for the approximate solutions obtained here, when compared with those obtained in [START_REF] Goodman | Application of integral methods to transient nonlinear heat transfer[END_REF] using the heat-balance integral, and in equation (5.28) of [START_REF] Eckert | Analysis of Heat and M ass Transfer[END_REF] by the help of a rather incorrect treatment of the problem, with a trial solution of the same form, is self-evident. We also conclude that the optimization with respect to temporal change of temperature field yields, in this case, much better results than the same with respect to spatial change of temperature.

(C) Semi-i'!finite body with an arbitrary heat flux input

As the last example, consider the case of the transient heat-conduction problem in the semi-infinite solid with constant thermal properties whose initial temperature is zero. In order to involve both linear and nonlinear boundary conditions the assumption will be made that the heat flux at the surface x = 0 is an arbitrary function of surface temperature T. = T(O, t) and time. Hence, the problem can be stated mathematically as follows. The governing equation is (13), with initial condition T(x, O) = 0 and boundary condition at

x = 0 aT ko -+ f(T. , t) = 0. ox (58) (59)
Following references [START_REF] Goodman | Application of integral methods to transient nonlinear heat transfer[END_REF] and [START_REF] Vujanovic | Heat transfer with nonlinear boundary conditions via a variational principle[END_REF] let us assume the cubic temperature profile in the form:

T = T.(t) (I-q~tJ (60)
where q = q(t) is the penetration depth and both q and T. are unknown functions of time.

From the boundary condition (59) it follows that penetration depth and surface temperature are dependent

3ko T. = qf (61)
which indicates that we are going to deal with a constrained optimization problem in the appliance of the above exposed direct method to find T. and q. In the sense of this fact, we will use the Lagrange's multiplier technique, when dealing with the minimization with respect to physical components of temporal change of temperature.

Let us start by substituting the trial function, equation ( 60 and performing the indicated integration, one obtains t/q T. f.q 3 T.2q2

Z = -+ --+ n -- 7 7 q • 2 2
12 T. T. 9 T. . 2 2 T.

---,CJ. ----;er: -q+l

C/. - q q2 q3 (63)
where the dot represents differentiation with respect to time. The next step is to recognize the physical components of temporal and spatial change of temperature in this expression so that the minimization procedures may be performed straightforwardly. Since aT t. CJ.-2 = 6 2 (q-x)er: (66) ax q has only one complex:

T, K=a-(

q In order to demonstrate two possible solutions we will solve the problem minimizing corresponding Gauss's constraint with respect to temporal and spatial change of the assumed temperature distribution respectively. (i) Minimization with respect to temporal change. Omitting last term in (63) as irrelevant since it does not contain components W1 and W2 defined by (65), and using (61) whence 3k0 f.-4! -q(f'f.+ /) = 0 (68) where/= of /o t andf' = of/8 4 , we have

3ko T. W? 4 W1 Wz 4fWl z = 7f + 7 + 35ko 4af af2 -! -W1 -! -2 Wz ko ko ( 3ko 4Wd' +J. 3k o W~-fWz- f 3kof4f)
where). is an unknown Lagrange's multiplier. Minimization with respect to oT /ot, i.e. -k5 kf, .

Assuming that, in second case, f = fo • tN where N is an integer, equation (77) becomes t(f.)2 + H N42 = * af~2 t2N +' (79) ko and reduces to the first case for N = 0. The solution to equation ( 79) is

fo ( 49at )t N 4 = k 0 45N +39 t • (80)
Upon separation of variables the differential equation for the third case (78) yields fr, 784/ 2 -66T/jf'

+ 184 3 !' 2 T. _ 7at Jo / 3 (7f -3T.f') d ' -kij • ( 81 )
Let us solve this equation for the case of power low heat transfer at the body surface into a medium at zero temperature, i.e.

f(T,) = H(T.-To)'"

where Hand mare constants. Introducing 11 = 1 + T. / To equation (81) becomes Upon elimination J. from equations (71) and (72) one is led to 45kij T. f.+ 11fko 4q-49af 2 -15k5 T. 2 f. j'

-6k0 4 2 4[' + 2lafT.f' = 0 (73)
which when combined with equations ( 68) and (61) gives the first order differential equation for the temperature at the surface (74)

The initial condition for the problem is T.(O) = 0. There are three cases in which equation ( 74) can be integrated analytically:

1. f = const. In this case equation (74) reduces to

(75) (83)
The method of handling (82) applies also to the case in which the surface flux is a sum of powers of the surface and environment temperatures, and thus to heat convection from bounding surface into a fluid at temperature T 1 , for which the boundary condition is (linear case) f=h(4-T 1 ), and also to the blackbody radiation into a medium at temperature

Tm :f = a(T, 4 -T;). The discussion of the results given above will take place after solving the same problems using the second possibility, viz. minimizing the constraint (63) with respect to the component of spatial change of temperature.

(

2. f = f(t) 3. f = j(T,) f T. = 1.127-(a.t)t ko 2 r. t. _ 6 T. 2 ! = 80 a.P s s 7 f 03" k6 ibf T• l4T,f-6Tlf' _ext ! 3 dT,-2• o ko (90) (91)
where ' 7 = 1 -T. / T 1 .

Let us now proceed with the discussion and the comparison of the results. One thing is obvious-the closed form solutions obtained by minimization of Gauss's constraint with respect to the component of spatial change of temperature distribution are more attractive since they are of a rather less complicated form. To evaluate the accuracy of the results we are going to compare all of them with the exact solutions for surface temperature.

In the case of constant heat flux the exact solution is known to be T. = 1.128(! /k0 )(rxt)+, which when compared with equations (76) and (91) indicates the improvement of the result obtained by the help of minimization with respect to K. But when compared with the result l.l5(f /k0 )(at)t of application of the integral method [START_REF] Goodman | Application of integral methods to transient nonlinear heat transfer[END_REF] to the same problem, both results are better.

To compare the results for the case when heat flux is the power function of time the corresponding ratios for the approximate solutions obtained by variational [START_REF] Vujanovic | Heat transfer with nonlinear boundary conditions via a variational principle[END_REF] and integral [START_REF] Goodman | Application of integral methods to transient nonlinear heat transfer[END_REF] method are also listed.

(f = fo • tN) we calculate T r(N+ 3 ) r(N) = .-approx. = 7 • R(N) (100)
Comparison of the approximate solutions with the numerical ones [START_REF] Jaeger | Conduction of heat in a solid with a power law of heat transfer at its surface[END_REF], for the case when the heat flux is a power low function of the surface temperature, is evident from Fig. 5. The dots represent approximation equations ( 96), (97) and (98) for m = 4, 2 a nd 5/4 respectively. If one superimposes the results of Goodman [START_REF] Goodman | Application of integral methods to transient nonlinear heat transfer[END_REF], or Vujanovic and Strauss [START_REF] Vujanovic | Heat transfer with nonlinear boundary conditions via a variational principle[END_REF] (available only for m = 4), or those obtained by minimization with respect to o Tjot in this paper, it may be concluded that all of them are less accurate.

-Exact When the solution of equation ( 99) is compared with the exact one, and the percentage error E = 100('7approx. -'lexactl/'lexact is calculated it turns out that E = 3%, for (T./ TJ)approx = 0.8, which reduces to 1% for (T./ T.)approx. = 0.5, while for (T./Tf)approx. ~ 0.3, E ~ 1%. Doing the same with equation ( 86): E = 6.4% at (T./ T,)approx = 0.8; E = 1% at (T./Tf)approx. = 0.4 and less than 0.5% for (T./TJ)approx. < ~.3.

Thus, in all cases the solutions here obtained differ only slightly from the exact ones. Further, the important conclusion is that the result obtained by minimization of Gauss's constraint with respect to the component of spatial change of temperature distribution are always more accurate, although of rather less complicated form.

DISCUSSION

In conclusion several remarks would be of interest:

1. The method we have presented is quite useful in the search for approximate solutions of linear and nonlinear heat-conduction problems. Usually, in the case of nonlinear analysis, it pertains to a high accuracy when compared with the other approximate methods.

2. A common feature of all approximate methods in heat transfer, including the method presented, is that the solution of a problem should be selected to some extent a priori. The choice of the form of a solution which contains some parameters that should be determined by the help of the Gauss's method depends upon all the information available from empirical, experimental, intuitive, etc. data.

Our aim has been to demonstrate two possibilities

for obtaining approximate solutions from the same Gauss's constraint. The question as to which of the two possible approximate solutions should be taken is to be decided by considering the mean square residual [START_REF] Finlayson | The Method of Weighted Residuals and Variational Principles[END_REF] (p. 388) in the form Z=f''f (X-Y) 2 dVdt.

(101)

'o Y
Naturally, the value of Z defined by ( 101) is equal to zero for the exact solution. Hence, the better solution is that one for which Z has the smallest value. It should be pointed out that the same criteria can be employed for evaluating the accuracy of any approximate solution obtained by some other approximate method.

4. The extension of the method presented in this paper on some more elaborated mathematical schemes as for example the method of finite elements is possible and will be reported elsewhere.

5. This method can be extended also in a straightforward way to handle the numerous problems arising in transport phenomena.

  Consider a free dynamical system of n particles, subject to impressed forces F;(i = I, 2, . . . , n). If m; are masses, r; position vectors and a; = d 2 r;/dt 2 accelerations, the differential equations of motion are (I) Let us introduce the quantity* Consider the nonlinear equation of heat conduction t' T div(k • grad T)-pc ---:. ; -= 0.

  n bZ =±I {[F;-m;a;+b(m;a;j]2-(F;-m;a;) 2 } i=l =!I [6(m;a;)] 2 + I (F;-m;a;)b(m;a;) (5) i=l i=l and the last term vanishes in virtue of (I). Thus bZ>O unless b(m;a;) is zero.

  ..,-=0 o<P yields the differential equation /1 rP-l2r:x¢>-MlJ¢> 2 = 0.

  the constant of integration. The equation oJ -= 0 oCz *This standard procedure for finding the initial condition of </>(t) was applied by many authors previously. See for example [ 4]. will yield 2 56 C 2 = -3[2a-421 2 To (37)

  38) (b) Rectangular cross-section. Let the boundaries of a rectangle be defined by x= ±l y =±b.

FIG. 2 .

 2 FIG. I. Temperature distribution in the triangular crosssection atFo = rxt/4/ 2 = 0.08 for linear (u = 0) and nonlinear (a= 0.5) case, (a) temperature distribution in the center line (x = y), (b) isothermals T = 0.2 and 0.4.

  (i) Minimization with respect to spatial change of temperature K. Condition oZjoK = 0 together with K = 11e yields the differential equation 2 ~~ = 24(A + 1) .

2 -A - 2 2~~=FIG. 3 .

 223 FIG. 3. Negative roots of equation (56).

8 FIG. 4 .

 84 FIG. 4. Percentual error in determination of melting constant by different approximate methods.

  ), into z = fq[(oT)2 -2CJ. oT o2T + CJ.2(o2T)2] dx (62)

It

  IS obvious that temporal change of temperature distribution possesses two natural components (65) On the other hand spatial change of temperature a 2 T r.

2 + 2 . 7 : 1 :

 2271 105J.k5f(1-T/) = 0 (71) 5kij T. f.+2fko T.q-7af 2 -35J.k5f = 0. (72) or 4 =f... (~at)t = 1.121£ (ar)t ko ko Iff = f(t), equation (74) becomes and 2 + 11 Ts 2 / = 49 ajl ' ' n f 3"9" k5 , 3. Whenf = f(4), equation (74) yields . 42 t. , 6 43 t.rz 244-HT f +TJ-1 -2- -~ af2-ti r4f'T.

Final 2 1 -

 21 ao = 3/7 a1 = -40/49 a2 = 107/252 as= -275/31104 a6 = -1375/ 248832 a1 = -6875/ 1492992 b = -34375/17915904 (84) a3 = -11/360 a4 = -55/3456 c = 376528644/ 13168189440 (b) form= To2H2 -rxt = -[i--H'7+ih2-fi'73] nln'l (86) where ' 7 = 1 -T./T 1 since the boundary condition is taken of the formf = h(T 1 -T,).

  the second case we choose f = fo • t N to get the solution of (92) T. _ fo ( 80a.t )t N s -ko 72N + 63 t • (94) Assuming againf(T,) = Trf • H • Y/m where Y/ = 1 + T./ To we obtain the solution for the third case in the form 9 f~(7-3m)'1 2 +(6m-7)1]-3m y 0 2m -zH2 case of heat convection boundary conditionf(T,) = h(T 1 -T.) from equation (93) one obtains

FIG. 5 .

 5 FIG. 5. Surface temperature of the semi-infinite solid with surface flux the mth power of the temperature. The dots represent approximation equations (96), (97) and (98).

  the help of several concrete examples the efficiency and accuracy of this new method is demonstrated.

		NOMENCLATURE	fJ,	temperature;
	a,	acceleration;	fJP,	temperature of phase change.
	b, c, E, F, Fo, H,	characteristic length; specific heat; percentage error; force; Fourier number rxt/1 2 ; coefficient;	Subscripts 0, refers to initial values; f, refers to the fluid; m, refers to the surrounding medium ; S, refers to the surface of the body.
	h,	convective heat-transfer coefficient;
	K,	complex representing the component of
		spatial change of temperature distribution;
	k,	thermal conductivity ;
	L,	latent heat of melting ;
	/,	characteristic length;
	m,	mass;
	q,	penetration depth;
	r,	position vector;
	t,	time;
	T,	nondimensional temperature;
	v,	velocity;
	W,	complex representing the component of
		temporal change of temperature distribution ;
	x,	length coordinate ;
	X ,	nondimensionallength xj l;
		spatial change of temperature distribution
		(10);
	y,	length coordinate ;
	Y,	nondimensionallength yjb;
		temporal change of temperature distribution
		(10);
	Z,	Gauss's constraint.
	Greek symbols
	rx,	thermal diffusivity ko/poco ;
	{3,	melting constant;
	A.,	Lagrange multiplier ;
	p,,

= 2(8o -lJp)co/L; e. location of melting line; p, density; a, constant specifying the dependence of k on T; Stefan-Boltzmann constant in example C; I. INTRODUCTION THE USE of variational methods for obtaining solutions of transport phenomena has become an indispensable tool in the domain of nonlinear conduction. Many publications dealing with variational principles of various kinds have appeared in the technical literature in the past 20 years. Roughly speaking, all variational

  ii) Minimization with respect to spatial change. This time omission of terms not containing component K

	defined by (67) reduces constraint (63) to the following
	expression of relevance:
	12 Z = -,KT,--s-K-q+12 -. • 9 T. . K 2 q q	(87)
	Optimization condition
	az aK=O	(88)
	together with (67) will yield, upon using (61) and (68),
	the differential equation
	. 63T.T.-27f(T.f'+ )=40 k6 T.2 . f afz	(89)
	determining surface temperature.
	If we consider, as above, the cases for which equation
	(89) is integrable, the following reductions are evident:
	1. f = const.
	2T. t. = 80 af 2 s s 03" k5
	with the solution