Yang Wang

Shanshan Fu

Enrico Zio

Di Zhang

Xinping P Yan

Modeling epistemic uncertainty in

whether they are published or not. The documents may come

INTRODUCTIONS

1.1 Resilience in water transportation systems Safety has been a major concern in water transportation due to the potentially hazardous and dangerous environment. The state-of-the-art technology has provided advanced equipment for waterborne navigation and shore-side surveillance, yet the inherent complexity and the uncertainty in water transportation still make it a safety-critical system for the researchers and practitioners.

Risk is strongly featured in water transportation, and the study on the maritime risk analysis or assessment has attracted sustained work for years. For this reason, risk analysis has attracted sustained work for years in the field of water transportation. Recently, resilience is receiving growing attention [START_REF] Zio | Challenges in the vulnerability and risk analysis of critical infrastructures[END_REF]. The resilience study analyzes the adverse situation with multiple facets rather than only from the angle of pre-event prediction. Strategic planning, pre-event preparedness and post-event recovery are all encompassed in the resilience study [START_REF] Zobel | Characterizing multi-event disaster resilience[END_REF]. Emphasis is placed on the capability of the system to withstand disruptions and the actions taken to hedge against the damages or losses by absorption, adaption and restoration.

This shift of paradigm in safety research is expected to bring forth profound influence upon the domain-specific research. Nevertheless, the study of resilience encounters similar difficulties as follows.

The first challenge comes from the setup of the problem, i.e., how to formulate resilience for water transportation, which are human-intensive, loosely coupled system, with autonomous bodies (the ships) distributedly moving a large area. Some recent research work has been carried out on resilience of water transportation. [START_REF] John | A risk assessment approach to improve the resilience of a seaport system using Bayesian networks[END_REF] studied the resilience of a seaport system by using Bayesian networks, with the CPT (conditional probability table) being estimated by a fuzzy analytical hierarchy process; [START_REF] Praetorius | Modelling vessel traffic service to understand resilience in everyday operations[END_REF] investigated the resilience aspect of the VTS (Vessel Traffic Service) operation in maritime administration by using FRAM (Functional Resonance Analysis Method).

The second challenge is the quantification of resilience. Results from theory exist [START_REF] Ferrario | Goal Tree Success Tree-Dynamic Master Logic Diagram and Monte Carlo simulation for the safety and resilience assessment of a multistate system of systems[END_REF] in the literature, yet the engineering application to the transportation domain calls for a practical, explainable and suitable quantitative model such that the process and behavior of the system is clearly reflected [START_REF] Gao | Universal resilience patterns in complex networks[END_REF]. Currently, the resilience quantification is conceptually prescribed in terms of the system performance curve, but to compute such curve in compliance with multiple and interrelated factors, with uncertainties and evolution over time, is quite challenging.

The third challenge relates to the pre-event assessment. Resilience analysis can be applied as both pre-event and post-event assessment in which ex ante knowledge is indispensable because there usually do not exist absolute laws with which the inducing factors can conform [START_REF] Kaplan | Fitting hierarchical holographic modeling into the theory of scenario structuring and a resulting refinement of the quantitative definition of risk[END_REF]. Prior knowledge from past experience, and subjective judgement from analysts should be highly accommodated in the evaluation model [START_REF] Creaco | Comparison between entropy and resilience as indirect measures of reliability in the framework of water distribution network design[END_REF].

To meet the above challenges, this paper proposes a network of factor model according to the function, component and structure analysis of water transportation systems. The contribution aims to build a way to formulate knowledge for obtaining resilience curve and its quantitative assessment.

THE BACKGROUND

The resilience and modeling of system evolution

The ex ante resilience assessment should be conducted in terms of more specific disruption. Figure 1 illustrates the conceptual resilience curve proposed by [START_REF] Henry | Generic metrics and quantitative approaches for system resilience as a function of time[END_REF]. A resilient system should show a curve with shallow drop and a steep rise after the point of minimum system resilience is reached. However, this quantitative description only conveys a conceptual basis for the analysis, and in order to obtain a more rigorous quantitative form, various aspects should be addressed, which including the system dynamics and modeling the uncertainty of system behavior under disruption.

Possibly relevant methodologies for quantification are DET (Dynamic Event Tree) [START_REF] Hofer | Dynamic event trees for probabil istic safety analysis[END_REF], and the DBN (Dynamic Bayesian Networks) [START_REF] Murphy | Dynamic Bayesian Networks: Representation, Inference and Learning[END_REF]. DET is utilized to analyze the stochastic events that occur during accidental evolution from an initial adverse event [START_REF] Karanki | A dynamic event tree informed approach to probabilistic accident sequence modeling: Dynamics and variabilities in medium LOCA[END_REF]. DBN is a time-sliced Bayesian network, in which the conditional probability among the nodes is dependent upon the correlation of the factors not only in-slice but also cross-slice. Both DET and DBN models have been extended from their original static versions to model the discrete event [START_REF] Mercurio | Identification and classification of dynamic event tree scenarios via possibilistic clustering: application to a steam generator tube rupture event[END_REF]. However, limitations exist for the two models. DET does not fully account for the inter-restrictions among factors, and DBN does not explicitly indicate the time scale. Besides, DBNs assume that the values of some of the nodes can be observed, which is not the case for an ex ante study.

The rest of the paper is organized as follows. Section 3 discusses the methodology and proposes a model in response to the requirements; section 4 introduces some of the properties of the model; section 5 conducts numerical study by Monte Carlo simulation and discusses the inclusion of knowledge; section 6 summarizes the work and suggests the future research.

MODELING THE CYCLE OF DISRUPTION IN WATER TRANSPORTATION

The generic model of incident disposal

Human-intensiveness in water transportation system implies that the interactions between the human and the software/hardware are closely interweaved, and that the decisions and actions from some key individuals can largely influence the overall situation. Moreover, the main decision/action entities in a context of water transportation accident are the ship and the shore-side administration. Figure 2 illustrates the simplified process of parallel decision and action dynamics between the two entities, when confronting an incident. The humans in the system mainly include the crews (at the ship-side) and the administrative personnel (at the shore-side). Decision objective refers to the aims of decision, and action plan refers to the implementation options. Although the diagram above is structured with some resemblance to Bayesian networks by a DAG (Directed Acyclic Graph), the process would hardly be straightforwardly modeled by it, as traditional Bayesian network model does not support the duration of a status. In this regard, the following two lay-er structure is proposed as a generic model, as is de-picted in Figure 3. With the two layer structure, the process of the accident development is divided into two planes: the decision plane (DP) and the evolution plane (EP). The decision plane is described by a series of Bayesian networks, which illustrate the logical dependence. The evolution plane is used to depict the evolving process by sequential state transitions.

The idea behind is to create a plane that can accommodate the system state transition over time, while the nodes in the original plane are regarded as transient nodes that do not reflect time factors. In Figure 3, the node S1 has logically direct interlace with T0. Obviously, there should be a time elapse between the time slice 0 and time slice 1, hence the two time slices are bridged by a continuum of state transitions. It is worth noting that, each time slice is a snapshot of the system, and that the time elapse is represented by the number of transitions in EP.

Quantification

The conventional BN is formulated as a DAG, in which each node stands for a variable. These variables represent the factors that make the building blocks of the whole reasoning map. We use a triplet

 

,,  G N E P to denote the BN, where

  12 , , ..., X X X  N N denotes the set of nodes,   , , , i j i j i j X X X X X X    EN
is the set of directed edges, and the directed pair , ij XX indicates that is the parent of j X , or j X is the child of i X . A node may have none or multi i X ple parents/children.

Denoted by  

j X  =   , i i j X X X  E the parent set of j X
, then is the CPT set of all the nodes. For a node j X , each of its parent nodes   j X  can have impact on j X , and this dependence is depicted by the CPT. For a BN, the foundational property is given by the factoring of the joint probability distribution as follows:

      12 1 , ,..., ii i p X X X p X X     N N
By placing a set of BNs in a sliced tandem, a general DBN is formed. In DBN, the directed edges can connect two nodes within a slice, or across two slices. Moreover, DBN allows a node at tth time slice to be conditionally dependent not only on its parents at the same time slice (intra-slice) but also on its parents and its own states at previous time slices (inter-slice) [START_REF] Khakzad | Application of dynamic Bayesian network to risk analysis of domino effects in chemical infrastructures[END_REF].

The proposed two layer model is based on the DBN with some modifications. For the decision plane, the inter-slice dependence is allowed only in two consecutive slices. This assumption does not lose practical value because the epistemic cognitive capability of an analyst will decay sharply when applied to longer term reasoning. For the evolution plane, the sequential nodes can be regarded as an inter-slice edge connecting two BNs.

Referring to Figure 3, we use the italic font to indicate the variables that are hosted in the corresponding nodes. In the conventional BN, each node stands for exactly one variable. However, in real problems, if a combination of variables can be treated as one, it will help to make the model more concise. Theoretically, there does exist an equivalent variable in the product space for the combination of the sub-variables, hence, we also use one italic capitalized letter to denote the variable when necessary and with little ambiguity.

IMPLEMENTATION OF THE MODEL

System status and scoring

The essential information that a resilience curve aims to convey is the quantified system status and its changes over time. In this sense, the status of a system is associated with a scored value. Denote by

 

12 , , , m s s s  the state set of the system in question,

      i i i p X X X   PN and by   [0,1] i Rs  1 im
 the scoring of each state. It can be assumed without loss of generality that be the one-step Markov state transition matrix over  . We assume that the state transitions take place at uniform step, so that the time gap between two consecutive transitions is a constant  . In other words, the number of transitions determines the evolution time. According to the Markov property, if the distribution of the system states is

        1 2 3 10 m R s R s R s R s       ,
() r  at the rth transition, then ( 1) r   = () rM  .
As is illustrated in Figure 3, after the action plans, when Ei and Fi take effect, the system enters the evolving process until another external event occurs. The action plans are modeled by the state transition matrix. Furthermore, the initial state and the initial event within the time slice on the DP are also the causal factors that lead to the beginning state of the evolution on the EP. The third key parameter is how long the evolution will last before renewed evolution starts.

