Stefan-Gabriel Chitic

Julien Ponge
email: julien.ponge@insa-lyon.fr

Olivier Simonin
email: olivier.simonin@insa-lyon.fr

SDfR -Service discovery for multi-robot systems

Keywords: Robotic fleet, Service Discovery, Service oriented Architecture

INTRODUCTION

Nowadays more and more robotic systems tend to be composed of several robots moving and cooperating, generally called fleets of robots. They are able to perform one or multiple tasks together and to share information about complex missions.

One of the problems in developing multi-robot applications is the communication inside the fleet. The robots need to know the reachable peers in any type of communication infrastructure. In an ad-hoc network, this problem becomes challenging because the robots are highly mobile. A fleet can easily partition or merge with another fleet. Also single robots can become isolated or join an existing fleet.

There is a great need of having a suitable middleware that can abstract the neighbor and service discovery layer and offer high-level application interfaces (APIs) to robotic applications. The challenge is to provide an adapted system that can offer a suitable mechanism for network configuration. Second, this mechanism needs to maintain a list of connected peers and their services in a highly mobile scenario.

In this paper we propose a new service discovery middleware, called SDfR (Service Discovery for Robots), an extension of the Simple Service Discovery Protocol (SSDP) used in Universal Plug and Play (UPnP). SDfR is able to make an auto-configuration for the connectivity to the fleet ad-hoc network, dis-cover reachable neighbors and their services. We also provide a ROS integration node, simplifying the usage of SDfR.

The paper is structured as follows: Section 2 presents the background of service orientated architecture and discovery in multi-robot systems, Section 3 discusses the service discovery approaches in middleware for robotics, Section 4 presents our proposal for service discovery in a fleet context, Section 5 presents the arhitecture and the implementation of SDfR, Section 6 presents the main results of overhead benchmarking and Section 7 describes a couple of use-cases. Section 8 concludes the paper.

MOTIVATION

One problem of robotic fleets is how functionalities can be applied to make the collaboration between peers easier. Robots in fleet system that perform a mission together need to communicate with other peers. The communication can be done using a centralized node or directly using an ad-hoc network.

Robotic applications as services In our vision, robots need to advertise their functionalities as services in order to allow other members of the fleet to interact with them. In network based applications, service-oriented programming is now a largely accepted principle [START_REF] Issarny | Service-oriented middleware for the future internet: state of the art and research directions[END_REF].

Service-oriented architecture greatly simplifies the implementation of highly-adaptive, constantlyevolving applications [START_REF] Frénot | Various Extensions for the Ambient OSGi framework[END_REF]. It also reduces the process of developing and deploying new robotic applications. This architecture is very suitable to quickly cope with new developing models and requirements.

In order to increase the mobility of the robots and to distribute the communication without having a central node, there is a need for the communication to be decentralized using ad-hoc networks. In this case, the robots do not have any image of their neighbors. Furthermore, the communication across peers is susceptible to route change and different peers can be used to relay a data package. The ad-hoc network becomes the sum of peer to peer network with 2 robots at least.

Fleet service discovery

In order to solve the problem of neighbors and service discovery in an ad-hoc network, a fleet robot needs a protocol that is able to constantly discover new robots in its coverage area, while maintaining a neighbor connectivity quality indicator. Since there is not any central node that can manage IP address allocation, the protocol should be able to negotiate an IP address inside the network and to have a conflict management tool in case of an IP collision. Our work focuses on ad-hoc multirobot systems communication, especially on a family of middleware that accelerates the development and the deployment of new robotic applications.

RELATED WORK

In this section we present some of the major middlewares in robotics as well as some of existing service discovery protocols used in the field.

Middleware for robotics

All aspects of communication, application deployment and configuration can be facilitated using a proper middleware. The biggest difference between a classical middleware that runs in a cloud centralized infrastructure and a robotic one, is the mobility of the fleet and the decentralization of its components. An exhaustive survey of the existing middlewares for robot contexts is clearly impossible because of the large number of existing middleware and frequent releases of new ones. A detailed survey can be found in [START_REF] Chitic | Are middlewares ready for multi-robots systems?[END_REF]. Most of the middlewares are designed for mono-robot with little applicability in multi-robot systems.

One of the first middleware that emerged for robotics is Player/Stage [START_REF] Kranz | A player/stage system for context-aware intelligent environments[END_REF] project which is designed to provide an infrastructure, drivers and a collection of dynamically loaded device-shared libraries for robotic applications. It neither offers a service orientated architecture, nor a discovery for neighbors services system.