All necessary parameters are obtained from the data in the knowledge, and are grouped in the node Ti in the model of Figure 3.

The nodes Ti and Si+1 are bridged by a series of interim states in EP. This network setup reflects the dependence relationship that the node Ti involves. In fact, the node Ti is a comprehensive variable that in fact consists of three components: M introduced in in the previous section. As also discussed in the previous section, def i T can be represented as the number of state transitions for the inter-slice evolution.

For the input dependence, the i T is jointly determined by the prior state ( i S ), the nature of the event ( 0 I for the initial and i U for the follow-up), and the actions taken to dispose the event from the ship-side ( i F ) and from the shore-side ( i E ).

THE TEST CASE

Application background and the states set

In this section, the proposed two-layer model is applied to the case of not under command (NUC) ships. NUC is a kind of accident with non-negligible occurring-rate in the middle and lower reaches of Yangtze River, and it will give raise to secondary accident if it is not tactically handled. The main causes of NUC include the failures of main engine/rudder, or black out. For the purpose of simplicity, we define six elements for the state set of the system, as is described in Table 1 below.

The envisioned event set 

Let  be the event set that makes the basic sample space for I0 and Ui . Table 2 lists the elements of the event set for the problem of interest. These events are considered within the scope of epistemic knowledge, and the sequence and inter-dependence are described by a series of Conditional Probability Tables related to actions and system states. A bridge appears at the downstream between the NUC ship and the nearest anchorage.

The pre-accident event.

u

There is dense vessel traffic between the NUC ship and the nearest anchorage.

u

The anchorage is at the downstream shortly ahead. The NUC ship anchors.

The accident is partially under control.

u

First/More Tug(s) or maritime cruiser is approaching.

u

The whole accident is completely settled. (There may be different extent of damage and loss.)

The final event.

Table 3. The decision objectives of the NUC incident.

Di (ship-side decision objectives) Ci (shore-side decision objectives)

Variable outcomes

Description

Variable outcomes

Description d1

Immediately dispose the incident at the cost of potential hindrance to the circumambient traffic or other ships.

c1

Protecting the hydraulic structures (e.g. the bridge) is of higher priority.

d2

Tolerate a certain amount of sacrifice to stabilize the situation quickly. c2

Protecting the NUC ship and the ships passing by is of higher priority.

d3

Take adventures to seek the most elegant settlement. d4

Scrambled choice is taken. Send one/more tug(s) to the incident site to facilitate further operation.

f2

Try to beach or anchor in outer limit of the channel. e2

Send one/more maritime cruiser(s) to the incident site.

f3

Try to traverse crowded waters or a bridge to find the anchorage. e3

Regulate the sailing restriction of all ships.

f4

Try to immediate anchoring in the channel. e4

Regulate the one-way traffic organization.

5.3

The decision objectives of the ship-side and the shore-side

The decision objectives reflect the preferences of the decision makers. In the resilience assessment context, the analysts do not know exactly which actual preference the decision makers will adopt, and the judgements made by the analysts come from the experience data they know about. In the same way, the quantitative formulation of the epistemic knowledge is presented by CPTs. Table 3 lists the possible outcomes of the variables Di and the Ci.

5.4

The action plans of the ship-side and the shoreside

The action plans belong to the set of control options to be taken by the decision-makers, and they are the implementation tactics in accordance with the given decision objectives. Although confronting the same scenario, the action plan chosen by the two sides(the ship-side and the shore-side) may not necessarily be the most compatible. This compatibility is usually due to the asymmetry of the scenario perception or the ambiguity in the ship-shore communication. Table 4 lists the action plan variables and the possible variable outcomes.

Building the conditional probability tables with knowledge

The CPT of Ci and Di

Considering that the case is focused on the evaluation of NUC ship accident starting from normal state, we set the prior probability of the variable S0 and I0 to be the value for the root node S0 and I0, such that   0 1,0,...,0 S  , I0=i0. With this regard, a sample CPT of Ci and Di is given below in Table 5 and Table 6. Table 5. The CPT of C0 and D0.

The CPT of Ei and Fi

The variables Ei and Fi depend upon Ci and Di, respectively. Table 7 depicts the dependence in terms of probabilities.

The component of Ti and the CPT

As has been introduced previously, the node Ti hosts three variables, namely def i T , inst i T and forward i T . In the context of the current problem, the number of state transition is 1 10 def i T  . In practice, if we assume that the time scale of transition is a constant, for instance, one minute, def i T defines the evolution length for the implementation of the decision. inst i T is a 1⨯m vector to indicate the instant distribution of system state after the event occurs, and forward i T is a m⨯m stochastic matrix to denote the state transition. Due to space limit, the details of these matrices are omitted. Table 8 shows a section of the CPT. Table 8. A section of the CPT for the component of Ti.

i E i F 0 / i IU i S def i T inst i T forword i T 1 i U  Probability e1 f1 i0   1,

Monte Carlo simulation results and analysis

Using the above setup, a series of Monte Carlo simulations are conducted to gain a family of resilience curves given the NUC as the initial event. Figure 5 depicts the sampling of one round of simulation.

The curves convey the rough resilience behavior under the epistemic estimation provided by the analysts. To characterize the resilience in a more straightforward way, we employ the resilience measurement for each resilience curve, as follows:

1 0 10 () t t f t dt tt    R (1)
This definition of resilience measure is very similar to what is given by [START_REF] Bruneau | A framework to quantitative ly assess and enhance the science the seismic resilience of communities[END_REF]. The main idea of (1) is to find the proportion of non-loss part during the disruption and recovery process. t0 stands for the starting point of disruptive event, and t1 stands for the time point of complete recovery. Figure 6 shows the interpretation of the measure based upon the resilience curve. Obviously, greater R implies a stronger capability to withstand the disruption. Figure 7 displays the resilience measures of the multiple simulation tests where the simulation is repeated for 1600 times to account for the complexity of the resilience process. From the result, it can be seen that the samples of full-cycle disruption handling show their resilience measure largely within the interval [0.7, 0.9], while there may be extremely negative cases with resilience as low as 0.4 or even less.

CONCLUSION AND FUTURE WORK

In this paper, we present a method for the epistemic assessment of system resilience by using a revised model of a dynamic Bayesian Network and a discrete Markov process. The key features of the work are outlined below. (i) Resilience focuses on the process and dynamics of system evolution under negative accident, and this is the very point of difference between the resilience study and the traditional risk study. The process-based view calls for a temporal change of system state instead of the cause-effect analysis, and this proves to be a challenge in modeling. In this paper, the revised DBN is used to represent the timesliced system performance. (ii) When we deal with subjective assessment, uncertainty is a major concern. There are many mathematical tools to treat the issue; this paper employs the Bayesian network as the reasoning engine due to the human-intensive, loosely coupled and autonomous nature of the application domain. As the knowledge of individuals is usually confined to local causal relationship, the DBN can mechanically extend the reasoning to longer spans of time and broader interlace of factors. (iii) There are some limitations in the presented work. For the utmost, the knowledge is finally represented by CPTs and state transition matrices. However, when the states (or outcomes) of the variables increase, the size of CPTs will grow exponentially and usually, we have to resort to some rulebased reasoning or merging to reduce the size of CPTs. This will hinder the practical assessment setup, and improvement to this problem will be of interests for future work.

(iv) Despite the aforementioned limitations, the presented work provides a useful method of ex ante assessment for system resilience with respect to a given disruptive scenario. The quantitative results can facilitate choices of system design. For instance, a possible improvement is to enhance the capability maturity of the ship-side and shore-side operators, such that the joint tactics can be reached with less time and more compatibility.

Figure 1 .

 1 Figure 1. The conceptual resilience curve (Henry et al. 2012).

Figure 2 .

 2 Figure 2. Accident management model of the ship-shore coordination.

Figure 3 .

 3 Figure 3. An intuitive two layer structure of accident management.

  meaning the sequential states of i s are of descending order in terms of the system score, with being normal and m s being complete failure. In this manner, the i S ( i ℕ∪{0}) of node Si in DP are probability vector variables over [0,1] m , where the sum of the vector components is 1. If one is concerned with the resilience curve from the very beginning, i.e., the normal state, he can set the initial distribution for the root node as

Figure 4 .

 4 Figure 4. A partial enlargement of the node Ti's dependence.

4 u

 4 There are relatively open waters around the NUC ship.

Fi

Figure 5 .

 5 Figure 5. The resilience curves generated by one round of simulations.

Figure 6 .

 6 Figure 6. The resilience measure values based on the calculation of the resilience curve.

Figure 7 .

 7 Figure 7. The resilience measure of multiple tests.

Table 1 .

 1 The state set of the NUC incident.

	State

1 s Normal state. All the post-accident operation is completed, and the loss and damage have been ascertained. 1 2 s If no action is taken, accident is unavoidable.

A major part of search and rescue (SaR) operation is completed.

Table 2 .

 2 The event set of the NUC accident.

	Event	Description	Remark
	i0	Ship not under command (NUC) is detected.	The initial event
	1 u

Table 4 .

 4 The option plans of the NUC incident.

Table 6 .

 6 A section of the CPT of Ci and Di (i⩾1).

	S	i	Ui	Di	Prob	Ci Prob
				d1	0.2	c1	0.7
	-	u1	d2 d3	0.3 0.4	c2	0.3
				d4	0.1
				d1	0.25	c1	0.3
	-	u2	d2 d3	0.25 0.4	c2	0.7
				d4	0.1

… …

Table 7 .

 7 The CPT of Ei and Fi.

	Ci	Ei	e1	e2	e3	e4
		c1	0.4	0.3	0.1	0.2
		c2	0.2	0.2	0.4	0.2
	Di	Fi	f1	f2	f3	f4
		d1	0.1	0.4	0.1	0.4
		d2	0.2	0.3	0.2	0.3
		d3	0.2	0.3	0.4	0.1
		d4	0.1	0.3	0.3	0.3

ACKNOWLEDGEMENT The research is supported by Key Project in the National Science & Technology Pillar Program (Grant No. 2015BAG20B05), and by National Natural Science Foundation of China (Grant No. 51479158).


Yang Wang

Enrico Zio

Shanshan Fu

Di Zhang

Xinping P Yan

Some reflections on pre-and post-accident analysis for water transport: A case study of the Eastern Star accident

Some reflections on pre-and post-accident analysis for

water transport: A case study of the Eastern Star accident

INTRODUCTION

Water transport is a traditional industry that can date back to the most ancient human activities. Despite the fast transition of the world's socio-technical being, waterway remains the essential means of transportation, accounting for over 90% of the world's trade [START_REF] Rodrigue | The geography of transport systems[END_REF]. During its long course of development, water transport has been always considered potentially hazardous and dangerous. Regardless of the truth or falseness of such consideration, records of numerous maritime accidents, such as Costa Concordia, Sewol and Eastern Star disasters, definitely show that waterway transportation is an industry intrinsically associated with risk.

In view of this, there has been a large volume of academic research work devoted to safety analysis for water transport [START_REF] Goerlandt | A probabilistic model for accidental cargo oil outflow from product tankers in a shipship collision[END_REF]. Apart from the risk analysis studies [START_REF] Li | An overview of maritime waterway quantitative risk assessment models[END_REF]), the study on accident analysis, as the other side of the coin, is equally receiving ardent interest from the researchers. Apparently, the study of risk and accident analysis in the field of maritime transportation mainly draws its ideas, concepts, models and approaches from the general study of risk and accident analysis [START_REF] Yang | Maritime safety analysis in retrospect[END_REF][START_REF] Goerlandt | A framework for risk analysis of maritime transportation systems: A case study for oil spill from tankers in a ship-ship collision[END_REF], which is within the domain of safety research. In another sense, the study of general risk and accident analysis can find its application and validation by the data and cases in the maritime transportation (Zhang et al. 2014a, b, Montewka et al. 2014).

In the literature, ad hoc research on pre-accident and post-accident analysis has appeared extensively [START_REF] Wang | A spatial-temporal forensic analysis for inland-water ship collisions using AIS data[END_REF]. For the majority of the preaccident study, the risk-based perspectives and methods dominate the perception [START_REF] Goerlandt | A framework for risk analysis of maritime transportation systems: A case study for oil spill from tankers in a ship-ship collision[END_REF]; for the typical post-accident study, statistics or data of the past accidents are exploited to find the pattern of accidents, and these refined findings can in turn cast risk indications for accident prevention.

The risk-related study encompasses a broad spectrum of academic perceptions. In its most naï ve form, risk can be interpreted as the product of likelihood and the consequences of adverse events (Rausand 2011). Nevertheless, there have been different understanding about the rigid definition of risk [START_REF] Aven | On some recent definitions and analysis frameworks for risk, vulnerability, and resilience[END_REF][START_REF] Haimes | Responses to Terje Aven's paper: on some recent definitions and analysis frameworks for risk, vulnerability, and resilience[END_REF]. A sophisticated definition may contain uncertainty, knowledge, or system state vectors etc [START_REF] Haimes | On the complex definition of risk: a systems-based approach[END_REF].

Accident models are the mainstream problems as well as methodologies in the ex-post study of accidents. The main categories of accident models include the sequential model, the epidemiological model and the systemic model [START_REF] Hollnagel | Cognitive systems engineering: new wine in new bottles[END_REF]. There have been numerous variants developed to interpret the causes, effects and reasoning of the accident process. AcciMap, systemtheoretic accident model and process (STAMP) and functional resonance accident model (FRAM) ABSTRACT: In the literature, research on pre-accident and post-accident analysis has been reported extensively. The former falls into the category of risk analysis to identify the potential hazards and the consequences, so that decisions can be made to take appropriate preventive and mitigation measures. The latter amounts to accident investigation for learning lessons and finding solutions. This paper uses the catastrophe of Cruiser Eastern Star in 2015 as a case study to analyze some typical accident models, including the Swiss Cheese Model, AcciMap, Stamp, FRAM, etc. By comparison, pros and cons of each perspective are identified. Overall common challenges for all methods are pointed out and a reflection is offered with respect to the added value of complementing the analysis under a resilience perspective, to enhance preparedness for unknown adverse events. emerge as the prominent methods to analyze accidents [START_REF] Underwood | Systems thinking, the Swiss Cheese Model and accident analysis: a comparative systemic analysis of the Grayrigg train derailment using the ATSB, AcciMap and STAMP models[END_REF]. It can be envisioned that the study of accident makes a stepping stone to the study of safety issue in water transport.

This paper conducts a case study on the catastrophe of Cruiser Eastern Star in 2015 to acquire a more comprehensive understanding about the domain-specific accident analysis. Typical accident models, including the Swiss Cheese Model (SCM), AcciMap, STAMP, FRAM and risk-based perspective are utilized to yield multi-angled views. By comparing these accident models, pros and cons of each perspective are identified, and challenges for them are pointed out. Based on these results, reflections are made with respect to the added value of complementing the analysis under a resilience perspective.

The remaining of the paper is organized as follows. Section 2 introduces the maritime accidents used for the case study. Section 3 analyzes the five relevant accident analysis methods. Section 4 discusses the strengths and weakness of the five methods. Section 5 concludes the research findings.

THE EASTERN STAR ACCIDENT

While the disasters in water transportation can always arouse the sensation of people, some of the stakeholders may show great scrupulousness to release the informative accident reports. Research on the foundering of Costa Concordia in 2012 began to appear shortly after the accident (Schröder-Hinrichs & Hollnagel 2012). Although the loss caused by Costa is much less than that of the TITANIC in 1912, it is somehow compared with the latter for introspection. Two years after the accident, the capsizing of Sewol in 2014 once again outreaches the scope of the experiential safety knowledge [START_REF] Hyungju | Assessment of accident theories for major accidents focusing on the MV SEWOL disaster: Similarities, differences, and discussion for a combined approach[END_REF].

Each occurrence of a distress will widen human's knowledge on safety/hazard, and correspondingly, the knowledge base of the researchers and the practitioners will be monotone increasing. The ideal hypothesis is that, with the ever growing knowledge base of accidents, the ability of foreseeing or coping with the accidents should be strengthened consequently, and that the tendency of accident occurrence will decline asymptotically.

Unfortunately, the actual situation is by far worse than the idealized envision. In 2015, the inland passenger ship, Eastern Star, once more interrogates the theory and practice of risk or resilience in water transport. Eastern Star is a cruise ship who operates the itineraries along the Yangtze River, featuring the Three Gorges sightseeing. The length of the ship is 76.5 meters, the beam 11 meters and the capacity 534 passengers.

The ship departed from Nanjing at the time 1300 hours (GMT+8, the same below) on 28th May, and the destination was Chongqing in the upper reaches of the Yangtze River. At the time 1930 hours on 1st June, 2015, the ship entered the JianLi waters, with heading nearly to the north. At this time of the trip, the regional waterway was being struck by a severe rainstorm. The ship rushed into the central area of the strong convection at the time 2119 hours, and was hit by the fierce northwest linear gale to the port. The ship swerved to the left during the time 2120 to 2121 hours, with heading jerked from 23° to 342°, aiming to counteract the gale. Approximately at the time 2132 hours, the ship capsized to the starboard within 1 minute, and was drifting downstream afterwards. The last AIS signal was broadcasted at the time 2131 hours at 29°42'N, 112°55'E.

The sinking of the Eastern Star cost 442 lives, making it the most fatal maritime accident in China since 1949. Only 12 persons on board survived from this accident, including the master and chief engineer officers.

ACCIDENT ANALYSIS VIA DIFFERENT MODLES

The barriers and defences

From the perspective of SCM [START_REF] Reason | The contribution of latent human failures to the breakdown of complex systems[END_REF]), several layers of barriers are placed to hedge an accident. As for the case of Eastern Star, the four defensive layers are described in Figure 1.  Layer 4: The maneuvering group underestimate the severity of the weather condition, and choose a venturesome strategy to break through, thus leaving the whole ship at stake.

Socio-technical system and AcciMap

To bring an accident in question into the sociotechnical context therein will offer a wider panorama for the analytic treatment. In his seminal paper, [START_REF] Rasmussen | Risk management in a dynamic society: a modelling problem[END_REF] introduces the dynamical boundary of the safe behavior. Inspired by this idea, the analysis of Eastern Star is illustrated in Figure 2.

In the case of waterway transportation, the accident analysis should specially consider the multiple actors within the scenrios whose actions can have joint impact on the evolution of accident.

Figure 2 presents an AcciMap-style analysis of the Eastern Star [START_REF] Rasmussen | Proactive risk management in a dynamic society[END_REF]. The original AcciMap model has six system levels, but in this analysis, the top level is omitted because the relevance is less direct and the omission will cause no loss of the essence. By AcciMap model, the accident analysis is placed in a panorama of sociotechnical context.

System-theoretic accident model and process

With the ideas of STAMP, systems are viewed as interrelated components that are kept in a state of dynamic equilibrium by feedback loops of information and control [START_REF] Leveson | A new accident model for engineering safer systems[END_REF][START_REF] Licu | EUROCONTROL-Systemic Occurrence Analysis Methodology (SOAM)-A "Reason" based organizational methodology for analyzing incidents and accidents[END_REF]. The three bedrocks of STAMP are safety constraints, hierarchical levels of control, and process loops.

Through the lens of STAMP, the major casual factor of the Eastern Stars accident is the constraint failure, which does not prevent Eastern Star from entering the wrong place at the wrong time. And this failure is mainly due to the block of some control loops that are supposed to take effects over time.

Figure 3 depicts the STAMP-style control structure of the accident scenario. The cruise company does not employ any surveillance infrastructure to trace the ship, which means no link with the ship; MSA does not obtain the adequate information about the meteorological conditions, and does not maintain a peer-to-peer supervision with the ships involved in the foreseen dangerous situation; the communication among the ship is not functioning well due to noise or disturbance of the storm.