Another highly used middleware for robotics is ROS (Robot operating system). It is a recent flexible middleware for robot applications [START_REF] Quigley | Ros: an open-source robot operating system[END_REF][START_REF] Cousins | Sharing software with ros [ros topics[END_REF][START_REF] Ros | Robot operating system[END_REF]. It is a collection of tools, libraries, and conventions that aim to simplify the task of creating complex and robust robot behavior across a wide variety of robotic platforms. It provides hardware abstraction, device drivers, visualizers, message-passing, package management.

ROS (ROS, 2014) seems to be the emerging middleware for mono-robot with the most potential to become the most used framework for robotic fleets. However, it has no multi-robot coordination system and no automated testing environment, but it has already the advantage of having a large community that develops new packages for it.

To the best of our knowledge, none of the existing middleware for robots provides a service discovery mechanism. ROS internal repository is represented as an internal node designed for internal service discovery, but it does not support neighbors discovery.

Service discovery and robotics

Classical protocols and middlewares for service discovery in distributed environments like data-grids, clouds or even smart environments have a centralized or decentralized registry that manages service description. Decentralized systems (e.g UPnP [START_REF] Ahn | Upnp approach for robot middleware[END_REF], Jini [START_REF] Pereira | Peer-to-peer Jini for truly service-oriented WSNs[END_REF] or SLP [START_REF] Romero | Service discovery in ubiquitous feedback control loops[END_REF]) can be a purely distributed solution where each node stores its own service repositories or a hybrid solution that includes super-nodes that aggregate information from other peers.

Another way to see a fleet of robots is like a service-oriented multi-agent system. Such environments like Peer-to-Peer (P2P), Multi-Agent Systems (MAS) or Service-Oriented Environments (SOE) tend to approach the problem of service discovery in a centralized, distributed or decentralized way. Centralized mechanisms like super-peers [START_REF] Gummadi | A measurement study of napster and gnutella as examples of peer-to-peer file sharing systems[END_REF], middle-agents [START_REF] Klusch | Automated semantic web service discovery with owls-mx[END_REF] or central registries [START_REF] Rompothong | A query federation of uddi registries[END_REF] are limited in number of agents in the system and in terms of the number of requests. They also use a centralized node which can have serious impact if the central point becomes unreachable. Distributed approaches such as Distributed Hash Tables (DHT) [START_REF] Maymounkov | Kademlia: A peer-to-peer information system based on the xor metric[END_REF] offer more scalability and robustness by having multiple specific nodes that can manage the resources. Decentralized systems consider all the nodes to be equal. This approach provides more flexibility, but it has its downsides, since each node has only a partial view of the entire system. As mentioned in [START_REF] Del Val | Enhancing decentralized service discovery in open serviceoriented multi-agent systems[END_REF], an interesting way to discover service inside a decentralized and self-organized multi-agent system is to use homogeneity between agents.

In the robotic world, an approach for service discovery in centralized networks is to use classical Universal Plug and Play (UPnP) protocol. Since the concept of having the robotic tasks and processes as services is not mature yet, the main focus of research on service discovery in robotics is oriented towards the integration with the environment, like smart homes or smart cities. Reference [START_REF] Borja | Integration of service robots in the smart home by means of UPnP: A surveillance robot case study[END_REF] provides a case study of integration of service robots and smarthomes via UPnP. In these cases, the authors are not referring to a robot as part of a specific fleet, but as part of an environment, in which the robot is considered as an entity that can offer services.

The solutions presented above have their downsides when applied to ad-hoc multi-robot systems. Firstly, due to the mobility of the robots, the network connection is highly instable and robots can disconnect and reconnect very often. Existing protocols do not perform the same way in a highly dynamic environment and in a static one. As mentioned in [START_REF] Issarny | Service-oriented middleware for the future internet: state of the art and research directions[END_REF], the challenge is to set the trade-off between physical mobility and scalability. Secondly, existing protocols are not very adaptive. The discovery protocol should be ready to be used at any time and track its usage and failures. Existing protocols like UPnP, have a limited memory factor and when the connection is timed-out, the discovery process reinitializes itself at reconnection.

Next sections present our extension of SDP used in UPnP in order to overcome the limitation of existing service discovery approaches.

SERVICE DISCOVERY FOR ROBOTICS

Our main goal is to propose a mechanism that allows highly mobile robots to keep track of the reachable peers inside a fleet while using an ad-hoc infrastructure. This mechanism should also be able to provide a list of services available on each peer. Another objective is to propose a network configuration negotiation protocol. Due to the mobility of robots, classical peer to peer network configuration techniques are not suitable.

In this section we present the general description of our service discovery protocol for robotic applications, called SDfR (Service Discovery for Robots).

SDfR vs SSDP

We propose a protocol that does not flood the network and provides an already seen memory feature built-in. SDfR protocol aims to be a dynamic, adaptive and scalable protocol. For this purpose we based our approach on Simple Service Discovery Protocol (Goland et al., ) (SSDP) from UPnP. SDfR also aims to provide service discovery with the smartenvironment in which the robots are being deployed.

In order to limit the network usage for the service discovery process, SDfR sends most of the internal messages in multicast, avoiding the overhead generated by unicast transmission in order to propagate the same message. In addition, to avoid failure in case of a disconnection due to the movement of the robots outside the coverage area, all the communications are done using UDP (User Datagram Protocol). Furthermore, to limit the network flooding when the protocol needs information from just one robot, a second transmission is enabled in unicast mode. SDfR does not need to reinitialize the entire discovery protocol when the connection is lost, because it disposes of a history map of all the already seen robots and their services.

In order to avoid services that are out of reach (e.g. service of robots that are present in the history map but are not present in the covered communication area), a connection indicator is computed for each robot represented by the success rate of pinging the connected peers. The main similarities between SSDP and SDfR are:

Multicast transmissions In order to avoid the overhead of retransmitting the same unicast message, most of the internal messages are multicast. HTTP-style messages The messages that are being sent use an HTTP-style structure composed of headers and a body.

The main differences between SSDP and SDfR are:

Limited multicast transmissions To avoid failure in case of a disconnection due to the movement of the robots outside the coverage area, all the communications in SDfR are done using UDP (User Datagram Protocol) and only in request model. Unicast transmissions To limit the network flooding when the protocol needs information from just one robot, a second transmission mechanism is enabled in unicast mode in SDfR.

History map SDfR does not need to reinitialize the entire discovery protocol when the connection is lost, because it disposes of a history map of all the already seen robots and their services. In order to avoid services that are out of reach (e.g service of robots that are present in the history map but are not present in the covered communication area), a connection indicator is computed for each robot.

Protocol Description

We designed our protocol as a finite state timed automaton. A detailed description of the design process and the model validation is out of the scope of this paper.

Each protocol method has at least one type of message that resides inside the request payload. A sample execution of the communications is illustrated in Figure 1. Two methods representing the desired action of a request are used in SDfR : M-SEARCH and NO-TIFY.

The M-SEARCH method is used for discovery requests to get the list of near by members and their services. The message type associated with this method is Discovery.

The other method, NOTIFY is used to respond to a Discovery request or to inform the others about changes in the current state of the robot. The message types associated with this method are: Update, Alive and Byebye.

• The Update message is sent as a response to a Discovery request or when the current services or capacities of the robot change.

• The Alive message is sent recurrently, as a beacon, in order to inform the others about the presence of the robot. The beacon sending rate can be set depending on the services need. The default value is at 10s.

• The Byebye message is sent when the robot stops gracefully, in order to inform the others about its disappearance.

ARCHITECTURE OF SDFR

In this section we zoom in on the architecture and features embedded in SDfR by presenting its implementation as a service.

SDfR Architecture

Service Discovery for Robots is developed as a service itself. The service-oriented architecture approach for robotic software development is not very widely used in the robotic community. The practice in this community is to develop built-in libraries in order to extend software features. The main advantage of having a service-oriented architecture is the compatibility with other robotic services developed in different programming languages and running over different operating systems. This is a critical feature for a heterogeneous robotic fleet.

Furthermore, SDfR service can run separately of the other processes on the robot and all the messages are consumed by instances of the service on multiple robots. If it fails, it would not affect the other services running on the robot. This sand-boxing also ensures that the information sent by the protocol is not corrupted by any other third-parties.

SDfR service is composed of two layers as shown in Figure 2: an API layer that communicates with other services and a Discovery Protocol layer. Each layer has an independent life-cycle and communicates internally via a shared memory. The API layer responds to requests independently from the lower layer, using the information from the shared memory. The lower discovery protocol layer is in charge of communication with the other SDfR nodes on an elastic number of robots in order to discover the reachable peers and their services.

RESTful communication API A Representational

State Transfer [START_REF] Fielding | Representational state transfer. Architectural Styles and the Design of Netowork-based Software Architecture[END_REF] (RESTful) webservice was chosen for the API. It lets other services communicate with SDfR service because it is based on normal HTTP requests which enables intent to be inferred from the type of request being made and Ad-hoc configuration management Since the fleet is operating in an ad-hoc infrastructure, the peers need to be able to negotiate and auto-configure their network configuration. A robotic fleet ad-hoc network is different from a classic ad-hoc hot-spot because the robots are highly dynamic and the network can be partitioned or merged. The mobility of the peers needs to be taken into consideration in the negotiation protocol of the configuration. SDfR service, based on a simple configuration file, is able to automatically connect to an ad-hoc network. The secured WiFi SSID network is composed using the fleet id. This mechanism allows to have multiple fleets of robots in the same network space. Moreover, the robots can auto assign IP addresses. The standard network space is 10.<fleet id>.<x>.<y>, where x and y are computed each robot from its internal MAC address in order to avoid IP conflict [START_REF] Thomson | Ipv6 stateless address autoconfiguration[END_REF]. Furthermore, if an IP conflict happens, the service has a mechanism to trigger an IP change on the robots. This mechanism is available all the time since an IP conflict can be triggered by a merge of 2 subnetworks.

ROS integration

In order to make SDfR service user friendly, we have created a ROS node that communicates with SDfR service and can be used by other nodes via topics and services. When the node starts, it launches a instance of SDfR if it is not running and then it provides support for other ROS nodes to publish or unpublish their services and their capacities. Furthermore, our ROS node provides the neighbors list of services and is capable of allowing other ROS applications to search for a specific service with a specific configuration.

A producer node can publish its services or capacities in an asynchronous way using ROS topics be-cause we consider that the registration is not highly bound to time. The same concept applies for the unpublish commands and for getting the list of neighbors. On the other hand, the search command for a specific neighbor and their services needs to be done in synchronous way using ROS services because the behaviors of the consumer node is depending on it.

EVALUATION

In the evaluation we want to measure what is the overload generated by the usage of SDfR in a robotic fleet context. Another objective is to see the impact of using text-plain protocol in the upper and the lower layer of the SDfR service. In this section we present the main results of resource overhead for SDfR. SDfR was benchmarked in two types of contexts:

• a static scenario where the robots do not move to evaluate the overhead in an ideal Wi-Fi communication scheme

• a real dynamic scenario (see Figure 3) where robots are moving and transmission can drop.

In both scenarios all peers should discover their neighbors and in the the second one, the neighbors discovery depends on the distance between peers.

The benchmarks were performed on Turtlebot 2 robots equipped with an Intel Core 2 Duo, 2.1 GHz CPU, 4Gb of Ram PC running on Ubuntu 13.04

Each test run was given 5 minutes to collect the data. In our test runs we used simulated services that try to register/subscribe into SDfR. We simulate three type of actions:

1. Publish. We simulate new service providers that try to publish with a delay time of 10 seconds. In order to simulate publishers, we used an Apache server on each robot that responds to the autodiscovery URL of each publisher.

2. Unpublish. Each of the already published service providers could be unpublished with a random delay between 5 seconds.

3. Subscribe. We generate separated threads for each consumer. Each thread will constantly request the table of neighbors from SDfR, in order to stress at maximum our protocol.

In the static scenario we considered different numbers of robots: 2, 4, 6. Each robot had a total number of service providers and service consumers equal to 100 simulated services. For each number of robots we used different ratios between providers and consumers: 30%, 50% and 70% publishers.

In the dynamic test-case we used for each variant of SDfR 6 robots with a number of 100 simulated services per robot. 70% of the services on each robot were publishers. The robots moved randomly in a 200 square meters room with poles and other obstacles. The room (see Figure 3) was exposed to Wi-Fi interferences from other networks that occupy all of the 2.4Ghz channels.

Overhead evaluation

A robotic application that provides a service for the fleet needs to register with SDfR. This must be done as fast as possible in order to avoid blocking the service when it starts. Another important overhead measure is the time of unregister request. This happens whenever a producer wants to remove itself from the SDfR registry. One of the most important metric for robotics application from a latency point of view is the time to request the list of neighbors and their services. In a real scenario, a producer registers once for its life-time cycle, but a consumer may request multiple times the list of reachable neighbors and services. The subscription response time can have a impact on the time to complete a fleet task.

Figure 4 reveals the resulting time consumed by a service to publish or unpublish its service or get the neighbors services in the static and dynamic scenar-ios. In the dynamic scenario (Figure 4), the latency is higher but the difference from the static scenario is less than 30%. This can be explained by the increase in computation load on robots generated by the mobility of the fleet. All the request times remain in a certain variation interval even if the number of robots increases. The publish request takes from 95ms to 130ms, the unpublish request varies between 98ms to 122ms and the subscribe request is less than 5ms.

CPU and Memory usage evaluation

In robotic applications, the computation power is critical. Fleet of robots are heterogeneous and can include different types of robots with different computational factor. Having a low CPU consumption discovery service, benefits the other processes involved in performing the fleet mission. Besides the CPU usage, another critical resource in fleets of robots is the memory. As an example, robotic fleets may include visual sensors like 3D cameras which are in high demand of memory. A service discovery protocol needs to have a low usage of the robot memory.

Figure 5 presents our results for % of CPU and Memory used. CPU consumption varies between 2.2% to 2.8% of CPU usage. We observe that the CPU usage is higher for the dynamic scenario because the CPU time is consumed by the mobility management and the number of CPU slots is less for other processes. We can also say that the CPU usage increases with the number of robots in the fleet with a rate of maximum 1% per robot. The memory usage is less than 0.3% (out of 4GB) and seems to remain constant even if the number of robots increases.

Latency evaluation

In a fleet context the communication between peers in ad-hoc network are very sensitive. The transmissions can be unreliable due to the mobility of the robots. This is why the network overhead needs to be as limited as possible in order to allow services to exchange information. Our measurements include the average of transmitted and received bytes per robot in a test run of 5 minutes. To increase the quality of our measurements we use intermediary check-points for each metric at each 10 seconds. Figure 6 presents the quantity of transmitted (TX) and received (RX) kilobytes per robot. The number of bytes varies from 220 kilobytes to 320 kilobytes for the transmitted bytes and from 540 kilobytes to 850 kilobytes for received bytes. This remains very limited considering the time of 5 minutes.

Conclusion

Our results show that SDfR is low latency brick, with a small overhead on CPU and memory that has a low network consumption. Moreover, the increase in number of robots in a fleet do not have a big impact on SDfR overhead. Being integrated with ROS makes SDfR a good asset in managing neighbors and service in multi-robot context.

USE-CASES

This section presents use-cases where SDfR and service oriented architecture were applied. These two use-cases are being studied in an ongoing project: 'CROME -Multi-view multi-robot scene understanding and fleet coordination 1 '. But the applicability of SoA and SDfR is larger than the examples presented below.

In general, we can use services in any robotic application that need to perform one or multiple tasks at 1 http://youtu.be/5kmxIgLEerQ once by having multiple nodes that manage different pars of hardware layers. In multi-robot context, SDfR can be used to simply discover the peers in our communication zone and avoid overloading the network with broadcasts or messages for unreachable robots. Furthermore SDfR provides a list of available services on each peer, thus a robot can decide if it can relay on a peer to perform together a task.

Peer-to-Peer monitoring in a fleet

This scenario consists in having a fleet of robots performing a mission in a decentralized network. The environment is big enough to make direct communication between all peers impossible.

In these condition, a monitor station cannot have a global image of what is happening in the fleet. In order to provide the monitor with information about each robot status, status messages and alerts need to be routed in a peer to peer context to the monitoring station.

In this case, on each robot inside the fleet, services are represented by the nodes that manage the hardware elements of each robot and that expose their status. In addition, a peer to peer monitoring service is present. It collects the data from other nodes and analyzes it in order to trigger alarms if necessary.

SDfR is used for all the hardware nodes to register in order to be discovered by the monitoring service. The monitoring service first registers itself into SDfR and then it gets constantly the lists of available service on the robot. In addition, this service uses SDfR to discover the peers that are in the communication area.

The peer-to-peer monitoring service searches for the neighbors with an instance of this same service. If it finds one, it routes in a state-less connection all the data from its nodes to the neighbor. It also tries to build a communication route for alert to get in a state-full environment to the monitoring station.

Complex scene observation by robotic fleet

This scenario consists in having a fleet of robots observing a complex scene using multiple sensors in order to have a better observation [START_REF] Matignon | Concentric and incremental multi-robot mapping to observe complex scenes[END_REF]. Combing the perception of each robot, the fleet has a better vision of the scene. The robots need to move in order to maximize the observation of the scene. In a particular case, the fleet observe the pose of a human and his joints.

In this case, on each robot we have a node that gets the sensor perception, a node that manages the movement of the robot and a upper AI layer service that decides if the robot needs to move in order to maximize the vision.

SDfR is used by the AI layer to discover the movement node and the perception node. Furthermore, SDfR is used to search for the other reachable peers in the fleet. Using SDfR allows an elastic number of robots to participate in the fleet in this type of mission.

CONCLUSION

In this paper we presented the challenges to make a service discovery protocol for robot fleet systems. We discussed the limited applicability of existing service discovery protocols in the context of robot fleets and then, we proposed a new protocol called SDfR that is suitable for service discovery inside an ad-hoc networked fleet.

The evaluation results shows that SDfR API generates a small latency and the SDfR service has a small impact on the CPU and memory used. Furthermore, SDfR has a low bandwidth consumption in both static and dynamic scenarios. Base on our benchmarking, we believe that SDfR can provide neighbor and service discovery for multi-robot applications without stressing the system. ROS integration and adhoc network auto-connect features facilitate the usage of SDfR.

We intend to further use SDfR in our current approach of improving the families of middleware to facilitate the development of multi-robot systems. SDfR can be used as part of bigger systems like remote application deployment environment in ad-hoc robotic network where robots can be dynamically provisioned with new applications while remaining in operation environment. SDfR can also be used to monitor the running services on each robot. We are also interested in opportunities for SDfR not just for robotics but also in other fields like mobile wireless sensor networks.

Figure 1 :

 1 Figure 1: Sample execution of SDfR communications

Figure 2 :

 2 Figure 2: SDfR service architecture.

Figure 3 :

 3 Figure 3: Turtlebots in the experimentation hall of CITI-Lab.

Figure 4 :

 4 Figure 4: Average request time in SDfR.

Figure 5 :

 5 Figure 5: CPU and Memory usage.

Figure 6 :

 6 Figure 6: Transmitted and Received bytes per robot in 5 min.


Stefan Bilbao

Olivier Thomas

Cyril Touzé

Michele Ducceschi

Conservative Numerical Methods for the Full von Kármán Plate Equations

This article is concerned with the numerical solution of the full dynamical von Kármán plate equations for geometrically nonlinear (large-amplitude) vibration in the simple case of a rectangular plate under periodic boundary conditions. This system is composed of three equations describing the time evolution of the transverse displacement field, as well as the two longitudinal displacements. Particular emphasis is put on developing a family of numerical schemes which, when losses are absent, are exactly energy conserving. The methodology thus extends previous work on the simple von Kármán system, for which longitudinal inertia effects are neglected, resulting in a set of two equations for the transverse displacement and an Airy stress function. Both the semi-discrete (in time) and fully discrete schemes are developed. From the numerical energy conservation property, it is possible to arrive at sufficient conditions for numerical stability, under strongly nonlinear conditions. Simulation results are presented, illustrating various

features of plate vibration at high amplitudes, as well as the numerical energy conservation property, using both simple finite difference as well as Fourier spectral discretisations. c ??? John Wiley & Sons, Inc.

I. INTRODUCTION

Thin plates and shells vibrating at large amplitude display geometric nonlinearity when the vibration amplitude is at least of the order of the thickness. Various interesting nonlinear phenomena are exhibited which can be classified according to the amplitude of the nonlinear vibration.

• For weakly nonlinear vibrations in the low-frequency part of the eigenspectrum, jump phenomena and hysteresis due to Duffing-like responses are well documented [START_REF] Yamaki | Influence of large amplitudes on flexural vibrations of elastic plates[END_REF][START_REF] Kung | Nonlinear flexural vibration of a clamped circular plate[END_REF][START_REF] Nayfeh | Nonlinear oscillations[END_REF]. Mode coupling and energy exchange between internally resonant modes is a second common feature, now well established in the literature [START_REF] Sridhar | Non-linear resonances in the forced responses of plates, part I: symmetric responses of circular plates[END_REF][START_REF] Hadian | Modal interaction in circular plates[END_REF][START_REF] Chang | Non-linear vibrations and chaos in harmonically excited rectangular plates with one-to-one internal resonance[END_REF][START_REF] Touzé | Asymmetric non-linear forced vibrations of free-edge circular plates, part I: theory[END_REF][START_REF] Thomas | Asymmetric non-linear forced vibrations of free-edge circular plates, part II: experiments[END_REF][START_REF] Amabili | Theory and experiments for large-amplitude vibrations of rectangular plates with geometric imperfections[END_REF][START_REF] Amabili | Nonlinear vibrations and stability of shells and plates[END_REF][START_REF] Ducceschi | Nonlinear dynamics of rectangular plates: investigation of modal interaction in free and forced vibrations[END_REF].

• Under strongly nonlinear conditions, a transition to a chaotic or wave turbulent regime is observed. This transition scenario is documented in [START_REF] Touzé | Transition to chaotic vibrations for harmonically forced perfect and imperfect circular plates[END_REF][START_REF] Touzé | Transition scenario to turbulence in thin vibrating plates[END_REF][START_REF] Awrejcewicz | Feigenbaum scenario exhibited by thin plate dynamics[END_REF], while spatio-temporal dynamics [START_REF] Awrejcewicz | Spatio-temporal chaos and solitons exhibited by von Kármán model[END_REF] and the existence of inertial manifolds and attractors are mathematically demonstrated in [START_REF] Chueshov | Inertial manifolds for von Kármán plate equations[END_REF][START_REF] Chueshov | Attractors and long time behavior of von Kármán thermoelastic plates[END_REF].

• In the strongly nonlinear or wave turbulent regime, an energy cascade from the injection to the disspative scale has been recently demonstrated theoretically [START_REF] Düring | Weak turbulence for a vibrating plate: Can one hear a Kolmogorov spectrum?[END_REF], and then studied experimentally [START_REF] Boudaoud | Observation of wave turbulence in vibrating plates[END_REF][START_REF] Mordant | Are there waves in elastic wave turbulence ?[END_REF][START_REF] Miquel | Nonlinear dynamics of flexural wave turbulence[END_REF] and numerically [START_REF] Düring | Weak turbulence for a vibrating plate: Can one hear a Kolmogorov spectrum?[END_REF][START_REF] Humbert | Wave turbulence in vibrating plates: The effect of damping[END_REF][START_REF] Yokoyama | Weak and strong wave turbulence spectra for elastic thin plate[END_REF][START_REF] Miquel | Transition from wave turbulence to dynamical crumpling in vibrated elastic plates[END_REF][START_REF] Ducceschi | Dynamics of the wave turbulence spectrum in vibrating plates: A numerical investigation using a conservative finite difference scheme[END_REF].

Most work on the nonlinear dynamics of thin plates uses the von Kármán assumptions, which introduce a particular truncation in the strain-displacement relationship [START_REF] Von Kármán | Festigkeitsprobleme im maschinenbau[END_REF][START_REF] Chu | Influence of large amplitudes on free flexural vibrations of rectangular elastic plates[END_REF][START_REF] Chia | Nonlinear Analysis of Plates[END_REF][START_REF] Thomas | Geometrically nonlinear flexural vibrations of plates: In-plane boundary conditions and some symmetry properties[END_REF]. More specifically, a simplified version of the von Kármán plate equations is used, where in-plane inertia is neglected, such that an Airy stress function can be used to describe in-plane motion. This approximation was first introduced by Föppl [START_REF] Föppl | Vorlesungen über technische Mechanik[END_REF][START_REF] Landau | Theory of Elasticity[END_REF][START_REF] Bazant | Stability of structures[END_REF] and the simplified model is often referred to as the Föppl-von Kármán system. In this contribution, this model will be named the simple von Kármán system. It is known to give good results in the moderately nonlinear regime, and in the low-frequency part of the eigenspectrum of the plate. However, in recent studies, numerical simulations are carried out up to a vibration amplitude that is 5 to sometimes 10 times the thickness, or, in the context of wave turbulence, down to small wavelengths where in-plane inertia may play a role [START_REF] Ducceschi | Nonlinear dynamics of rectangular plates: investigation of modal interaction in free and forced vibrations[END_REF][START_REF] Yokoyama | Weak and strong wave turbulence spectra for elastic thin plate[END_REF][START_REF] Miquel | Transition from wave turbulence to dynamical crumpling in vibrated elastic plates[END_REF]. Very few studies are available for checking the validity limits of the simple von Kármán plate model versus a full von Kármán system where in-plane inertia is retained, or versus more complex models including for example rotatory inertia or shear effects, following e.g. Reissner-Mindlin kinematics. To our knowledge one article covers this subject for plates [START_REF] Yosibash | Dynamic response of various von Kármán non-linear plate models and their 3-d counterparts[END_REF], where the dynamic response of simple, full and fully threedimensional plate equations are compared in particular cases. The paper showed clearly that the simple von Kármán model is a very good approximation when the thickness is smaller than 1/20 of other typical dimensions (such as a side length) [START_REF] Yosibash | Dynamic response of various von Kármán non-linear plate models and their 3-d counterparts[END_REF]. However the number of dynamical solutions simulated is scarce, and the high-frequency range was not tested so that other simulations are still needed in order to have a better picture of the validity limit of the model. Note that for shells, where a variety of simplified models are also available, a comparative study is provided in [START_REF] Amabili | A comparison of shell theories for large-amplitude vibrations of circular cylindrical shells: Lagrangian approach[END_REF], but only in the vicinity of the first eigenfrequencies.

Numerical simulation approaches to the study of plate vibrations have taken on a variety of different forms. Some have involved nonlinear coupled modal descriptions, using a Galerkin-based formulation with eigenmodes as a functional basis [START_REF] Sridhar | Non-linear resonances in the forced responses of plates, part I: symmetric responses of circular plates[END_REF][START_REF] Touzé | Asymmetric non-linear forced vibrations of free-edge circular plates, part I: theory[END_REF][START_REF] Ducceschi | Nonlinear dynamics of rectangular plates: investigation of modal interaction in free and forced vibrations[END_REF] or an ad-hoc functional basis having good properties [START_REF] Amabili | Theory and experiments for large-amplitude vibrations of rectangular plates with geometric imperfections[END_REF]. Finite difference methods have been successfully used in [START_REF] Bilbao | A family of conservative finite difference schemes for the dynamical von Kármán plate equations[END_REF]. Finite element approaches have also been used (see e.g. [START_REF] Ribeiro | Nonlinear vibration of plates by the hierarchical finite element and continuation methods[END_REF][START_REF] Ribeiro | Non-linear forced vibrations of thin/thick beams and plates by the finite element and shooting methods[END_REF][START_REF] Boumediene | Nonlinear forced vibration of damped plates by an asymptotic numerical method[END_REF][START_REF] Touzé | Direct finite element computation of non-linear modal coupling coefficients for reduced-order shell models[END_REF]), as well as spectral and pseudo-spectral methods [START_REF] Kirby | Solution of von Kármán dynamic non-linear plate equations using a pseudo-spectral method[END_REF][START_REF] Yosibash | Dynamic response of various von Kármán non-linear plate models and their 3-d counterparts[END_REF][START_REF] Düring | Weak turbulence for a vibrating plate: Can one hear a Kolmogorov spectrum?[END_REF][START_REF] Yokoyama | Weak and strong wave turbulence spectra for elastic thin plate[END_REF][START_REF] Miquel | Transition from wave turbulence to dynamical crumpling in vibrated elastic plates[END_REF]. In most cases, the spatial part of the problem is solved first, and then the solution advanced in time through a standard time discretization procedure, such as, e.g., leap frog, or one of the Runge Kutta families of schemes [START_REF] Awrejcewicz | Spatio-temporal chaos and solitons exhibited by von Kármán model[END_REF][START_REF] Kirby | Solution of von Kármán dynamic non-linear plate equations using a pseudo-spectral method[END_REF][START_REF] Düring | Weak turbulence for a vibrating plate: Can one hear a Kolmogorov spectrum?[END_REF].

In parallel, the last two decades have seen an increasing interest in defining conservative schemes, i.e. time-stepping methods where the conservation of invariants such as, e.g., the discrete energy is an in-built property [START_REF] Simo | Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics[END_REF][START_REF] Simo | A new energy and momentum conserving algorithm for the non-linear dynamics of shells[END_REF][START_REF] Hairer | Geometric Numerical Integration: Structurepreserving Algorithms for Ordinary Differential Equations[END_REF]. General procedures for Hamiltonian wave equations by McLachlan [START_REF] Mclachlan | Symplectic Integration of Hamiltonian Wave Equations[END_REF]. These methods are of great utility as a means of ensuring stability even under highly nonlinear conditions, and particularly in examining the transition to the wave turbulent regime, and indeed spectral characteristics of the fully wave turbulent regime. For the simple von Kármán plate model, a family of energy-conserving schemes has been derived in [START_REF] Bilbao | A family of conservative finite difference schemes for the dynamical von Kármán plate equations[END_REF].

The goal here is to introduce conservative schemes for the full von Kármán plate model where in-plane inertia is retained, in a comparative manner to the developments for the simplified system reported in [START_REF] Bilbao | A family of conservative finite difference schemes for the dynamical von Kármán plate equations[END_REF]. This paper is organized as follows: In Section II., the full von Kármán system is presented, in the case of a rectangular plate under periodic boundary conditions, accompanied by a brief energy analysis. Semi-discrete parameterized families of time stepping methods, with an inherent energy conservation property are introduced in Section III.; such methods are suitable for full discretization by a variety of means. To this end, in Section IV., several examples of full discrete methods, including both simple finite difference schemes, and Fourier spectral discretisation, are discussed accompanied by stability analysis, following from an energy conservation property in the fully discrete case. Such methods are distinct from those arrived at through separate spatial and temporal discretisation. Simulation results are presented in Section V. For the sake of completeness, the addition of rotatory inertia and loss terms is covered in the Appendix.

II. MODEL: FULL VON K ÁRM ÁN SYSTEM

The full von Kármán system, describing transverse and in-plane large amplitude vibration of a flat plate, has been presented by various authors [START_REF] Nayfeh | Linear and Nonlinear Structural Mechanics[END_REF][START_REF] Lukasiewicz | Local loads in plates an shells[END_REF][START_REF] Reddy | Mechanics of Laminated Composite Plates[END_REF][START_REF] Chia | Nonlinear Analysis of Plates[END_REF]. In nondimensional form, it may be written as

ü -div N = 0 (2.1a) ẅ -div div M -div (N∇w) = 0 (2.1b) M = -φ(∇∇w) (2.1c) N = φ(ǫ) (2.1d) ǫ = 1 2 ∇u + ∇u T + ∇w ⊗ ∇w (2.1e)
where here, w(x, y, t) is transverse displacement and u(x, y, t) = [u(x, y, t), v(x, y, t)] T is the in-plane displacement vector; both are functions of spatial coordinates x and y and time t. Dots indicate time differentation, and the operator ∇ is the gradient in 2D. The problem is assumed defined for t ∈ R + , and over a domain (x, y) ∈ Ω ⊂ R 2 . M and N are the bending moment and membrane force tensors, respectively. The independent and dependent variables may be related to their dimensional counterparts (with overbars) through the scalings

(x, ȳ) = r 0 (x, y) ū = r 0 u w = r 0 w t = r 0 ρH A t M = AHM N = AN (2.
2) Here, ρ is density and H is thickness, both assumed constant here. The membrane stiffness A is defined as A = EH/(1 -ν 2 ) in terms of H, Young's modulus E, and Poisson's ratio ν. The constant r 0 is defined as r 0 = H/ √ 12. The tensors M and N are defined in terms of a general 2×2 tensor mapping F → φ(F):

φ(F) = νtr (F) 1 + (1 -ν) F (2.3)
where 1 indicates the 2×2 identity matrix and tr (•) is the trace operation. See, e.g., [START_REF] Tataru | On the Cauchy problem for the full von Kármán system[END_REF].

The tensor gradient ∇ and product ⊗ are defined, in terms of functions f

= [f x , f y ] T and g = [g x , g y ] T as ∇f = ∂fx ∂x ∂fx ∂y ∂fy ∂x ∂fy ∂y f ⊗ g = f x g x f x g y f y g x f y g y (2.4)
ǫ is the mid-plane strain tensor, which is symmetric. Note in particular that div div M = -∆∆w, leading to the usual linear term involving the biharmonic operator in thin linear (Kirchhoff) models of plate vibration [START_REF] Graff | Wave Motion in Elastic Solids[END_REF].

The assumptions underlying this model are that:

• a Kirchhoff-Love kinematic is assumed, where any normal to the plate mid-surface before deformation remains normal to the deformed mid-surface. The transverse shear stresses are thus neglected;

• normal stresses along the transverse directions are neglected;

• the von Kármán-like strain-displacement law of (2.1d) is used, by neglecting nonlinear terms of higher order in the plane part of the Green-Lagrange strain tensor;

• the material is linear, homogeneous and isotropic;

• rotatory inertia terms are neglected.

Under further assumptions, namely that in-plane inertia terms may be neglected, it is possible to arrive at a simplified form of the von Kármán system, written in terms of transverse displacement w alone and a scalar function (often referred to as the Airy stress function). Conservative methods for this simplified system have been presented previously by one of the authors [START_REF] Bilbao | A family of conservative finite difference schemes for the dynamical von Kármán plate equations[END_REF]. The rotatory inertia term may be reintroduced here with little difficulty, but is neglected here for simplicity-see the Appendix for a discussion of the addition of this effect along with loss terms.

A. Boundary Conditions

In this article, concerned with general strategies for the design of time-dependent update schemes for system (2.1), the system is assumed defined over a closed square region of nondimensional side-length L:

Ω = {(x, y) ∈ [0, L] × [0, L]} (2.5)
and the boundary conditions are assumed to be of periodic type, by imposing that all displacements, slopes, force resultants and moments at the x = 0 (y = 0) edge are the same than those at the edge x = L (y = L). Such periodic boundary conditions are often used in studies of wave turbulence [START_REF] Düring | Weak turbulence for a vibrating plate: Can one hear a Kolmogorov spectrum?[END_REF][START_REF] Humbert | Wave turbulence in vibrating plates: The effect of damping[END_REF][START_REF] Yokoyama | Weak and strong wave turbulence spectra for elastic thin plate[END_REF]. In the discrete setting, over regular grids, such boundary conditions are easily implemented using periodic difference operations, which reduce to circulant matrix forms when such schemes are vectorized.

B. Inner Products

For two functions F(x, y, t) and G(x, y, t) representing general scalar/vector/tensor fields, an L 2 spatial inner product over Ω is defined as

F, G = Ω F ⋄ GdΩ, (2.6)
where ⋄ is the standard product in the case of scalar fields, the standard dot-product in the case of vector fields and F ⋄ G = tr F T G in the case of tensor fields (tr(•) indicates the trace operation and T the transpose operation). The norm of such a field F is defined as

F = F, F 1/2 (2.7)
C. Identities: Time Differentiation

For tensor (or scalar or vector) F and G,

d dt F, G = Ḟ, G + F, Ġ and d dt 1 2 F 2 = F, Ḟ (2.8)
where d/dt indicates a total derivative of a scalar quantity, and where the dot notation indicates partial differentiation with respect to t.

D. Identities: Integration by Parts

For scalar functions f and g defined over the periodic domain Ω, integration by parts holds as follows:

f, ∆g = -∇f, ∇g (2.9a) f, ∆∆g = ∆f, ∆g (2.9b)
Here, ∆ (•) = div (∇•) is the Laplacian operator in 2D.

For vector (resp. scalar) f and tensor (resp. vector) G, integration by parts may be written as

f , div G = -∇f , G (2.10)

E. Identities: Tensors

For vectors f and h and tensor G, it is true that

f ⋄ (Gh) = (f ⊗ h) ⋄ G =⇒ f , Gh = f ⊗ h, G (2.11)
For a tensor F and a symmetric tensor G, it is true that

F ⋄ G = 1 2 F + F T ⋄ G =⇒ F, G = 1 2 F + F T , G (2.12)
F. Properties of the mapping φ(F)

The mapping φ(F) preserves symmetry:

F = F T =⇒ φ (F) = φ (F) T (2.13)
This property follows from the definition of the mapping φ(F) from (2.3) in terms of the identity and F itself.

The mapping φ(F) satisfies a non-negativity property:

φ(F) ⋄ F = ν (tr (F)) 2 + (1 -ν) F ⋄ F ≥ 0 when 0 ≤ ν ≤ 1 (2.14)
(Note that in the present application, ν, representing Poisson's ratio, satisfies 0 ≤ ν ≤ 1/2.) This property extends to the inner product as

F, φ(F) = ν tr (F) 2 + (1 -ν) F 2 (2.15)
Using the identity F ⋄ φ (G) = φ (F) ⋄ G, the following time differentiation identity also holds:

Ḟ, φ(F) = d dt 1 2 F, φ(F) (2.16)

G. Energy

Taking the vector inner product of (2.1a) with u gives:

u, ü -u, div N = 0 (2.17a) u, ü + ∇ u, N (2.10) = 0 (2.17b) u, ü + 1 2 ∇ u + ∇ uT , N (2.12)
= 0 (2.17c)

d dt T m + 1 2 ∇ u + ∇ uT , N (2.8) = 0 (2.17d)
where the in-plane kinetic energy T m is defined as

T m = 1 2 u 2 ≥ 0 (2.18)
Similarly, taking the scalar inner product of (2.1b) with ẇ gives

ẇ, ẅ -ẇ, div divM -ẇ, div (N∇w) = 0 (2.19a) ẇ, ẅ -∇∇ ẇ, M + ∇ ẇ, N∇w (2.10) = 0 (2.19b) ẇ, ẅ -∇∇ ẇ, M + 1 2 (∇ ẇ ⊗ ∇w + ∇w ⊗ ∇ ẇ) , N (2.11)(2.12) = 0 (2.19c) d dt T b + ∇∇ ẇ, φ(∇∇w) + 1 2 (∇ ẇ ⊗ ∇w + ∇w ⊗ ∇ ẇ) , N (2.8)(2.1c) = 0 (2.19d)
where, using (2.9), the bending kinetic energy T b is defined as

T b = 1 2 ẇ 2 ≥ 0 (2.20)
Adding (2.17d) and (2.19d) gives

d dt (T m + T b ) + ∇∇ ẇ, φ(∇∇w) + ǫ, φ(ǫ) (2.1d) = 0 (2.21)
Using identity (2.16), an energy balance results, i.e.,

d dt H = 0 (2.22)
where the Hamiltonian or total energy H is defined as

H = T m + T b + V b + V m (2.23)
where the in-plane potential energy V m and bending potential energy V b are defined as

V m = 1 2 ǫ, φ(ǫ) V b = 1 2 ∇∇w, φ(∇∇w) (2.24)
From (2.14), it follows that

H ≥ 0 -→ H(t) = H(0) ≥ 0 for t ≥ 0 (2.25)

III. SEMI DISCRETE APPROXIMATIONS

In this section, a semi-discrete (in time) approximation to the full von Kármán system (2.1) is introduced, so as to show Hamiltonian construction techniques, regardless of the form of spatial discretization. In particular, care must be taken to split the nonlinear terms, as shown below.

Suppose that the problem is discretized in time, so that functions w n = w n (x, y), u n = u n (x, y), are to be calculated at integer time step n. If the time step is k, then such functions are approximations to the solution at times t = nk. All inner product and norm definitions and spatial differentiation identities hold as in the continuous case.

A. Time difference and Averaging Operators

For a scalar, vector or tensor f = f n , the following shift operators and the identity may be defined as

e t+ f n = f n+1 e t-f n = f n-1 1f n = f n (3.1)
Forward, backward and centered difference operators, all approximations to a first time derivative may be defined in terms of unit shifts and the identity as

δ t+ = 1 k (e t+ -1) δ t-= 1 k (1 -e t-) δ t• = 1 2k (e t+ -e t-) (3.2)
A centered approximation to a second time derivative may be defined as

δ tt = δ t+ δ t-= 1 k 2 (e t+ -2 + e t-) (3.3)
Averaging operators may be defined as

µ t+ = 1 2 (e t+ + 1) µ t-= 1 2 (1 + e t-) µ t• = 1 2 (e t+ + e t-) µ tt = µ t+ µ t- (3.4)

B. Identities

The following identities hold, for any

f = f n δ t• f, δ tt f = δ t+ 1 2 δ t-f 2 (3.5a) δ t• f, f = δ t+ 1 2 f, e t-f (3.5b)
Concerning the mapping φ, the following identities hold, for tensors F and G:

δ t• F, φ(F) = δ t+ 1 2 F, e t-φ(F) (3.6a) δ t+ F, φ(G) = δ t+ F, µ t+ φ(G) + µ t+ F, δ t+ φ(G) (3.6b) F, φ(e t-F) = µ t-F, φ(µ t-F) - k 2 4 δ t-F, φ(δ t-F) (3.6c) C. Splitting the Tensor N
The membrane strain tensor ǫ, now also semi-discrete in time, with ǫ = ǫ n , can be decomposed classically as

ǫ = ω + γ with ω = 1 2 ∇u + ∇u T γ = 1 2 ∇w ⊗ ∇w (3.7)
where ω and γ denote respectively the linear and the quadratic part of the strain tensor. Note that ω is proportional to the in-plane displacement vector u only, whereas the quadratic part involves the transverse displacement w only. Accordingly, the in-plane force tensor N can de decomposed as:

N = N u + N w (3.8)
where N u groups the linear terms in the in-plane displacement and N w the quadratic terms in w.

Analyzing the linear and nonlinear terms appearing in the full von Kármán model, one can note the following features:

• Transverse and in-plane motions are linearly decoupled, and nonlinearly coupled.

• The nonlinear coupling for the in-plane motions is quadratic and only due to the term div N w . Hence in-plane motions are proportional to quadratic terms in the transverse displacement only, formally of the form w 2 .

• For the transverse motions the nonlinear coupling terms involve a quadratic coupling of the form uw (coming from the term div (N u ∇w)), and a cubic term of the form w 3 (provided by the term div (N w ∇w)).

D. Parameterized Approximations and Identities

For a scalar w(x, y, t), a three-level symmetric discrete approximation w n ξ to w(x, y, t) at t = nk, parameterized by the real number ξ is

w n ξ = ξw n + (1 -ξ)µ t• w n (3.9)
For the quantity γ = 1 2 (∇w ⊗ ∇w), a parameterized approximation, again over three time levels, is

γ ξ• = ξ 4 (µ t• ∇w ⊗ ∇w + ∇w ⊗ µ t• ∇w) + 1 -ξ 2 µ tt (∇w ⊗ ∇w) (3.10)
A parameterized (backwards) approximation over two time levels is

γ ξ-= ξ 4 (e t-∇w ⊗ ∇w + ∇w ⊗ e t-∇w) + 1 -ξ 2 µ t-(∇w ⊗ ∇w) (3.11)
The following identities hold:

γ ξ• = µ t+ γ ξ- (3.12a) δ t+ γ ξ-= 1 2 (δ t• ∇w ⊗ ∇w ξ + ∇w ξ ⊗ δ t• ∇w) (3.12b)
A backwards approximation to ω is given by

ω -= µ t-ω (3.13)

E. Semi-discrete Schemes for the von Kármán System

Consider two step schemes of the following form, which follow from (2.1a) and (2.1b):

δ tt u -div Nu + Nw = 0 (3.14a) δ tt w -div div M -div Nu ∇ w + Nw ∇ w = 0 (3.14b)
Here, Nu = Nn u and Nu = Nn u represent approximations (possibly distinct) to N u at t = nk, and similarly, Nw = Nn w and Nw = Nn w approximations (possibly distinct) to N w at t = nk. w = wn and w = wn are approximations to w, again at t = nk. The forms of these approximations will be specified subsequently. M = M n is defined as in (2.1c).

Similarly to the fully continuous case, one may take inner products of (3.14a) and (3.14b) with δ t• u and δ t• w to get

δ t• u, δ tt u -δ t• u, div Nu + Nw = 0 (3.15a) δ t• w, δ tt w -δ t• w, div div M -δ t• w, div Nu ∇ w + Nw ∇ w = 0 (3.15b)
Using identities (3.5a) and (3.5b), and integration by parts for the second term in (3.15b), and adding the two equations, one may arrive at

δ t+ (T m + T b + V b )-δ t• u, div Nu + Nw -δ t• w, div Nu ∇ w + Nw ∇ w = 0 (3.16)
where

T m = 1 2 δ t-u 2 T b = 1 2 δ t-w 2 V b = 1 2 ∇∇w, e t-φ(∇∇w) (3.17)
Furthermore, using integration by parts (2.10), one may continue to

δ t+ (T m + T b + V b ) + δ t• ∇u, Nu + Nw + δ t• ∇w, Nu ∇ w + Nw ∇ w = 0 (3.18)
Restricting Nu , Nw , Nu and Nw to be symmetric (as is the case for N itself), one may go further, using identities (2.11) and (2.12), and definition (3.7), to

δ t+ (T m + T b + V b ) + δ t• ω, Nu + Nw + 1 2 δ t• ∇w ⊗ ∇ w + ∇ w ⊗ δ t• ∇w, Nu (3.19) + 1 2 δ t• ∇w ⊗ ∇ w + ∇ w ⊗ δ t• ∇w, Nw = 0
At this point definitions of the approximations Nu , Nu , Nw , Nw , w and w are necessary. The goal is to find such choices such the remaining terms in the energy balance above coalesce into a single time difference of a scalar quantity, V m .

F. A Conservative Family of Schemes

There are three sets of terms in (3.19) which need to be examined separately:

• δ t• ω, Nu • δ t• ω, Nw + 1 2 δ t• ∇w ⊗ ∇ w + ∇ w ⊗ δ t• ∇w, Nu • 1 2 δ t• ∇w ⊗ ∇ w + ∇ w ⊗ δ t• ∇w, Nw
Consider the following choices of approximation:

Nu =φ(ω) (3.20a) Nu =µ tt φ(ω) (3.20b) w =w α (3.20c) w =w β (3.20d) Nw =φ(γ α• ) (3.20e) Nw =φ(γ β• ) (3.20f)
The latter four approximations are parameterized by the real numbers α and β, as per the discussion in Section IIID.

Concerning the first term above, and using the definition (3.20a) and identity (3.6a):

δ t• ω, Nu = δ t• ω, φ(ω) = δ t+ 1 2 ω, e t-φ(ω) (3.21)
Concerning the second term above, using definitions (3.20b), (3.20d) and (3.20f), as well as identity (3.12b), it may be rewritten as

δ t• ω, Nw + 1 2 δ t• ∇w ⊗ ∇ w + ∇ w ⊗ δ t• ∇w, Nu = δ t• ω, φ(γ β• ) + δ t+ γ β-, µ tt φ(ω) (3.22
) or, using definition (3.13) and identity (3.12a), as

δ t+ ω -, µ t+ φ(γ β-) + δ t+ γ β-, µ t+ φ(ω -) (3.23)
Finally, using identity (3.6b), this may be rewritten as

δ t+ ω -, µ t+ φ(γ β-) + µ t+ ω -, δ t+ φ(γ β-) = δ t+ ω -, φ(γ β-) (3.24)
The third term may be written, using the definitions (3.20e) and (3.20c), and identity (3.6b), as

δ t+ γ α-, φ(γ α• ) = δ t+ γ α-, µ t+ φ(γ α-) = δ t+ 1 2 γ α-, φ(γ α-) (3.25)
This family of schemes, dependent on the parameters α and β, is thus conservative:

δ t+ H = 0 H = T m + T b + V b + V m (3.26)
where V m given by

V m = 1 2 ω, e t-φ(ω) + ω -, φ(γ β-) + 1 2 γ α-, φ(γ α-) (3.27)
A further manipulation leads to great simplification in subsequent analysis. The expressions for V m in (3.17) and the first term in V m in (3.27) above may be written, using identity (3.6c) as

1 2 ω, e t-φ(ω) = 1 2 ω -, φ(ω -) - k 2 8 δ t-ω, φ(δ t-ω) (3.28) 1 2 ∇∇w, e t-φ(∇∇w) = 1 2 µ t-∇∇w, µ t-φ(∇∇w) - k 2 8 δ t-∇∇w, δ t-φ(∇∇w)
The conserved energy H may then be rewritten as

H = T ′ m + T ′ b + V ′ b + V ′ m (3.29)
where

T ′ m = 1 2 δ t-u 2 - k 2 8 δ t-ω, φ(δ t-ω) (3.30a) T ′ b = 1 2 δ t-w 2 - k 2 8 δ t-∇∇w, δ t-φ(∇∇w) (3.30b) V ′ b = 1 2 µ t-∇∇w, µ t-φ(∇∇w) (3.30c) V ′ m = 1 2 ω -, φ(ω -) + ω -, φ(γ β-) + 1 2 γ α-, φ(γ α-) (3.30d)
In this form, it is particularly simple to arrive at bounds on solution size in the fully discrete case. See Section IV D.

G. Special Case: α = β

When α = β, the energy term V ′ m may be simplified to

V ′ m = 1 2 ω -+ γ α-, φ(ω -+ γ α-) = 1 2 ǫ α-, φ(ǫ α-) ≥ 0 where ǫ α-= ω -+ γ α- (3.31
) Here, contributions to the total energy due to nonlinear effects are assured to be nonnegative (note also that V ′ m ≥ 0), leading to great simplification in the subsequent analysis of stability in fully discrete schemes. In essence, stability analysis reduces to the analysis of the linear problem-i.e., determining conditions for non-negativity of the kinetic energy terms T ′ m and T ′ b . There are more general conditions on the parameters α and β such that V ′ m ≥ 0, but given the practical importance of the case α = β = 1 (see below), such conditions will not be explored further here.

H. Special Case

: α = β = 1
The scheme (3.14) is necessarily implicit. When either α = 1 or β = 1, however, at a given time step, the determination of unknown values in terms of previously computed values requires the solution of nonlinear equations, leading to new concerns regarding existence and uniqueness.

When α = β = 1, however, the unknowns appear linearly, and thus may be solved for uniquely. Expanding the various time difference and averaging operators in (3.14) explicitly at time step n leads to update equations of the form

A (u n , w n ) u n+1 w n+1 = b u n , w n , u n-1 , w n-1 (3.32)
Here, [u n+1 , w n+1 ] T represents the unknown state at the next time step (n + 1), and the operators A and b depend only on previously computed values at the time steps n and n -1. In particular, in the fully discrete case, the operator A may be written as a square matrix dependent on previously computed values, and a linear system solution will be required for the update.

IV. FULLY DISCRETE SCHEMES

The family of schemes described in the previous section is semi-discrete in time, and possesses an energy conservation property. Spatial discretization is as yet unspecified, and may be carried out in a variety of ways-finite element methods are often used, but pseudospectral methods have also been employed for a variant of this system [START_REF] Kirby | Solution of von Kármán dynamic non-linear plate equations using a pseudo-spectral method[END_REF][START_REF] Yosibash | Dynamic response of various von Kármán non-linear plate models and their 3-d counterparts[END_REF].

For illustrative purposes, in this section, a simple discretization over a regular grid will be employed. The grid functions w and u,

w n l,m u n l,m (l, m) ∈ Ω Q = {(l, m)|0 ≤ l, m ≤ Q -1} (4.1)
are approximations to w n (x, y) and u n (x, y) at x = lh and y = mh over the periodic domain Ω Q corresponding to a square of side length L, where h, the grid spacing is defined as h = L/Q for some positive integer Q.

A. Inner Products, Difference and Fourier Operators and Summation by Parts

For any two grid functions F l,m and G l,m (again scalars, vectors, or general rectangular matrices), the l 2 spatial inner product and norm over Ω Q are defined as

F, G = h 2 (l,m)∈ΩQ F l,m ⋄ G l,m F = F, F 1/2 (4.2)
where the ⋄ operation is as defined in (2.6). Consider now a generic difference operator D = [D x D y ] approximating the gradient ∇ over the grid Ω Q . Similarly to the case of integration by parts, such an operator will satisfy a summation by parts identity with an associated operator, D † , representing an approximation to the divergence operation:

f, D † g = -Df, g (4.3)
There are obviously many choices of such an operator; the simplest possible choices are

D = D + = [D x+ D y+ ] T and D = D -= [D x-D y-]
T , representing forward and backward difference operations. In terms of (say) a vector grid function f n l,m , the operations D x+ and D x-are defined as:

D x+ f n l,m = 1 h f n l+1,m -f n l,m D x-f n l,m = 1 h f n l,m -f n l-1,m (4.4)
where spatial indeces l are taken modulo Q. A similar definition holds for D y+ and D y-.

The associated divergence operations are

D † + = D T -and D † -= D T + .
Another useful choice of approximation (used often in studies of wave turbulence) is the Fourier spectral gradient D = D spec = [D spec,x D spec,y ] T . When applied, e.g., to a vector function f n l,m , the operator D spec,x may be written as

(D spec,x f ) n l,m = π N -1 l ′ =0,l ′ =l (-1) l-l ′ cot π (l -l ′ ) Q f n l,m (4.5)
with a similar definition holding for D spec,y . Such operators may be efficiently applied using the Fast Fourier Transform [START_REF] Trefethen | Spectral Methods in Matlab[END_REF].

B. Operator Bounds

For a given difference approximation acting over a grid function defined over Ω Q , it is assumed here that the operation is bounded. For a gradient approximation D applied to a vector or scalar grid function f n l,m , for instance, one has

Df ≤ D max f and 1 2 Df + (Df ) T ≤ D max f (4.6)
For the simple forward and backward difference operations D + and D -, for example, one has D ±,max = 2 √ 2/h. For the spectral operator D spec , the bound exhibits a slightly more complex dependence on h, but in the limit of small time steps, one has approximately D spec,max = √ 2π/h. The associated divergence operators D † are subject to the same bound, i.e., D † max = D max .

C. Discrete Method and Conserved Energy

One may proceed directly from the semi-discrete system (3.14) to a fully discrete method

δ tt u -D † a
Nu + Nw = 0 (4.7a)

δ tt w -D † b D † c M -D † a Nu D a w + Nw D a w = 0 (4.7b)
where here, all variables are assumed defined over Ω Q , and where the various approximations to the gradient operation D a , D b and D c (possibly distinct) have been employed. M is approximated as

M = -φ(D c D b w) (4.8)
and the tensors Nu , Nw , Nu and Nw and the grid functions w and w as in (3.20). The quantity ω is as defined in (3.7), and γ ξ• and γ ξ-are as defined in (3.10) and (3.11), with the operator D a replacing ∇.

Consider the scheme (4.7) with α = β. The energy analysis of the semi-discrete case follows through as before, where now, the various components are defined as

T ′ m = 1 2 δ t-u 2 - k 2 8 δ t-ω, φ(δ t-ω) (4.9a) T ′ b = 1 2 δ t-w 2 - k 2 8 δ t-D c D b w, δ t-φ(D c D b w) (4.9b) V ′ b = 1 2 µ t-D c D b w, µ t-φ(D c D b w) (4.9c) V ′ m = 1 2 ǫ α-, φ(ǫ α-) (4.9d) D. Non-negativity Conditions
Note that by (2.14), V ′ m ≥ 0, and V ′ b ≥ 0. What remains is to bound the terms T ′ m and T ′ b away from zero. To this end, note that from (2.15), one has

δ t-ω, φ(δ t-ω) = ν tr(δ t-ω) 2 + (1 -ν) δ t-ω 2 (4.10a) = ν D † a (δ t-u) 2 + (1 -ν) δ t- 1 2 D a u + (D a u) T 2 (4.10b) ≤ D 2 a,max δ t-u 2 (4.10c)
and

δ t-D c D b w, δ t-φ(D c D b w) = ν tr (δ t-D c D b w) 2 + (1 -ν) δ t-D c D b w 2 (4.11a) ≤ D 2 c,max D 2 b,max δ t-w 2 (4.11b)
Thus under the conditions

k ≤ 2 D a,max k ≤ 2 D b,max D c,max (4.12)
the terms T ′ m and T ′ b are non-negative. For simple difference operations, these bounds are equivalent to those arrived at through von Neumann analysis of the linear part of the von Kármán system alone. For example, if

D a = D + , D b = D + and D c = D -, the conditions above reduce to k ≤ h/ √ 2 and k ≤ h 2 /4
, for a grid spacing h. For the spectral scheme, with

D a = D b = D c = D spec , then k ≤ √ 2h/π and k ≤ h 2 /π 2 .

E. Summary

In sum, the method proposed here in (3.14), defined over a discrete periodic domain is energy conserving, and dependent on two free parameters α and β. Under the choices α = β = 1, the discrete energy function can be shown to be non-negative under simple conditions on the grid spacing in terms of the time step, as for the underlying linear problem. The property is independent of the particular type of spatial discretization, and both simple finite difference and Fourier based methods have been described here.

For readers wishing to examine the implementation of such an algorithm, Matlab code is available as supplementary material. See reference [52] below.

V. SIMULATION RESULTS

Simple examples of the use of the fully discrete method described above are presented in this section. As a useful test configuration, consider an initialisation, over the periodic square domain Ω Q with the first sinusoidal mode shape in the transverse displacement:

w 0 l,m = w 1 l,m = a cos (2πl/Q) cos (2πm/Q) u 0 l,m = u 1 l,m = 0 (5.1)
where here, a is a nondimensional amplitude. (Recall that amplitudes have been scaled with respect to the quantity H/ √ 12, where H is the physical plate thickness.) The nondimensional side length L is chosen as L = 20 in all cases, and Poisson's ratio ν is set to ν = 0.3.

A. Time Evolution of Plate Profile under Different Initial Condition Amplitudes

As a first example, consider the use of the scheme (3.14), with α = β = 1, and using a simple finite difference discretization as given in Section IV A, and employing varying initial condition amplitudes a.

Under very small amplitude conditions (here, a = 0.01), the plate undergoes transverse vibration in its first linear mode, at nondimensional frequency 4π/L 2 . See Figure 1. As the initial amplitude is increased, particularly beyond approximately a = 1, the gross frequency of vibration of the plate increases, and as a is increased further, beyond about a = 3, one begins to see transfer of energy to other modes, and the vibration is in general aperiodic. See Figure 2. See also Figure 3, showing the displacement of the corner value of the displacement w n 0,0 as a function of time, for various different initial condition amplitudes. The time step here is chosen as k = 1/50, and the simple finite difference discretization described in Section IV A is employed, with h chosen to satisfy the stability bounds (4.12) as close to equality as possible. In this case, the grid size is given by Q = 70.

B. Fourier vs. Finite Difference Approximations

Particularly at high initial condition amplitudes, calculated solutions using a simple spatial discretization method such as finite differences are subject to substantial drift. See Figure 4, showing a comparison between a high accuracy solution, calculated using a small time step, and one calculated using a coarse spatial grid with a larger time step. As might be expected, the Fourier discretization does a much better job, and solutions calculated using a large time step over a coarse spatial grid match well with those from an accurate finite difference scheme. See Figure 4. Particularly when it is of interest to examine delicate phenomena such as e.g., mode bifurcations under forced conditions, such a method (which is similar to modal methods used previously [START_REF] Ducceschi | Nonlinear dynamics of rectangular plates: investigation of modal interaction in free and forced vibrations[END_REF]) will be preferred. Note, however, that in a conservative scheme which can be written in an update form such as (3.32), the difference matrices will no longer be sparse, so that the computational effort associated with a time update will be considerably larger than in the case of locally defined difference operators.

C. Energy Conservation

In this section, the energy conservation property of the schemes presented in Section IV. is demonstrated, for various choices of the initial condition amplitude a. Calculated values of the total conserved energy H are given in Table I, for the case of a finite difference discretization, and in Table II for a Fourier discretization.

In both cases, even at high amplitudes, energy is conserved to many places-in this case, in double precision floating point, near the level of machine accuracy. For more comments on accuracy in finite precision, see Section VI.

Also of interest is the partition of the energy into the various potential and kinetic components as defined in (4.9).

VI. CONCLUDING REMARKS

This article has explored a family of energy conserving numerical methods for the full von Kármán system, in both semi-discrete and fully discrete forms, and, more generally, with construction techniques for implicit numerical time-stepping methods for nonlinear systems. Though, for simplicity, only the doubly periodic domain has been discussed here, along with relatively simple spatial discretization strategies, it is anticipated that nontrivial boundary conditions and more elaborate discretization methods should fit into the framework described here. There are, however, various issues which remain unresolved.

The two-parameter family of schemes given in (3.14), along with the split-tensor time discretizations given in (3.20) leads to a semi-discrete update with a conserved energy. Under appropriate choices of the spatial discretizations, the conservation of energy then follows through to the fully discrete case. The determination of numerical stability reduces to conditions under which the conserved energy is a non-negative function of the state. Under particular choices of the scheme parameters (α = β), such conditions are easily arrived at, and reduce to stability conditions for the underlying linear systemthat is, the additional potential energy associated with the nonlinearity is non-negative. When α = β, however, despite the fact that the scheme is conservative, such nonnegativity conditions are not immediately forthcoming. A related issue is the existence and uniqueness of numerical solutions-these properties can be ensured when α = β = 1, as can easily be deduced from the update form (3.32), but otherwise are not immediately apparent.

At the level of the von Kármán system itself, due to the disparity between wave speeds, even in the linear case, between the longitudinal and transverse motions, the variables u should be treated distinctly (using either distinct time steps or distinct grids) in order to achieve the best numerical behaviour (i.e., the least anomalous numerical dispersion). This issue has not been addressed in this article, but it is hoped that the construction techniques employed here can be extended to handle such cases.

Finally, all the schemes presented here are implicit, including the useful special case of α = β = 1, and will require the solution of linear systems, which form the largest part of the computational load. Under small amplitude vibration conditions, the linear systems to be solved are diagonally dominant, and relatively simple and efficient methods (such as, e.g., Jacobi methods [START_REF] Strikwerda | Finite Difference Schemes and Partial Differential Equations[END_REF]) may be used; at higher amplitudes, however, the diagonal dominance property is lost, and one must revert to more general (and computationally costly) methods, such as, e.g., conjugate gradient techniques. Simple conditions on the range of vibration amplitudes which guarantee diagonal dominance would be of great interest. On top of this, though numerical energy is conserved in infinite precision, finite wordlength effects lead to a deterioration of the energy conservation property. For very large vibration amplitudes, one should expect that the linear systems to be solved will be poorly-conditioned, and thus one will see greater variation in the numerical energy, leading, ultimately, to instability due to truncation effects. In the present case, however, it appears that such vibration amplitudes greatly exceed the conditions for validity of the von Kármán model itself, and thus the family of schemes presented here should be safe to use under realistic conditions.

VII. ACKNOWLEDGMENTS

This work was supported by the European Research Council, under grant StG-2011-279068-NESS.

APPENDIX: ROTATORY INERTIA AND LOSS

If effects of rotatory inertia and loss are included, the dynamic equations (2.1a) and (2.1b), in dimensionless form, may be augmented to ü -div N + σ u u = 0 (7.1a) ẅ -∆ ẅ -div div M -div (N∇w) + σ w ẇ = 0 (7.1b) where the terms involving σ u ≥ 0 and σ w ≥ 0 allow for simple loss.

In this case, the energy conservation equation (7.2) becomes

d dt H ′ = -Q (7.2)
where

H ′ = H + 1 2 ∇ ẇ 2 ≥ 0 Q = σ u u 2 + σ w ẇ 2 ≥ 0 (7.
3

)
The system is thus dissipative, so that 0 ≤ H ′ (t) ≤ H ′ (0), for t ≥ 0. Similarly, scheme (3.14) may be augmented as δ tt u -div Nu + Nw + σ u δ t• u = 0 (7.4a) δ tt w -δ tt ∆w -div div M -div Nu ∇ w + Nw ∇ w + σ w δ t• w = 0 (7.4b)

where the centered time difference operator δ t• , as defined in (3.2) has been employed. The semi-discrete conservation equation (3.26) then becomes

δ t+ H ′ = -Q (7.5)
where

H ′ = H + 1 2 δ t-∇w 2 ≥ 0 Q = σ u δ t• u 2 + σ w δ t• w 2 ≥ 0 (7.6)
The semi-discrete system is thus also dissipative, so that 0 ≤ H ′n ≤ H ′0 , for time steps n ≥ 0.

FIG. 1 .

 1 FIG. 1. Snapshots of plate displacement, under an initial condition of the first linear mode shape, with a = 0.01, at times as indicated. The transverse displacement field w is shown along the top row, and one longitudinal displacement (the first component of u) on the bottom row. Light and dark indicate regions of positive and negative displacement, respectively.

FIG. 2 .

 2 FIG. 2. Snapshots of plate displacement, under an initial condition of the first linear mode shape, with a = 5, at times as indicated. The transverse displacement field w is shown along the top row, and one longitudinal displacement (the first component of u) on the bottom row. Light and dark indicate regions of positive and negative displacement, respectively.

FIG. 3 .

 3 FIG.3. Plate response, for a plate with L = 20, and ν = 0.3, under an initial condition of the first linear mode shape, for different amplitudes a, as indicated. The response, taken as w n 0,0 , is plotted against time. The time step is taken as k = 1/50.

FIG. 4 .

 4 FIG.4. Plate response, for a plate with L = 20, and ν = 0.3, under an initial condition of the first linear mode shape, for amplitude a = 5. The response, taken as w n 0,0 , is plotted against time, for three different simulations: Solid line: an accurate response using a finite difference discretization, for k = 1/50. Dotted line: coarse response using a finite difference discretization, with k = 1/10. Dashed line: coarse response using a Fourier discretization, with k = 1/10.

FIG. 5 .

 5 FIG. 5. Partition of numerical energy as a function of time, for a finite difference scheme, under different values of the initial amplitude a, as indicated. The various components of the energy, T ′ b (solid black), V ′ b (solid grey), T ′ m (dashed black) and V ′ m (dashed grey), as given in (4.9) are shown. The total energy H is indicated by a thick black line. The time step is again k = 1/20.

TABLE I .

 I Calculated numerical energy values, for different values of the initial condition amplitude a, and for particular values of the time step n, using a finite difference discretization. Here, the time step is chosen as k = 1/20.

	n	a = 0.01	a = 1	a = 5
	1	1.941538876811134×10 -4 2.093504501191837 1.435261060075711×10 2
	2	1.941538876811129×10 -4 2.093504501191842 1.435261060075714×10 2
	3	1.941538876811128×10 -4 2.093504501191841 1.435261060075717×10 2
	4	1.941538876811132×10 -4 2.093504501191845 1.435261060075721×10 2
	1000	1.941538876811936×10 -4 2.093504501194523 1.435261060076516×10 2

TABLE II .

 II Calculated numerical energy values, for different values of the initial condition amplitude a, and for particular values of the time step n, using a Fourier discretization. Here, the time step is chosen as k = 1/20.

	n	a = 0.01	a = 1	a = 5
	1	1.948149596712886×10 -4 2.100321257296917 1.438201589955981×10 2
	2	1.948149596712884×10 -4 2.100321257296906 1.438201589955980×10 2
	3	1.948149596712904×10 -4 2.100321257296910 1.438201589955976×10 2
	4	1.948149596712887×10 -4 2.100321257296918 1.438201589955981×10 2
	1000	1.948149596712743×10 -4 2.100321257296243 1.438201589955047×10 2