From Figure 3, we can see that the Eastern Start bears few control loops from other components of the system. This suggests that the ship is subject to no constraint or control in the extreme environment. And finally, the unfortunate encounter of the traversing squall line becomes the last straw.

For simplicity, the further qualitative process delineation in line with the STAMP template is omitted.

Functional resonance accident model

Using the language of FRAM, the system is viewed through a functional abstraction rather than structural decomposition. The interaction among the functional entities governs the behaviors the system. In this regard, the context of Eastern Star is expressed by three functional entities as can be demonstrated with the three hexagons in Figure 4. The relations among the functional entities depicted by the arrows reflect the normal behaviors and performance variability in a routine sense.

In FRAM, function resonance refers to the high performance variability of a functional entity, which is deemed as a sign of a potential accident. Setting the damping factors is the practical way to curb the function resonance, and the accident in the end. However, the FRAM shows certain drawbacks in analyzing accident due to extreme events, as in the case of Eastern Star, where the extreme event usually cannot be manifestly reflected by a functional entity in the FRAM diagram.

The risk analysis in hindsight

Risk analysis is normally conducted a priori to prevent or mitigate risks [START_REF] Zio | An introduction to the basics of reliability and risk analysis[END_REF][START_REF] Zio | Challenges in the vulnerability and risk analysis of critical infrastructures[END_REF]. As for the case of Eastern Star, we explain to what extent the risk analysis can yield meaningful results. This can help us understand the gap between the current theoretical presentation and the accident-prone activities in the real world. Here, again, the attempt is to set out from the practical and engineering angle.

We follow the widely accepted conceptual risk description (C, Q, K) [START_REF] Aven | A new perspective on how to understand, assess and manage risk and the unforeseen[END_REF], with the three elements standing for the event/ consequence (C), quantitative measure of uncertainty (Q) and background knowledge (K), respectively. We assume that the following analysis is carried out as if we were at the time prior to the Eastern Star accident.

A broad spectrum of hazard events can be involved and investigated in the navigation of Eastern Star, at the time she starts the journey. For instance, windstorm, collision, contact/grounding, fire/explosion are main hazard sources [START_REF] Mazaheri | Modeling the risk of ship grounding-a literature review from a risk management perspective[END_REF]. The variety of root causes makes it difficult with selecting or exhausting the "relevant" events. Even focusing on a certain type of hazard event, the variance in time and space will still cause a considerable variety of events.

The measure of uncertainty is largely a technical problem, in which there are alternative tools, such as probability theory, Bayesian Network, D-S theory and possibility theory [START_REF] Zio | Industrial disasters: extreme events, extremely rare. Some reflections on the treatment of uncertainties in the assessment of the associated risks[END_REF]. These tools are backed by their own theories with concrete mathematical foundations. It is worth noting that the measure of uncertainty is the direct reflection of the background knowledge. In most quantitative risk assessment (QRA), the measure of uncertainty and the knowledge share the same theoretical kernel, once the knowledge is represented with a specific formalized model. Although the Q and K have very close affiliation, the risk description regards them as separated input to better denote the various knowledge from different sources.

In the case of Eastern Star, if we could apply the knowledge-featured risk description to the ship prior to her last journey, we shall have to confront the additional challenges as follows:

 The actual accident is unexampled in several aspects with respect to the accident data base during the past 15 years, including the location (Damazhou waterway), the root cause (entering the center of a storm), the ship type (the cruise ship) and the master's competence (highly experienced master). If the knowledge relies too much on the frequency or probability of similar accidents based on the historical dataset, the risk analysis may give false negative results.

 The knowledge about the Eastern Star journey can be classified into open knowledge and private knowledge. The former represents information that is publicly available, such as the meteorological and hydrological parameters; the latter represents the information that is only available for a few insiders or professionals, such as the foreseen squall line, the past rebuilding record and the SCNWP(Stability Criterion Numeral of Wind Pressure) of the ship. This disparity in knowledge will result in disparity in risk assessment from different analysts.

 For the cost-effective purpose, the risk analysis would be conducted with less sufficient knowledge and for more abstract objective, e.g., at a higher scale. For instance, it is more likely for MSA to analyze the risk of navigation of a specific range of waterway during a specific period, instead of the risk of a specific ship like Eastern Star, and this leads to the lacks of discriminatory information for individual ships.

In this manner, the risk analysis is often reduced to a less informative form, i.e., the early warning message.

DISCUSSION

Based on the above conventional treatments of the Eastern Star analysis, two questions are of special interest for the discussion:

 Is the Eastern Star case truly an inevitable accident? Or otherwise, where is the least modification that can prevent the misfortune if everything could be repeated?

 Apart from the latter-wit, how can we distill something from these conventional methods that can make the system more immune to accidents, even though we are blind about when, where and, how the next accident will occur?

The first question is equivalent to finding the cause of accident. The root cause can be ascribed to the irresistible natural force. However, it is more pertinent to categorize this accident into human factor and organizational factor. It is due to human factor that the crew fails to make the right decision to prevent the accident; it is due to organizational factor that the navigational information and warning from the shore-side administration is inadequate and untimely. Other factors, such as the degrading of the ship's resistance against wind pressure, are not the dominating factors.

The second question is rather ambitious, and there may be no off-the-shelf answers. Table 1 summarizes the limitations of the aforementioned treatments for accident analysis. It has to be acknowledged that most of the existing analytical approaches are post-hoc study, meaning that we have already known the incident and try to make a backward analysis. However, it is purely a paradox that we could know something about the future accident and take precautions.

As has been mentioned previously, one of the key characteristics of accidents is the unpredictability. The unpredictability can relate to the time, the space/location and the causation etc. Although human's knowledge is broadened with each occurrence of accident, and we endeavor to take preventive measures based on the accumulating knowledge, we still can see no extinction of accidents.

The reason can be interpreted by Figure 5, which illustrates the space of causal element. With each occurrence of accident, the number of hollow dot increases and the circled area around them expands.

However, with the ever-growing precision, compactness and complexity of large systems, the number of unknown causal elements (the solid dots) may be rising at even faster speed over time. In other words, what we know about the casual elements of future accidents is much less than the potential casual elements which we don't know or are out of our anticipation. Returning to the Eastern Star case, despite the limitations listed in Table 1, we are in the position of making the most out of the existing theories and methods. To this end, we resort to a resilience view to find a partial solution to the second question [START_REF] Hosseini | Modeling infrastructure resilience using bayesian networks: a case study of inland waterway ports[END_REF]. Regarding resilience, the most salient shift from traditional risk analysis (or assessment, management) lies in the holistic thread for adverse event responding [START_REF] Baraldi | A combined Monte Carlo and possibilistic approach to uncertainty propagation in event tree analysis[END_REF]. The main limit for risk study, as we can think of, is the relatively short run-length of consideration span along the time dimension for the adverse possibility. The reaction and recovery cycle of the adverse situation disposition, is less covered in the conventional risk study while it is the cycle that the resilient study places emphasis on. It is difficult to define and identify the abnormal performance variability for the functional entity 3

The spatial and temporal interdependence cannot be manifestly described

(C,Q,K) Risk Analysis 1
Disruptive events can hardly be exhaustively enumerated 2

Subjective knowledge differs significantly among multiple analysts 3

Statistics on historical data may lose effect for the foresight of extreme events such as catastrophes due to their low frequencies 4

Post-disruption situation is usually beyond the boundary of the current risk study

CONCLUSION

The diversity of accident models provide multiple angles to gain insights into the cause path and mechanism of the accident, whereas the risk analysis is supposed to help stakeholders evaluate the potential hazard before accident.

The innate difficulty for these methodologies is also very straightforward: what will these methodology bring forth for the decision-makers when no accident occurs?

Obviously none of the above methodology can act as a fortune-teller by foretelling when, where and which type of accident will occur. Nevertheless, these methodologies can show their values by enhancing the consciousness of safety, including the following aspects:

(1) Understanding the structure, function and interplay of the components to know the trajectory of any potential accident in terms of time, space and logic dependence.

(2) Being aware of the potential hazards according to the dynamic status of the system;

(3)Paying more attention to the human factors and the organizational factors in human-intensive systems such as waterway transportation.

With the growing complexity of the modern socio-technical development, it is ever difficult to exhaust the potential casual factors to guarantee an accident-free system. This paper uses the catastrophe of Cruiser Eastern Star in 2015 as a case study to analyze five typical accident models, including the SCM, AcciMap, STAMP, FRAM and risk-based perspective for accidents analysis. These methods generally fall into the categories of pre-and post-accident analysis. Based on the comparison, it is shown that all the conventional methods demonstrate some limitations in dealing with extreme events. For the safety of a complex system like the water transport, the utmost difficulty to tackle is the unknown or the unpredictable hazard, and this makes a paradox for researchers. With this understanding, this paper finally tends to favor a more practical paradigm of thoughts, i.e., resilience, as a promising perspective. More work is anticipated to incorporate the traditional point of view into the new mindset, and this will be the direction of our future research.

Figure 1 .

 1 Figure 1. A SCM-style analysis. The realization of the accident means that the holes in the barriers do exist, and cause the failures of the defensive layers, as follows:  Layer 1: The cruise company tacitly consents to adventurous sailing in treacherous environment at night.  Layer 2: The early warning messages dispatched by the maritime safety administration (MSA) come out too late, and there is a lack of decisive and strong intervention in the vessel traffic at the time of hazardous weather condition.

Figure 2 .

 2 Figure 2. An AcciMap-style analysis.  Layer 3: Poor information about the hazardous weather condition is shared over the VHF communication among the ships near the area. Six ships are involved in the storm simultaneously, but the reciprocal communications and behavior is scarce.

Figure 3 .

 3 Figure 3. A STAMP-style safety control structure analysis.

Figure 4 .

 4 Figure 4. A FRAM-style function entity diagram.

Figure 5 .

 5 Figure 5. A space of causal element.

Table 1 .

 1 Comparison of the aforementioned treatments for accident analysis.

	Method

ACKNOWLEDGEMENT The research is supported by Key Project in the National Science & Technology Pillar Program (Grant No. 2015BAG20B05), and by National Natural Science Foundation of China (Grant No. 51479158).


Yu V Mikhlin

B I Morgunov

Normal Vibrations in Near-Conservative Self-Excited and Viscoelastic Nonlinear Systems

Keywords: Self-excited nonlinear systems, nonlinear normal modes (NNMs), viscoelastic nonlinear systems, power series, iterations

A perturbation methodology and power series are utilized to the analysis of nonlinear normal vibration modes in broad classes of finite-dimensional self-excited nonlinear systems close to conservative systems taking into account similar nonlinear normal modes. The analytical construction is presented for some concrete systems. Namely, two linearly connected Van der Pol oscillators with nonlinear elastic characteristics and a simplest twodegrees-of-freedom nonlinear model of plate vibrations in a gas flow are considered.

Periodical quasinormal solutions of integro-differential equations corresponding to viscoelastic mechanical systems are constructed using a convergent iteration process. One assumes that conservative systems appropriate for the dominant elastic interactions admit similar nonlinear normal modes.

Introduction

Lyapunov proved that nonlinear finite-dimensional systems with an analytical first integral allow a one-parameter family of periodic solutions which tend towards linear normal vibration modes as amplitudes tend to zero. Natural frequencies of corresponding linearized systems must not be integrally related [START_REF] Lyapunov | The General Problem of the Stability of Notion[END_REF]. The Lyapunov's solutions possess all the properties of linear normal modes [START_REF] Mikhlin | Normal vibrations of a general class of conservative oscillators[END_REF].

Nonlinear normal vibrations modes (NNMs) are a generalization of the normal (principal) vibrations of linear systems. In the normal mode, a finite-dimensional system behaves like a conservative one in having a single degree of freedom, and all position coordinates can be analytically parametrized by any one of them.

Rosenberg, in a series of papers [START_REF] Rosenberg | The normal modes of nonlinear n-degrees-of-freedom systems[END_REF][START_REF] Rosenberg | Nonlinear vibrations of systems with many degrees of freedom[END_REF][START_REF] Rosenberg | Nonsimilar normal mode vibrations of nonlinear systems having two degrees of freedom[END_REF] provided analytical methods for computing NNMs. He defined them as 'vibrations in unison' and introduced a broad class of essentially nonlinear conservative systems allowing for NNMs with rectilinear trajectories in a configuration space ('similar' NNMs).

For example, 'homogeneous systems' whose potential is an even homogeneous function of the variables belong to such a class. It is interesting to note that the number of the NNMs may exceed the number of degrees of freedom of the oscillator. In general, the NNMs trajectories are curvilinear instead of straight lines in linear systems.

For some particular cases, curvilinear trajectories were analyzed by Rosenberg and Kuo [START_REF] Rosenberg | Nonsimilar normal mode vibrations of nonlinear systems having two degrees of freedom[END_REF] and by Rand [START_REF] Rand | A higher order approximation for nonlinear normal modes in two degrees of freedom systems[END_REF]. The power series method was proposed by Manevich and Mikhlin [START_REF] Manevich | Periodic solutions close to rectilinear normal vibration modes[END_REF] for the construction of the curvilinear trajectories of NNMs.

Later, new results concerning NNMs of conservative systems were performed by several authors, including Vedenova et al. [START_REF] Vedenova | Normal oscillations of a string with concentrated masses on non-linearly elastic support[END_REF], Vakakis [START_REF] Vakakis | Non-similar normal oscillations in a strongly non-linear discrete system[END_REF], Vakakis and Rand [START_REF] Vakakis | Normal modes and global dynamics of a two-degree-of freedom non-linear system. I. Low energies. II. High energies[END_REF], Shaw and Pierre [START_REF] Shaw | Nonlinear normal modes and invariant manifolds[END_REF][START_REF] Shaw | Normal modes for nonlinear vibratory systems[END_REF], Nayfeh and Nayfeh [START_REF] Nayfeh | On nonlinear modes of continuous systems[END_REF], and Mikhlin [START_REF] Mikhlin | Matching of local expansions in the theory of nonlinear vibrations[END_REF].

In [START_REF] Vedenova | Normal oscillations of a string with concentrated masses on non-linearly elastic support[END_REF], nonlinear normal mode localization is studied in non linear systems with impact nonlinearities. In [START_REF] Vakakis | Non-similar normal oscillations in a strongly non-linear discrete system[END_REF], for an oscillator with weak coupling stiffness, both localized and nonlocalized modes are detected using an asymptotic methodology. In [START_REF] Vakakis | Normal modes and global dynamics of a two-degree-of freedom non-linear system. I. Low energies. II. High energies[END_REF], the NNMs and global dynamics of nonlinear systems are analyzed by means of Poincaré maps. In [START_REF] Shaw | Nonlinear normal modes and invariant manifolds[END_REF][START_REF] Shaw | Normal modes for nonlinear vibratory systems[END_REF], the authors reformulated the concept of NNMs for a general class of nonlinear discrete oscillators. The analysis is based on the computation of invariant manifolds of motion on which the NNM oscillations take place. In [START_REF] Nayfeh | On nonlinear modes of continuous systems[END_REF], a complex invariant manifold formulation is proposed. In [START_REF] Mikhlin | Matching of local expansions in the theory of nonlinear vibrations[END_REF], Padé approximations are used for an analysis of NNMs with large amplitudes.

Note that publications on the generalization of NNMs to nonconservative systems are not numerous.

It is well known that forced vibrations in nonlinear systems with one degree of freedom and small periodic disturbances are similar to those in an undisturbed conservative system (in the resonance region).

It should be noted that periodic solutions in nonautonomous systems close to Lyapunov systems were thoroughly investigated by Malkin [START_REF] Malkin | Certain Problems of the Theory of Nonlinear Vibrations[END_REF].

The Rauscher method was used for analyses of normal vibrations in nonautonomous systems by Kinney and Rosenberg [START_REF] Kinney | On the steady state vibrations of nonlinear systems with many degrees of freedom[END_REF]. To solve the problem, one assumes that the external periodic force is proportional to some chosen positional coordinate raised to such a power that the resultant autonomous system is homogeneous.

NNMs in general finite-dimensional nonautonomous systems close to conservative systems allowing similar NNMs were considered by Mikhlin [START_REF] Mikhlin | Resonance modes of near-conservative nonlinear systems[END_REF]. Rauscher's ideas and the power-series method for trajectories in a configurational space are used in the construction of resonance solutions.

Some other works on NNMs in nonautonomous systems, have been published by Yang and Rosenberg [START_REF] Yang | On forced vibrations of a particle in the plane[END_REF], Szemplinska-Stupnicka [START_REF] Szemplinska-Stupnicka | The modified single mode method in the investigation of the resonant vibration of nonlinear system[END_REF], and Vakakis and Caughey [START_REF] Vakakis | A theorem on the exact nonsimilar steady state motions of a nonlinear oscillator[END_REF].

Of particular note is the work by Rand and Holmes [START_REF] Vakakis | Normal Modes and Localization in Nonlinear Systems[END_REF], where periodic solutions of two weakly coupled Van der Pol oscillators are analyzed. Note that NNMs in a model of two strongly coupled Van der Pol oscillators are considered in this work.

Basic and new results on NNMs are presented in the book by Vakakis et al. [START_REF] Manevich | The Method of Normal Oscillation for Essentially Nonlinear Systems[END_REF] which describes quantitative and qualitative analyses of NNMs in conservative and nonautonomous systems, including localized modes, an analysis of stability, a generalization, and an analysis of NNMs in distributed systems.

Here, the perturbation methodology is utilized for the analysis of NNMs in broad classes of finite-dimensional self-excited nonlinear systems close to conservative systems allowing normal vibrations with rectilinear trajectories in a configuration space (Section 1). The analytical construction of NNMs is presented in Section 2 for some concrete systems. The first of them represents two linearly connected Van der Pol oscillators with nonlinear elastic characteristics. Another example refers to the problem of plate vibrations in the flow of gas. A simple, two-degrees-of-freedom nonlinear model is considered. In Section 3, periodical quasinormal solutions of integro-differential equations corresponding to viscoelastic mechanical systems are constructed. One assumes that elastic forces are dominant and conservative systems appropriate for the elastic interactions allow similar NNMs which are selected as generative solutions. A convergent iteration process is used here.

Normal Vibrations in Near-Conservative, Self-Excited Nonlinear Systems

EQUATIONS DESCRIBING THE TRAJECTORIES OF NNMS

We will demonstrate that the approaches for analysis of NNMs in conservative systems are applicable to near-conservative autonomous systems with small self-excited perturbations. We will also consider modes when all positional coordinates of the finite-dimensional, self-excited system are linked. In relation to these modes, the system behaves like a single-degree-offreedom one. The periodic solutions close to NNMs in the generating conservative system could be called nonlinear normal modes (NNMs) of the self-excited nonlinear system.

Consider the following near-conservative system:

ẍi + x i (x 1 , x 2 , . . . , x n ) + εg i (x 1 , ẋ1 , . . . , x n , ẋn ) = 0, (2.1)
where ε is a small parameter and the functions g i may be nonlinear with respect to ẋi . The functions f i and are assumed to be analytical in x i , ẋi ; a potential energy is subject to some limitations which will be shown later. The system may involve friction of any physical nature, such as viscous, dry, or turbulent. The conditions that ensure the occurrence of selfexcited vibrations (limiting cycles) are discussed later. Assume that a conservative system (c = 0) permits similar NNMs, i.e., normal modes with rectilinear trajectories in a configuration space: x i0 = k i x i0 (i = 2, . . . , n; k i are constants). Any terms may be regarded as a perturbation.

Consider one of these NNMs as a generating one. Rotating the coordinate axes so that the new x-axis is directed along the rectilinear trajectory, one obtains a generating solution in the new coordinates as

x i0 (i = 2, 3, . . . , n), x i0 ≡ x = x(t).
(2.2)

Select some positional variable x ≡ x i .

Let us find a solution that all phase coordinates are defined as single-valued and analytical functions of x:

x i = x i (x, ε), ẋ = ẋ(x, ε), ẋi = ẋi (x, ε). (2.3)
Introducing a new independent variable x instead of t, one obtains from (2.1) the equations describing the trajectories of NNMs:

x i ( ẋ(x)) 2 + x i [-x (x, x 2 (x), . . . , x n (x)) -εg 1 (x, ẋ(x), x 2 (x), ẋ2 (x), . . .)] + x i (x, x 2 (x), . . . , x n (x)) + εg i (x, ẋ(x), x 2 (x), ẋ2 (x), . . .) = 0 (i = 2, 3, . . . , n). (2.4)
Here and henceforth, a prime denotes differentiation with respect to x; it is clear that ẋi (x, ε)

= x i (x, ε) ẋ.
Select the trajectory return points where all velocities are equal to zero. When x = X j (j = 1, 2) and ẋ = 0, we obtain the additional conditions from (2.4):

x i [x (x, x 2 (x), . . . , x n (x))εg i (x, ẋ(x), x 2 (x), ẋ2 (x) . . .)] + x i (x, x 2 (x), . . . , x n (x)) + εg i (x, ẋ(x), x 2 (x), ẋ2 (x) . . .) x=X j = 0 (i = 2, 3, . . . , n; j = 1, 2).

(2.5)

Here (X, x 2 (X), . . . , X n (X)) being the trajectory return points. Equations (2.4) in combination with conditions (2.5) allow a univalent determination of x i (x). Note that equations similar in form and additional boundary conditions were previously considered in the case of conservative systems [2-4, 22, 23].

In the zero approximation

(ε = 0) ẋ2 0 = 2[h 0 -(x, 0, . . . , 0)],
where h 0 is the energy of the nonperturbed conservative system. We shall restrict our consideration of (2.4) to the first approximation in ε:

2x i1 [h 0 -(x, 0, . . . , 0)] + x i1 [-x (x, 0, . . . , 0)] + n k=2 x k1 x i x k (x, 0, . . . , 0) + g i (x, ẋ0 (x), 0, . . . , 0) = 0 (i = 2, 3, . . . , n).
(2.6)

In this approximation, the boundary conditions are obtained from (2.4) at ẋ0 = 0, x = X j (j = 1, 2):

x i1 [-x (x, 0, . . . , 0)] + n k=2
x k1 x i x k (x, 0, . . . , 0)

+ g(x, ẋ(x), 0, . . . , 0)

x=X 1,2 = 0 (i = 2, 3, . . . , n). (2.7)
The power series method may be applied to find x i1 (x).

THE POWER SERIES METHOD

It is assumed, for example, that the generating conservative system is homogeneous, i.e. is an even homogeneous function of the power of r + 1 in all the variables (r may take the values r = 1, 3, 5, . . .). Since the unperturbed system is homogeneous, the matrix B = x i x k (x, 0, . . . , 0) may be written as B = b ik X r-1 . Note that owing to the conservative nature of the system, b ik = b ki , and the symmetric matrix b ik is reduced to a diagonal form by nondegenerate linear transformation of the positional coordinates [START_REF] Poincaré | Les Méthodes Nouvelles de la Mécanique Celeste[END_REF]. The transformation can be chosen in order to preserve a generating solution in the form (2.1). Therefore, without loss of generality, one can assume that in (2.6) and (2.7) the function x i x k (x, 0, . . . , 0) = 0 for i = k. Hence, the set (2.6),together with the boundary conditions (2.7) is 'split' in the variables x i1 .

Having substituted the series

x i1 = ∞ j =0 a ij x j (2.8)
into (2.6), one finds that the coefficients a ij are interrelated by the following infinite set of recurrent relationships (it is written out coefficients at x r+j ):

2h(r + j + 2)(r + j + 1)a i,r+j +2 -j (j + 1)2 (1, 0, . . . , 0)a i,j +1 -(j + 1) x (1, 0, . . . , 0)a i,j +1 + x i x i (1, 0, . . . , 0)a i,j +1 = φ i,r+j (j = 0, 1, 2, . . .), (2.9)
where φ i,r+j denotes the corresponding terms in the Taylor-series expansions of the functions g i . The recurrent relationships can be used to express all coefficients a ij in terms of the leading coefficients a i0 and a i1 (i = 2, 3, . . . , n). Substituting the series x i1 (x) into the boundary conditions (2.7), and in view of the recurrent relationships (2.9), one obtains the equations governing a i0 , a i1 :

n i=2 (R 0i a i0 + R 1i a i1 ) = R 2i (i = 2, 3, . . . , n). (2.10)
The determinant of the system may be represented as products of an infinite number of factors [START_REF] Mikhlin | Normal vibrations of a general class of conservative oscillators[END_REF][START_REF] Manevich | The Method of Normal Oscillation for Essentially Nonlinear Systems[END_REF][START_REF] Bellman | Introduction to Matrix Analysts[END_REF]] (2.11) where δ j i are Kronecker's delta, p = 0, 1, 2, 3, . . . . When the generating system is linear, conditions (2.11) mean an absence of internal resonances in the generating system.

K p = |q ij |, q ij = δ j i [p(p -1)2 (1, 0, . . . , 0) + p x (1, 0, . . . , 0) -x i x i (1, 0, . . . , 0)],
The procedure presented here is similar to that applied in a conservative case [START_REF] Mikhlin | Normal vibrations of a general class of conservative oscillators[END_REF][START_REF] Manevich | The Method of Normal Oscillation for Essentially Nonlinear Systems[END_REF][START_REF] Bellman | Introduction to Matrix Analysts[END_REF], and it is possible to make use of analogous arguments. In particular, a convergence of the series (2.8) can be proved by L'Hospital's rule.

It follows from Poincaré's theorem on the small parameter series expansion [START_REF] Ilushin | Foundations of the Mathematical Theory of the Thermoviscouselasticity[END_REF] that there is a value of ε 0 > 0 such that, for all |ε| < ε 0 , the series in powers of ε of the functions (2.3) converge in the domain; as ε = 0 this solution becomes the trivial generating solution x j 0 = 0 (j = 2, 3, . . . , n) (see the details in [START_REF] Mikhlin | Normal vibrations of a general class of conservative oscillators[END_REF][START_REF] Manevich | The Method of Normal Oscillation for Essentially Nonlinear Systems[END_REF][START_REF] Bellman | Introduction to Matrix Analysts[END_REF]).

Note that the trajectories x j (x) can be derived not only in terms of power series in x, but also by the method of iteration [START_REF] Mikhlin | Normal vibrations of a general class of conservative oscillators[END_REF][START_REF] Bellman | Introduction to Matrix Analysts[END_REF].

SINGLE-DEGREE-OF-FREEDOM AUTONOMOUS SYSTEM AND A POTENTIALITY CONDITION

After having obtained a smooth NNM trajectory, the problem of finding a periodic solution reduces to the integration of the following single-degree-of-freedom autonomous system:

ẍ + x (x, x 2 (x), . . . , x n (x)) + εg 1 ( ẋ, x(x), x 2 (x), ẋ2 (x), . . .) = 0.
The solution x(t) is determined as an inversion of the following quadrature:

t + φ = 1 √ 2 x x(0) dξ (h -(ξ, x 2 (ξ ), . . .) -εg 1 (ξ, ẋ(ξ ), x 2 (ξ ), . . .)) 1/2 , (2.12)
where h is the integration constant which has the meaning of an energy constant of the conditional single-degree-of-freedom system realizable along the NNM analytical trajectory; φ is a phase of the solution which can be taken to be equal to zero because the original system (2.1) is autonomous.

Next, the equation of phase trajectories is of the form ( ẋ = v)

dv dx = - x (x, x 2 (x), . . . , x n (x)) + εg 1 (x, v, x 2 (x), x 2 (x)v, . . .) v .
In the general case, one has to employ approximate methods to find all phase trajectories and due allowance is to be made for the small values of ε. However, this equation allows exact integration if g 1 is a linear function with respect to w = v 2 .

At the analytical closed phase trajectory of the limiting cycle, the system's behavior is similar to that of a conservative system with one degree of freedom. Therefore, the condition that the work of all forces over the period is equal to zero should hold (the condition can be called a potentiality condition):

g 1 (x, ẋ(x), x 2 (x), ẋ2 (x), . . .) dx = 0 or T 0 g 1 (x, ẋ(x), x 2 (x), ẋ2 (x), . . .) ẋ dt = 0.
(2.13)

The integral is taken over the period of vibrations. Equation (2.8) is used to find the value of the energy h 0 (and, accordingly, the amplitude values of X j ) for the generating solution of the conservative system, i.e. the limiting solution (at ε tending to zero) for the required NNMs of the self-excited system (2.1).

Self-Excited Systems: Examples

TWO LINEARLY CONNECTED VAN DER POL OSCILLATORS

As a simple example, let us consider two linearly connected Van der Pol oscillators with nonlinear elastic characteristics:

q1 + q 1 + αq 3 1 + β(q 1 -q 2 ) -ε(1 -q 2 1 ) q1 = 0, q2 + q 2 + αq 3 2 + β(q 2 -q 1 ) -εγ (1 -q 2 2 ) q2 = 0. (3.1)
Here ε is a small parameter. Two vibration modes may obtained in the zero approximation by ε: q 2 = ±q 1 . Let us introduce new coordinates x, y as follows:

q 1 + q 2 = x, q 1 -q 2 = y.
In these coordinates, one of the NNMs, the so-called in-phase mode is governed by y 0 = 0, x 0 = x 0 (t) in the zeroth approximation by ε, and

ẍ0 = -x 0 - α 4 x 3 0 , ẋ2 0 = 2 h - x 2 0 2 - αx 4 0 16 , (3.2)
where h is the system energy constant. An equation describing the trajectory y(x) in the first approximation with respect to ε has the form

y 1 ẋ2 0 + y 1 ẍ0 + y 1 1 + 2β + 3α 4 x 2 0 -1 - x 2 0 4 ẋ0 2 (1 -γ ) = 0. (3.3)
The equation should be solved together with the boundary conditions of the form (2.5) (at ẋ0 = 0, x 0 = Q):

y 1 (Q) ẍ0 (Q) + y 1 (Q) 1 + 2β + 3α 4 Q 2 = 0. (3.4)
Substituting the solution y(x) as a power series into (3.3), we obtain a set of recurrent equations in the expansion coefficients. (Note that ẋ0 must be presented as a power series by x 0 too). One must also use the condition (3.4). Let α = 1.5; β = 1.5; γ = 0.5, ε = 0.2. In the first approximation, one obtains the trajectory of the in-phase vibration mode:

y 1 -0.07203 + 0.01330x 2 -0.000077x 4 .
Here x(0) = Q = 4. The last value was obtained from the potentiality condition (2.8) in the first approximation with respect to ε in a quasiharmonic approximation.

The second of the NNMs, the so-called out-of-phase mode, is governed by x 0 = 0, y 0 = y 0 (t) in the zeroth approximation by ε; respectively, The equation describing the trajectory x(y) approximation with respect to ε has the form

ÿ0 = -(1 + 2β)y 0 - α 4 y 3 0 ; ẏ2 0 = 2 h -(1 + 2β)
x 1 ẏ2 0 + x 1 ÿ0 + x 1 1 + 3α 4 y 2 0 -1 - y 2 0 4 ẏ0 2 (1 -γ ) = 0. ( 3.6)
The equation should be solved together with the boundary conditions of the form (2.7) (at ẏ0 = 0, y 0 = Q):

x 1 (Q) ÿ0 (Q) + x 1 (Q) 1 + 3α 4 Q 2 = 0. (3.7)
Substituting the solution x(y) as a power series into (3.6) and (3.7), one obtains, in the first approximation, the following trajectory of the out-phase vibration mode:

x 1 -0.10179 + 0.00881y 2 -0.000386y 4 .

Here y(0) = Q = 4. The value was obtained from the potentiality condition as above.

There is good agreement between the approximation solution and the numerical calculations checked by a computer.

PLATE VIBRATIONS IN THE FLOW OF GAS

Another example refers to the problem of plate vibrations in the flow of gas. A simple, twodegrees-of-freedom model is chosen.

Next, q 1 is the vertical displacement of the bend of the center axis of the plate; q 2 is the angle of rotation around the axis.

Both the restoring force P (q 1 ) and the restoring moment M(q 2 ) are nonlinear:

P (q 1 ) = (d 1 + d 2 q 2 1 )q 1 , M(q 2 ) = (d 3 + d 4 q 2 2 )q 2 ,
where d 1 , d 2 , d 3 , d 4 are coefficients depending on the plate's elastic characteristics. The aerodynamic lift is

A(q 2 ) = σ v 2 sin 2q 2 σ v 2 2q 2 - 4q 3 2 3 ,
where v is the flow velocity, and σ is dependent on the flow density and the airfoil aerodynamic chord.

Let ρ be the radius of inertia with respect to the bend (peck) center axis, l is the distance from the center of rigidity to the center of mass, and m is the mass of the airfoil per unit length. The equations of motion for the system are the following:

m q1 -ml q2 + (d 1 + d 2 q 2 1 )q 1 = σ v 2 2q 2 - 4q 3 2 3 -ml q1 + m(ρ 2 + l 2 ) q2 + (d 3 + d 4 q 2 2 )q 2 = 0. ( 3.8)
The analysis of the NNMs of the system (3.8) is based on the assumption that the values of d 3 and d 4 are close to those of d 1 and d 2 respectively, the sum ρ 2 + l 2 is close to unity and l 2 is much smaller than unity.

Make the coordinate transformation,

q 1 + q 2 = x √ 2, q 1 -q 2 = y √ 2.
One obtains the following equations with respect to variables x, y in place of the system (3.8):

m ẍ + d 1 x + 0.5d 2 (x 3 + 3xy 2 ) = ε ml ẍ + σ v 2 [(x -y) -(x -y) 3 /3] + 0.5m(1 -ρ 2 -l 2 )( ẍ -ÿ) + [(d 1 -d 3 ) + (d 2 -d 4 )(x -y) 2 /2](x -y)/2 , m ÿ + d 1 y + 0.5d 2 (y 3 + 3yx 2 ) = ε -ml ÿ + σ v 2 [(x -y) -(x -y) 3 /3] -0.5m(1 -ρ 2 -l 2 )( ẍ -ÿ) -[(d 1 -d 3 ) + (d 2 -d 4 )(x -y) 2 /2](x -y)/2 , (3.9)
where ε is a formal small parameter characterizing the relative smallness of the right-hand sides of Equations (3.9). Consider the in-phase NNM close to the mode q 1 = q 2 or x = x(t), y = 0 in the zeroth approximation with respect to ε. Make use of the relations (2.4) and (2.5).

The corresponding equation in the first approximation with respect to ε for the mode y(x) is

2y h 0 - d 1 2 x - d 2 8 x -y d 1 x + d 2 2 x 3 + d 1 + 3 2 d 2 x 2 y = σ v 2 x - x 3 3 + (1 -ρ 2 -1) 2 2 d 1 x + d 2 2 x 3 -(d 1 -d 3 ) + (d 2 -d 4 ) x 2 2 x 2 ,
while the boundary conditions ( ẋ = 0, x = Q) may be written as

y -d 1 x - d 2 2 x 3 + d 1 + 3 2 d 2 x 2 y -σ v 2 x - x 3 3 + (1 -ρ 2 -l 2 ) 2 d 1 x + d 2 2 x 3 -(d 1 -d 3 ) + (d 2 -d 4 ) x 2 2
x 2 x=Q , h 0 being the energy of the generating system. Numerical calculations were performed for m = 1, l = 0.2,

d 1 = 2, d 2 = 1, d 3 = 2.3, d 4 = 1.1, σ = 0.5, ρ = 1,
and v 2 = {0; 0.2; 0.4}.

The following expressions were obtained for the NNM y = y(x), x(0

) = Q = 1 at ε = 1: v 2 = 0 : y 0.0374x + 0.0018x 3 -0.00014x 5 , v 2 = 0.2 : y 0.040x + 0.009x 3 -0.0005x 5 , v 2 = 0.4 : y 0.125x + 0.022x 3 -0.001x 5 .
The potentiality condition (2.8) for the system at hand is identically satisfied for any vibration mode defined by an analytical function y(x) (or x(y)), and therefore the NNMs are dependent on two free parameters, namely the vibration amplitude Q and the phase φ.

Iterative Computation of Quasinormal Vibrations in Nonlinear Viscoelastic

Mechanical Systems

Here we will construct periodical solutions of integro-differential equations corresponding to some viscoelastic mechanical systems. One assumes that elastic forces are dominant. A conservative system corresponding to the elastic interactions is selected as generative and the NNMs of the system are selected as generative solutions.

CONSTRUCTION OF PERIODICAL SOLUTIONS IN TWO-DEGREES-OF-FREEDOM NONLINEAR MODEL

Analysis of vibrations in nonlinear viscoelastic mechanical systems is associated with studies into integro-differential partial equations [START_REF] Morgunov | Asymptotic methods for solution of some problems of nonlinear mechanics[END_REF][START_REF] Mikhlin | An application of the Ince algebraization to the stability of nonlinear normal vibration modes[END_REF]. By representing solutions of such equations as linear combinations of coordinate functions with time-dependent coefficients and by using the Bubnov-Galerkin method, sets of nonlinear integro-differential equations are obtained.

Restricting ourselves to two principal coordinate functions, we obtain the following set of equations with respect to the required coefficients of the linear combinations x 1 (t), x 2 (t):

ẍi + x i (x 1 , x 2 ) = εg i (i = 1, 2), ( 4.1)
where = (x 1 , x 2 ) is the potential energy that corresponds to the principal elastic interactions and is assumed to be an analytical function of its arguments; ε is a small parameter. Here and later on, the case of forced resonance will be called Case A, that of self-excited vibrations, Case B. In Case A

g i = g i1 (t, x 1 , x 2 , ẋ1 , ẋ2 ) + t -∞ R i (t -s)g i2 (x 1 (s), x 2 (s)) ds, while in Case B g i = g i1 (x 1 , x 2 , ẋ1 , ẋ2 ) + t -∞ R i (t -s)g i2 (x 1 (s), x 2 (s)) ds,
Here, g i1 , g i2 are analytical functions which, in Case A, are continuous and periodical in t with a period of T ; R i are relaxation kernels that satisfy the following conditions:

dR i (τ ) dτ < 0, ∞ -∞ R i (τ ) dτ = 1.
We assume that a generating system (at ε = 0), ẍi + x i (x 1 , x 2 ) = 0, allows similar NNMs of the type

x 2 = kx 1 , (4.2)
wherein the constants k are defined by the algebraic equation

k x i (x i , kx i ) = x 2 (x 1 , kx 1 ).
Without loss of generality, we assume that after some rotation of the coordinate system (an axis of a new coordinate system is directed along the rectilinear trajectory of solution (4.2)), the generating system takes on the form ẍi0 + x 1 (x 10 , x 10 ) = 0, x 20 = 0, i.e., k = 0 for the similar NNMs under condideration and the generating solution in Case A is

x 10 = x 10 (F, , t + φ), x 20 = 0, (4.3)
in Case B,

x 10 = x 10 (F, , t), x 20 = 0, (4.4)
where F and are two amplitudes of vibrations (two values of x 10 that correspond to the cusp, ẋ10 = 0); φ is the generating solution phase which, in Case B, is arbitrary due to the system being autonomous and therefore is not explicitly introduced into the solution. The amplitudes F and are related to the period of the generating solutions (4.3) and (4.4) by the equality

T (F, ) = √ 2 F dx (h -(x, 0)) 1/2 ,
where h is the energy of the generating solution, h = (F, 0) = ( , 0). Without loss of generality, one assumes that F is greater or equal to . Let us look for a solution of the original perturbed system (4.1) in the form

x 1 = x 10 + ξ, x 2 = η.
In order to find ξ and η, we obtain from (4.1), selecting a linearized part of the system with respect to ξ and η:

ξ + p 1 (x 10 )ξ = G 1 , η + p 2 (x 10 )η = G 2 , ( 4.5)
where

p 1 (x 10 ) = x 1 x i (x 10 , 0), p 2 (x 10 ) = x 2 x 2 (x 10 , 0); G 1 = εg 1 -x 1 -p 1 (x 10 )ξ, G 2 = εg 2 -x 2 -p 2 (x 10 )η, here g 1 , g 2 , x 1 , x 2 are calculated as x 1 = x 10 + ξ , x 2 = η.
It is clear that the first equation of (4.5), linearized with respect to ξ and η, is a variational equation for the generating solution. Therefore, fundamental solutions of the linearized system can be obtained by the derivation of the generating solution with respect to the arbitrary parameters, namely the amplitude F and the phase φ [START_REF] Malkin | Certain Problems of the Theory of Nonlinear Vibrations[END_REF]:

w 11 = 1 ω dx 10 dφ = γ ẋ10 (γ = const.), w 12 = 1 ω dx 10 dF ,
where

ω = 2π T (F ) , ω = dω dF .
The first fundamental solution is periodical, while the second one involves a product of some periodical function by t.

A fundamental solutions of the second linearized equation in (4.5) may be obtained for certain classes of potential functions in the form of hypergeometric functions or Lamé functions [START_REF] Manevich | The Method of Normal Oscillation for Essentially Nonlinear Systems[END_REF][START_REF] Bellman | Introduction to Matrix Analysts[END_REF][START_REF] Dmitrenko | Iterative calculation of self-induced oscillations of essentially non-linear viscoelastic systems[END_REF]. In other cases, it may be constructed by a power series expansion in x 10 or by other methods, see [START_REF] Mikhlin | Normal vibrations of a general class of conservative oscillators[END_REF][START_REF] Bellman | Introduction to Matrix Analysts[END_REF]. Let us denote the fundamental solutions by w 21 , w 22 .

ITERATION PROCESS

In order to find periodical quasinormal solutions of Equations (4.5), we employ here the following iteration process: Besides: the periodicity conditions for ξ k+1 and ξk+1 must be written in Case A:

ξ k+1 = 1 1 t 0 G τ 1k (
Here i = w i1 w i2 ẇi1 ẇi2 = const.; f θ ij (z) = f ij (z)| z=θ ; ξ (n) k (n = 1, . . . ,
T 0 G τ 1k w τ 11 dτ = 0, T = T (F, ) = √ w F dx ( (F, 0) -(x, 0)) 1/2 (4.8)
(the period T is known); (F, 0) = ( , 0). From this set, the phase φ may be found in the corresponding approximation of the perturbation method, together with the amplitudes F , of the generating solution.

In Case B, the periodicity conditions for ξ k+1 and ξk+1 have a slightly different form:

T k 0 G τ 1k w τ 11 dτ = 0, T = T k = √ w F k k dx ( (F k , 0) -(x, 0)) 1/2
(4.9)

(the period T k is not known); (F k , 0) = ( k , 0). Equations (4.9) may evaluate the period T k and the amplitudes F k , k of generating solution in the corresponding approximation with respect to ε.

It is assumed that all variables and parameters vary over a limited region, all functions in (4.1) and their derivatives, inclusive of those of third-order, are bounded, satisfy the Lipschitz condition, and w ij , ẇij , ij (i, j = 1, 2) are also bounded; the remaining conditions are given in [START_REF] Mikhlin | Normal vibrations of finite-dimensional non-linear systems[END_REF]30].

For the iterations of (4.6) to converge to the solution of the initial integro-differential equations,the following principal conditions are necessary: the variables ξ k+1 , η k+1 , ξk+1 , ηk+1 should remain in a bounded region as should the phase φ k and the amplitudes F and in Case A or the period T k and the amplitudes F k , k in Case B. The operators that determine the transformation (4.6) should be contracting; the periodicity conditions (4.8) or (4.9) should be uniquely solvable with respect to φ k , F , or T k , F k , k ; and the corresponding operators determining the phase, amplitude, and period of the vibrations should be contracting.

These requirements lead to a set of inequalities that make the small parameter ε bounded above [START_REF] Mikhlin | Normal vibrations of finite-dimensional non-linear systems[END_REF]30]; the inequalities are cumbersome and therefore are not shown here.

In the first approximation of the iteration process, the number of periodic solutions of the set of integro-differential equations is determined by the number of roots of Equations (4.8) in Case A or the number of roots F k of Equations (4.9) in Case B. When ε tends towards zero, all these solutions tend towards NNMs of the generating conservative (elastic) system.

RESONANCE NONLINEAR TRANSVERSAL VIBRATIONS OF A VISCOELASTIC ROD

By way of example, let us consider resonance nonlinear transversal vibrations of a viscoelastic rod. Within the framework of the well-known Kirchhof hypothesis, the equations of motion are taken in the form

u tt + EI u xxxx - EI 2F I 0 (u ξ ) 2 ]rmdξ u xx = ελ sin pt ∞ k=1 a k sin kπ x 1 + ε t -∞ R(t -s)u 3 (s) ds,
where u is the transversal displacement; x is the longitudinal independent variable; E is the Young modulus; I is the moment of inertia, l is a distance between end points, R(z) is the relaxation kernel, F is an area of the transversal cross-section of the rod, and ε is a small parameter. The rod ends are assumed to be restrained.

Setting the spatial wave formation in the form (m and n are integers)

u = x 1 (t) sin mπ x 1 + x 2 (t) sin nπ x 1 ,
by using the standard Bubnov-Galerkin method one obtains:

ẍ1 + EF 4 mπ 1 4 x 3 1 + mπ 1 2 nπ 1 2 x 1 x 2 2 = -EI mπ 1 4 x 1 + [sin pt + J (x 1 , x 2 )], ẍ2 + EF 4 nπ 1 4 x 3 2 + mπ 1 2 nπ 1 2 x 2 x 2 1 = -EI nπ 1 4 x 2 + ε[sin pt + J (x 2 , x 1 )].
Here

J (u, v) = t -∞ R(t -s) 3 4 u 3 (s) + 3 2 uv 2 ds, λa n = λa m = 1.
In the zeroth approximation, ε = 0, there are two nonlinear normal vibration modes with rectilinear trajectories (similar NNMs), x 1 = 0 and x 2 = 0. A numerical computation using the presented analytical approach was performed by Petrov for a rod of acrylic plastic ST-1 having E = 6.3 × 10 9 N m -2 and l = 0.4 m; the rod had a square cross-section with a side of 0.02 m; the representation R(τ ) = A e -βτ /τ 1-α , with A = 0.0286, β = 0.05 and α = 0.075 was assumed for the relaxation kernel; n = 1, m = 2. Parameters of the resonance modes close to the mode x 2 = 0 were calculated in the zeroth and the first approximations with respect to ε. Some variation of the analytical procedure was introduced here, namely, the amplitude F 0 of the zeroth approximation was given and the period T was calculated.

When the amplitude F 0 of the zeroth approximation is equal to 5, the period T = 1.916 × 10 -3 s and the phase shift between the zeroth and the first approximation solutions, x 10 and x 11 is: π/2 < φ < π. When F 0 = 8, the calculations yield T = 1.197 × 10 -3 s and the phase shift between x 10 and x 11 lies within the above range.

Conclusions and Discussion

At the present time, it is known that NNMs are typical periodical solutions in n-degreesof-freedom nonlinear conservative systems. Moreover, normal or quasinormal vibrations exist in broad classes of nonlinear, near-conservative systems, as was indicated above, in nonautonomous systems (forced resonances), self-excited systems and systems containing viscoelastic interactions which are presented in the corresponding equations in a form of integral operators. It is found that approaches for the analysis of NNMs in conservative systems are applicable to near-conservative nonlinear finite-dimensional systems.

It is fruitful to extend the concept of NNMs to distributed nonlinear systems. One of possible extensions of the concept to continuous oscillators was performed by Vakakis [START_REF] Manevich | The Method of Normal Oscillation for Essentially Nonlinear Systems[END_REF]. In the spirit of Rosenberg [START_REF] Rosenberg | The normal modes of nonlinear n-degrees-of-freedom systems[END_REF][START_REF] Rosenberg | Nonlinear vibrations of systems with many degrees of freedom[END_REF] and Shaw and Pierre [START_REF] Shaw | Nonlinear normal modes and invariant manifolds[END_REF], the continuous nonlinear normal modes were defined as motions during which all material points of the system vary equiperiodically, vanishing or reaching their extremal values at the same instant in time. This way leads to equations which are similar to the equations considered previously in the case of conservative or near-conservative systems.

Other generalizations of the NNMs to distributed systems are possible too:

1. It is possible to determine NNMs in distributed systems as solutions in separated variables with a uniform time function. The solutions correspond to similar normal modes in finitedimensional systems with rectilinear trajectories in the configuration space. 2. During NNM, a finite-dimensional system behaves like a single-DOF conservative one.

Therefore it is advisable to introduce nonlinear stationary traveling waves with a single phase (these are so-called simple waves) as a new generalization of NNMs to the distribution case. 3. There is a coupling of the existence of NNMs and symmetry properties of finitedimensional systems [START_REF] Manevich | The Method of Normal Oscillation for Essentially Nonlinear Systems[END_REF]. It is possible to determine NNMs in a distributed system as solutions which are invariant with respect to the symmetry groups of the system.

Note, finally, that problems of NNMs extending to distributed nonlinear systems as a development of approaches for the analysis of the solutions, are interesting, but no details are available at present.

4 )

 4 are substituted into the right-hand sides of the above relationships (4.6) together with φ k (Case A) or F k (Case B); D ik are arbitrary constants. Let us use the conditions of T -periodicity of the functions η k+1 , ηk+1 , to find the unknown constants D 1k , D 2k (the period T is given in Case A and unknown in Case B): D 1k (w T 21w 0 21 ) + D 2k (w T 22 -