Anthony Surleraux

Romain Lepert

Jean-Philippe Phillipe Pernot

Samuel Bigot
email: bigots@cardiff.ac.uk

Computer-aided Micro-EDM die-sinking tool design optimisation

come   L'archive ouverte pluridisciplinaire

Introduction

Micro and nano manufacturing technologies have been recognized as a promising new source of innovation and growth for decades to come [START_REF] Islam | Nanotechnology innovation system: understanding hidden dynamics of nanoscience fusion trajectories[END_REF]. The microelectronics and IT revolution which benefited from such technologies started many decades ago and was built on silicon-based IC technologies. However, capabilities for micro and nano processing on a wider range of materials, including the structuring of true 3D-forms, is becoming increasingly important for the development of innovative new applications. Developing and improving such capabilities to facilitate their integration in the next generation of factories for high value products is crucial and would enable more flexible and cost-effective manufacturing of multifunctional products made of different materials in practically all manufacturing sectors, including biotechnology (e.g. biosensors, micro fluidics) and ICT (e.g. optoelectronics) [START_REF]IPMMAN roadmap for micro-and nano-manufacturing[END_REF]. A particular issue with the manufacturing of 3D forms at the micro scale is the access to reliable and efficient measurement techniques allowing for appropriate part verification. This, in combination with the relative lack of maturity of some of these micro manufacturing technologies, makes the design process difficult as assessing the manufacturability of new features is far from straightforward. At the same time, with the constant improvement in computer graphic visualization and information management systems, new methods to support design and manufacture are emerging, such as the Digital Enterprise Technology (DET) which was defined as "the collection of systems and methods for the digital modelling of the global product development and realization process, in the context of lifecycle management" [START_REF] Maropoulos | Digital enterprise technology-defining perspectives and research priorities[END_REF]. This promises highly positive impact on product development and realization strategies in particular through the use of computer-aided manufacturing risk mitigation approaches. Method such as DET require the development of accurate manufacturing process modelling, which could then be used for reliability modelling and verification [ 4 ] and consequently for product development risk reduction Thus, for micro manufacturing techniques to be fully integrated in these new digital modelling techniques, it is essential to develop advanced modelling and control of manufacturing chains involving micro manufacturing technologies in order to facilitate their consideration in the design of new innovative products. In addition, new developments in high value products factories should allow for rapidly reconfigurable manufacturing chains. This paper is focusing on the integration of an efficient micro machining simulation technique, for the micro Electro Discharge Manufacturing (µEDM) process, which would support the development of such advanced control and facilitate the reconfigurability of manufacturing chains using the µEDM process and in particular the µEDM die-sinking. From any newly designed features to be produced by µEDM die-sinking, the proposed manufacturing simulation technique would enable designers to perform digital optimisation and manufacturability verification through a digital optimisation of µEDM die-sinking tool electrodes design. This will support the development of a sustainable manufacturing operation as, in addition to the use of an IT based design manufacturability check, the tool design optimisation is expected to enhance the efficiency of the µEDM die-sinking, in particular by reducing the number of EDM tools required thus further reducing waste and energy consumption. The sections two and three of this paper give an introduction to the µEDM process and to the developed µEDM simulation technique. This simulation technique is used to develop a new EDM tool electrode optimisation technique presented in section four.

µEDM Introduction

Electrical Discharge Machining (EDM) is a manufacturing process where material removal is due to successive electrical discharges. Its main interest lies in the possibility of machining any conductive material regardless of hardness. A few variants of EDM are in use but share an underlying principle: the tool and the workpiece (or electrodes) are immerged in a dielectric fluid and submitted to an electrical current. The gap between them is reduced until the dielectric reaches its breakdown voltage. The current is then free to flow from one electrode to the other creating a plasma channel. The energy dissipated partly evaporates the tool and workpiece leading to the formation of craters. Micro-EDM (µEDM) deals with micrometric dimensions. The polarity is usually inverted in order to alleviate the effects of tool wear: as the process goes on, material is removed from the tool, modifying its shape. While of manageable consequences in classic EDM, the influence of tool wear is most notable in its micro counterpart as can be seen in Figure 1. While proven methods exist for µEDM milling to compensate for the tool wear, diesinking µEDM would require the use of multiple tools to obtain the desired geometries. The ability to predict the location and intensity of wear in die-sinking µEDM would enable the design of optimal tools able to compensate for the upcoming wear and reduce the number of electrodes required for a successful machining. The optimisation of tool electrode to compensate for tool wear is used regularly in drilling and in micro EDM milling. When producing blind holes using micro-EDM drilling without wear compensation, because of tool wear, the real depth of the hole will be significantly smaller [START_REF] Pham | MicroEDM-recent developments and research issues[END_REF]. A method to achieve a specific depth in this case is to compensate for wear by constant electrode feeding in the Z-axis [START_REF] Bleys | Real time tool wear compensation in milling EDM[END_REF] working length of the tool using the expected wear ratio. In micro EDM milling, the use of simple shaped electrode (generally cylindrical) makes the wear prediction relatively straightforward [START_REF] Bissacco | Feasibility of wear compensation in micro EDM milling based on discharge counting and discharge population characterization[END_REF]. It uses specific layer by layer machining strategies, such as the Uniform Wear Method [START_REF] Yu | 3D micro-EDM with simple shape electrode[END_REF], which ensures that the machining occurs only on the bottom surface of the electrode and results in a relatively linear wear on the z-axis. Therefore, only feeding the electrode in z-axis can do the wear compensation, which is again equivalent to extending the working length of the tool. In the case of die sinking micro EDM, where more complex shapes are used as electrodes, the compensation is not as straightforward as the way electrode geometries are affected by the wear is difficult to predict. A way forward is to consider the use of simulation. This was the case in a paper [START_REF] Kunieda | Reverse simulation of sinking EDM applicable to large curvatures[END_REF] where the final shape of tool electrode can be predicted in one step using a "reverse simulation" technique. In this simulation, the removal of volume on both workpiece and electrode is achieved using voxels' representations and appears to performed per specific time periods on large areas rather than discharge by discharge. It is not clear if this method can be applied to micro-EDM due to resolution that can be applied to the voxels' representation used.

In this paper, an efficient simulation method using voxels in an octree data structure, which allows for high resolution simulation, was used as a building block for the tool electrode optimisation and is described in the next section.

µEDM Simulation Tool

Simulation Overview

The simulation algorithm used is an updated version of an algorithm developed previously by the authors [START_REF] Bigot | Micro-EDM numerical simulation and experimental Validation[END_REF]. Using 3D models of a tool and a workpiece, this algorithm iteratively simulate the EDM process by removing volumes in a manner described in Figure 2: Overall simulation diagram.Figure 2. The 3D models can be imported from any CAD format such as those emerging from 3D scans of real parts, consequently simulations could be performed on models derived from components measured in the middle of a manufacturing chain, prior to an EDM die sinking step, thus allowing for EDM process re-configurability if required.

The following parameter are fixed by the user before each simulation:

-Machining gap: the distance between the workpiece and the tool below which it is possible for sparks to take place. -Workpiece Crater Volume: the volume of "material" removed from the workpiece by each spark. -Tool Crater Volume: the volume of "material" removed from the tool by each spark. The first step consists in checking if the current depth of the tool is smaller than the objective depth. If this is true, the minimum distance d between the tool and the workpiece is computed as well as the corresponding locations. Should that distance d be smaller than the machining gap, craters are added through volume removal on the electrodes (tool and workpiece).

If d is greater than the machining gap, the tool electrode is moved down by a small increment.

The algorithm then compares the current depth against the objective again and a new iteration begins. Should the current depth meet the targeted one, the process is stopped. The simulation tool makes use of voxels embedded in a voxel octree to describe the geometries (the tool and the workpiece) involved in the process. As a result the modelling is purely volumetric due to the use of voxels while octrees significantly decrease memory usage by representing large volumes as a single and larger voxel as illustrated in Figure 3.

The octree data structure is a tree structure in which each node can have up to eight children (hence the name). A node without any children is called a leaf node. The topmost node (containing all others) is called the root.

Figure 3 : Cross-sectional view of a sphere represented as voxels in an octree tree structure.

Algorithms

Minimum distance search

The minimum distance search makes use of the hierarchical data structure available with octrees. The algorithm starts with the root nodes of each electrode and computes the minimum distance between all the possible couples of their children:

-If both children are leaves, the minimum distance between them can be immediately found.

-

If not, the minimum distance can be bounded between a maximum dmax and a minimum value dmin.

All the couple of nodes that verify the following relation are kept:

Where dsmallestMax is the smallest of all the dmax. All the remaining couples' children are then compared in an identical manner until only couples of leaves remain.

A fast exit condition is used in the case where a couple of nodes have an upper bound dmax that is smaller than the machining gap.

Crater insertion

The crater insertion algorithm also uses the hierarchical nature of octrees. At the relevant location, the intersections between a sphere and the root node's children are computed.

-If a child node doesn't intersect, nothing happens. -If a node is completely inside the sphere, it is deleted. -Nodes on the boundary of the sphere are kept for the next iterations.

The algorithm is then reapplied on the new candidate nodes until a certain resolution is met. A cross-sectional view of a crater can be seen in Figure 4.

The result of a full simulation on a workpiece can be seen in Figure 5.

As mentioned previously, this simulation tool can be used for the design of optimal tool electrodes that can take into account the upcoming wear, thus reducing the number of electrodes required for a successful machining. This can be achieved through the shape optimisation tool described in the next section.

Optimisation tool

Overview

What is proposed in this paper is to develop an iterative simulation process that adjusts the shape of the tool electrode in order to obtain an ultimate shape potentially able to produce in one single step the targeted cavity on the workpiece or its closest producible cavity. Although the tool presented in the previous section can provide 3D simulations, to demonstrate and validate the optimisation concept more efficiently it was decided to perform optimisation using the simulation with 2D profiles. However, the method can be easily adapted to 3D shapes. Thus, considering a targeted cavity profile TWp to be produced on a workpiece, each iteration of the optimisation process consists in five main steps:  Compute the difference or mismatch DWp between the TWp and the previous resulting cavity profile RWp. On the first iteration, the RWp is a virgin workpiece since no machining has been done yet. This means that the first DWp represents the actual overall target area to be removed (or volume in 3D).  Using the mismatch, adapt the virtual cavity profile VWp.  Using the VWp, generate a virtual tool profile VT, which would produce the selected VWp if no wear occurred on the tool.  Run the EDM simulation with VT on a new workpiece (not machined). This provides a new resulting cavity profile RWp to be analysed for the next iteration. The accuracy of the optimisation process is measured with the following metric:

The process is repeated until the accuracy no longer improves more than 0.5%, meaning that the resulting cavity is as close to the targeted workpiece cavity as can be produced with this process. As a consequence, the optimal tool for the desired workpiece is acquired. An example of the expected capabilities of the tool shape optimisation process is given in Figure 7. Elements a) and b) represent respectively the targeted cavity profile TWp to be produced on the workpiece and the initial tool shape (virtual tool profile VT at the first iteration). Elements c and d represent respectively the optimized tool shape (final virtual tool profile VT) after five iterations and the resulting workpiece (final cavity profile RWp). In this example, after four iterations, the final shape of the workpiece (fig. 7.d) is similar to the desired shape (fig. 7.a) with an accuracy of 98.2%.

Virtual cavity profile VWp and tool generation VT

A key aspect of this procedure is the way the mismatches occurring after each iteration are used to propose new tool electrodes which would eventually lead to an optimum. At the end of each iteration, an intuitive reaction would be to attempt to use the area difference DWp between the resulting cavity profile RWp and the targeted profile TWp to directly add "material" to the tool. However, due to the machining gap, this is not a straightforward approach. This is the reason why changes and optimisations are applied on the workpiece, rather than on the tool, thus producing a virtual entity called the virtual cavity profile VWp. The virtual cavity's VWp only purpose is to generate its conjugate tool VT, a tool that would produce the virtual cavity in a simulation without tool wear. This conjugate tool is obtained by applying the machining gap on the surface of VWp, as shown in Figure 8. The machining gap being the minimum distance at which sparks occurs and therefore the closest distance a potential tool can get before wearing out on both sides. On each iteration, the VWp is optimized according to the given mismatch DWp, thus a new tool VT is generated and produces a new workpiece shape RWp.

Assessing differences between resulting and targeted profiles

At the end of each iteration, the difference between resulting and targeted workpiece profiles is assessed by doing a symmetry difference. Workpiece profiles are converted into 2D Boolean matrixes where each element of the matrix represents one pixel. These matrices can take two values (Figure 9): 0 for blank pixels, 1 for colored pixels.

The mismatch matrix is then computed by doing the difference between the matrix equivalents of the targeted workpiece profile A (matrix equivalent of TWp) and the resulting workpiece profile B (matrix equivalent of RWp). This difference is performed following the table of truth depicted on Figure 10. This computation outputs two distinctive areas:  positive area: represents the excess of volume removed.  negative area: represents the lack of volume removed. It is important to make this distinction for corrective purpose. Indeed, these areas pinpoint the locations where the tool was under/over effective, as described in Figure 11. With this kind of information, it is possible to modify the virtual cavity and a fortiori the tool in order to compensate for its lack or excess of material removal. Analytics on the mismatch matrix also give valuable information regarding the advancement of the optimisation process.

Optimisation of the virtual cavity profile

Based on the assessed difference between the targeted workpiece profile and the result of an iteration, the optimisation objective is to make modifications on the virtual cavity profile VWp so that the virtual tool profile VT, generated in the next iteration, reduces the difference. The shape corrective function to apply to the virtual cavity is computed by taking the vertical summation of each column of the mismatch matrix (Figure 12). The resulting discrete function represents the vertical difference between the objective and the current result. Positive values referring to an excess of removal, and negative values referring to a lack of removal. The shape corrective function is then applied onto the surface of the virtual cavity, whose shape should evolve to improve the efficiency of the next tool conjugate generated, as described in Error! Reference source not found.. Looking to compensate the difference from a vertical point of view is simpler for data processing but also relevant in the context of die-sinking EDM, because the movement of the tool is restricted to the vertical axis and thus most interactions happen in this direction. It should also be mentioned that the shape corrective function is altered by the wear ratio in order to take into account the future expected wear on the tool, and thus improve the efficiency of each iteration of shape the optimisation process. Thus in cases of excess of volume to be removed, following approaches developed for drilling and EDM milling, one could suggest that the shape corrective factor obtained for each vertical summation should be altered by the wear ratio.

Figure 13 : Optimisation process on iteration (i)

The wear ratio is simply defined by the ratio between the wear occurring on the tool electrode and the wear occurring on the workpiece, which for the simulation would be equivalent to the ratio between the Tool Crater Volume and the Workpiece Crater Volume.

Experiments and discussion

In order to demonstrate the applicability of the tool shape optimisation method, three different targeted profiles have been used. The following parameters were used for all simulations:

 Machining gap: 5μm  Tool wear ratio : TWR = 0.5 To achieve a tool wear ratio of 0.5, the volume ratio is defined by the parameters of the craters which apply respectively on the workpiece and the tool as follows:

 Workpiece: crater radius = 3.00μm and crater depth: 2.25μm.  Tool: crater radius= 2.25μm and crater depth= 1.50μm. The three different target profiles tested are defined as such:

1. Circle (Figure 14 The figure 14 represents the evolution of the accuracy for each iteration until the process stop. The accuracy is the measure of the overall completeness of the machining in regard to the objective profile. This metric stabilizes after a few iterations for each tests subjects and the process stops when accuracy's progression over one iteration becomes less than 0.5%. It means that the following optimisation iterations will not improve the tool significantly anymore. Therefore, the optimisation process stops and the tool giving the highest accuracy is kept as the optimal tool. As can be seen in Figure 15, in all cases the optimal tool was obtained after 3 to 4 iterations with an accuracy rate of over 96%. This certifies that the optimal tools are very effective given that the remaining mismatch is mostly due to the imperfection of craters in regard to the target geometry as shown in figure 16.c) 17.c) an 18.c). 15 : Optimisation process data for all three target profiles.

Iteration number of optimal tool

Accuracy of optimal workpiece result

Mean duration of optimisation loops

It is important to notice that the cavity of the rectangular target profile is significantly bigger than the others, with the cavity of the triangular profile being the smallest. This means that there will be more material to remove. As a consequence, the duration of each iteration for the optimisation of the rectangular profile are the longest and those of the triangular profile are the shortest. The remaining mismatch on the triangular profile (Figure 16.c) peaks at 5.75μm around the tip of the isosceles triangle, questioning whether the sharpness objective of the profile is achievable. This highlights the fact that not all target profiles are achievable with the die-sinking EDM process, a fortiori using a single tool. The final mismatch can be considered as a machinability check metric that indicates to designers if the resulting workpiece will match their dimension tolerances. Consequently, this can be used as a designer's tool to check if the EDM process is adequate for the production of specific new features or, as mentioned previously, in the middle of a manufacturing chain using scans of produced components to check if a re-configurability of the EDM process is required and if it could lead to the desired features' tolerances. In order to test the limits of the tool shape optimisation method, the same three target profiles underwent tool optimisation method with three different machining gaps:

 Machining gap: 5μm, 20μm and 40μm Figures 19, 20 and 21 represent evolution of accuracy during the optimisation process of all three profiles for each machining gap. What is worth highlighting here is that, except for the circular profile, the machining gap significantly worsens the accuracy of the optimal tool as it gets bigger. This is due to the more pointy shapes of the triangular and rectangular profiles. These shapes are more difficult to acquire when the machining gap increase because the vertical compensation applied on the cavity tool is then nullified by the large machining gap as the conjugate tool is generated. As shown in Figure 22, a significant compensation applied on the virtual cavity does not affect the tool much in return. Figure 22 : virtual cavity and respective tool on different iterations with a machining gap of 40μm.

Conclusion and Future works

This paper presented a new method for computer-aided optimisations of micro EDM die sinking tools, which can be used by product designers for digital optimisation and manufacturability verification of micro features to be produced by micro EDM. Its applicability to simple shaped cavities was demonstrated. With further developments, this optimisation tool should enable rapid and cost effective manufacturability checks at product design stages or for manufacturing chains reconfigurability when the micro-EDM die sinking process needs to be considered for the production of new structures. This will be the focus of future experiments, with real machining tests using proposed optimum EDM tools. Further improvements for the optimisation process are also being considered. In particular, to achieve a more efficient optimisation process it could be useful to add other stopping criterion, such as with the achievement of specific tolerances on the end product features. Also, adapting the shape corrective function in cases where the machining gap is a handicap would be important to ensure the achievement of the optimal tool. Finally, the proposed framework enables the optimisation of the shape of a single tool, future works should also concentrate on the optimisation of a number of tools used in a sequence when a cavity cannot be achieved with as single one or when roughing and finishing strategies should be considered.

Figure 1 :

 1 Figure 1: Tool wear influence in micro-EDM. a) Desired result. b) Obtained result.

Figure 2 :

 2 Figure 2: Overall simulation diagram.

Figure 4 :

 4 Figure 4 : Cross-sectional view of a crater at a resolution of one voxel per 125 nm.

Figure 5 :

 5 Figure 5 : Voxel workpiece after simulation (left) and details of the craters (right).

Figure 6 :

 6 Figure 6 : Iterative optimisation process.

Figure 7 :

 7 Figure 7 : Tool shape optimisation capabilities

Figure 8 :

 8 Figure 8 : Virtual cavities VWp and their respective conjugate tool VT over different iterations.

Figure 9 :Figure 10 :

 910 Figure 9 : Conversion between bitmap image and 2D Boolean matrix.

Figure 11 :

 11 Figure 11 : Computation of the mismatch matrix.

Figure 12 :

 12 Figure 12 : Computation of the shape corrective function.

  Figure 14: Accuracy of each iteration during the optimisation process for all three target profiles

Figures 16 ,

 16 Figures 16, 17 and 18 represent the important geometries of the optimisation process respectively for the circular profile, the triangular profile and the rectangular profile. Figures a) are the geometries of the target profiles while figures b) are the geometries of the optimal tools resulting from the optimisation process. Finally, figures c) show the mismatch area remaining from the machining of the optimal tool.

Figure 16 :

 16 Figure 16: Tool optimisation results for the circular profile

Figure 17 :

 17 Figure 17: Tool optimisation results for the triangular profile

Figure 18 :

 18 Figure 18: Tool optimisation results for the rectangular profile

Figure 19 :

 19 Figure 19 : Evolution of accuracy during the optimisation process of the circular profile for different machining gaps.

Figure 20 :

 20 Figure 20 : Evolution of accuracy during the optimisation process of the triangular profile for different machining gaps.

Figure 21 :

 21 Figure 21 : Evolution of accuracy during the optimisation process of the rectangular profile for different machining gaps.


Ali Gharaviri

Mark Potse

Sander Verheule

Rolf Krause

Angelo Auricchio

Ulrich Schotten

Epicardial Fibrosis explains Increased Transmural Conduction in a Computer Model of Atrial Fibrillation

Recent work has shown that the transition from persistent to permanent AF in goats coincides with an increase in fibrosis in the outer millimeter of the atrial wall. Macroscopically this leads to reduced electrical conductivity orthogonal to the dominant fiber orientation.

We constructed a detailed geometry of the human atria including epicardial layer and all major endocardial bundle structures. The model also includes realistic one to three layers of fiber orientations, corresponding to their location in the atrium.

The numbers of waves, phase singularities, and breakthroughs (BTs) were quantified at different degrees of fibrotic tissue. Increase in the "fibrotic" volume from zero (Control) to moderate (50% Fibrotic), and severe (70% Fibrotic) increased both the number of waves and the number of phase singularities. Along with the increase in fibrosis, the endo-epicardial electrical activity dyssynchrony increased.

Introduction

The perpetuation and progression of atrial fibrillation (AF) is mediated by electrical and structural remodeling the atria [START_REF] Schotten | Pathophysiological mechanisms of atrial fibrillation: a translational appraisal[END_REF]. Electrical remodeling is a relatively fast process (1-2 days) that includes action potential shortening and reduction of effective refractory period. Structural remodeling develops much more slowly (months-years) and is considered as the main mechanism leading to AF stability. Several studies have demonstrated that structural remodeling is associated with a decrease in the size of fibrillatory waves and an increased incidence of transmural conduction ("breakthroughs") in both patients and goat models [START_REF] Verheule | Fibrillatory conduction in the atrial free walls of goats in persistent and permanent atrial fibrillation[END_REF][START_REF] De Groot | Electropathological substrate of longstanding persistent atrial fibrillation in patients with structural heart disease: epicardial breakthrough[END_REF][START_REF] Allessie | Electropathological substrate of long-standing persistent atrial fibrillation in patients with structural heart disease: longitudinal dissociation[END_REF]. Verheule et al. demonstrated that transition from persistent to permanent AF in goats is histologically characterized by an increase in endomysial fibrosis, particularly in the outer millimeter of the atrial wall [START_REF] Verheule | Fibrillatory conduction in the atrial free walls of goats in persistent and permanent atrial fibrillation[END_REF]. In this way, endo-epicardial transmural connectivity and electrical synchronization between endo-and epicardium would be hampered. Using simultaneous endo-epicardial mapping, Eckstein et al. showed a significant increase in endo-epicardial electrical dissociation (EED) as well as breakthroughs (transmurally propagating fibrillation waves) in goats in later stages of AF [START_REF] Eckstein | Time course and mechanisms of endoepicardial electrical dissociation during atrial fibrillation in the goat[END_REF][START_REF] Eckstein | Transmural conduction is the predominant mechanism of breakthrough during atrial fibrillation: evidence from simultaneous endo-epicardial high-density activation mapping[END_REF].

While these experimental studies show a clear association between structural remodeling, EED, and breakthrough incidence, they cannot establish a causal relationship or elucidate which aspects of the structural remodeling are responsible for the observed changes. The purpose of this study was to assess the effects of loss of side-to-side connectivity in the epicardial layer on EED and breakthrough incidence, using a numerical model in which all other factors were unchanged.

Method

3D model structure

An anatomical 3D model of human atria including fiber orientation was built manually using Blender (The Blender Foundation, Amsterdam, The Netherlands) as a 3D editing tool. The global shape was extracted from magnetic resonance imaging data of a subject with normal atrial anatomy. Endocardial and epicardial surfaces were created and combined into a 'myocardial envelope': a closed surface containing the atrial myocardium. 20 Pectinate muscles (PM) in the right atrial (RA) free wall, Bachmann's Bundle (BB), inter atrial bundles, and the crista terminalis (CT) were added based on published anatomical descriptions [START_REF] Ho | Atrial structure and fibres: morphologic bases of atrial conduction[END_REF][START_REF] Ho | The importance of atrial structure and fibers[END_REF][START_REF] Anderson | Cardiac Anatomy: An Integrated Text and Colour Atlas[END_REF]. The model contained four pulmonary veins, caval veins, left and right atrial appendages (LAA and RAA), as well as the coronary sinus, as illustrated in figure1A, B and C.

Dedicated software was used to create a computational mesh at 0.2-mm resolution with different fiber orientations from the manually-edited objects. Within the myocardial envelope a 1-mm layer of myocardium was placed to represent the thin right-atrial wall. A 3-mm layer was used for the thicker left atrium. The bundles consisted of myocardium with different conduction properties.

Fiber orientation and anisotropy

One of the most important properties of this model is the incorporation of realistic fiber orientation based on histological observations (see figure 1C andD) [START_REF] Ho | Atrial structure and fibres: morphologic bases of atrial conduction[END_REF]. Ho et al have described the left-atrial wall as 2 layers with different fiber orientations [START_REF] Ho | Atrial structure and fibres: morphologic bases of atrial conduction[END_REF]. Less information is available for the right atrium.

In our model, elements within the bundles received a fiber orientation parallel to the local bundle orientation.

Electrophysiological model

Simulations were performed with the propag-5 software and run on 2560 cores of a Cray XE6 supercomputer [START_REF] Potse | A comparison of monodomain and bidomain reactiondiffusion models for action potential propagation in the human heart[END_REF][START_REF] Krause | Hybrid Parallelization of a Large-Scale Heart Model[END_REF]. Electrical activity was simulated with a mono-domain reaction-diffusion equation. Ionic currents and calcium handling for each model node was described by Courtemanche-Ramirez-Nattel model [START_REF] Courtemanche | Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model[END_REF].

To reproduce electrical remodeling as observed in atrial fibrillation, conductivities for Ito, IcaL, and Ik1 were set at 15%, 15%, and 200% of their normal values respectively [START_REF] Schotten | Pathophysiological mechanisms of atrial fibrillation: a translational appraisal[END_REF].

Simulation protocol

A spiral was initiated via a single pacing site in combination with an associated temporary block line.

Implementation of fibrosis

Fibrosis segments were simulated by randomly removing transverse conductivity between elements. Three different degrees of fibrosis (slight (50%), moderate (60%), and severe (70%)) were selected and simulated to be compared with no fibrosis (control).

Analysis

Detection of phase singularities

In two dimensions, the tip of a reentrant wave forms a phase singularity (PS); a point where all phases of the action potential simultaneously occur. The threedimensional equivalent of a PS is a vortex filament. Model elements that were part of vortex filaments were detected using the method [START_REF] Fenton | Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation[END_REF]. A filament was defined by the line of intersection of the two surfaces 𝑉𝑉(𝑟𝑟, 𝑡𝑡) = 𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉, where

𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉 is an arbitrary selected potential value which used to distinguish polarized and depolarized regions of the tissue (in this study -49 mV), and 𝜕𝜕 𝑡𝑡 𝑉𝑉(𝑟𝑟, 𝑡𝑡) = 0. Since in our simulations the time is discrete we could approximate the 𝜕𝜕 𝑡𝑡 𝑉𝑉(𝑟𝑟, 𝑡𝑡) = 0 by 𝜕𝜕 𝑡𝑡 𝑉𝑉(𝑟𝑟, 𝑡𝑡) ≈ 𝑉𝑉(𝑟𝑟, (𝑛𝑛 + 1)∆𝑡𝑡) -𝑉𝑉(𝑟𝑟, 𝑛𝑛∆𝑡𝑡). Therefore, the filament position could be identified using formula [START_REF] Schotten | Pathophysiological mechanisms of atrial fibrillation: a translational appraisal[END_REF].

𝑉𝑉(𝑟𝑟, 𝑛𝑛∆𝑡𝑡) = 𝑉𝑉(𝑟𝑟, (𝑛𝑛 + 1)∆𝑡𝑡) = 𝑉𝑉 𝑖𝑖𝑖𝑖𝑖𝑖

(1)

In order to describe this method more practically we represent the filament line by the set of intersection points of this line with the lattice planes defined 𝑥𝑥 = ∆𝑉𝑉 , y= ∆𝑗𝑗 , 𝑧𝑧 = 𝑘𝑘∆𝑧𝑧, which they are all the y-z, x-z, and x-y planes that contain lattice points, respectively. Since each plane is made up of squares with one lattice point at each of the four corners, the algorithm proceeds by finding which of the 3NxNyNz squares contained in these planes are intersected by the line defined by Eq. ( 1) and what are the coordinates of these intersection [START_REF] Fenton | Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation[END_REF]. Finally, all of these points are sorted in order of increasing arclength S along a given filament, the arclength distance between two nearby points (𝑝𝑝(𝑝𝑝 1 , 𝑝𝑝 2 , 𝑝𝑝 3 ) and 𝑞𝑞(𝑞𝑞 1 , 𝑞𝑞 2 , 𝑞𝑞 3 )) being simply defined by the Euclidean distance (𝑑𝑑(𝑝𝑝, 𝑞𝑞) = �(𝑝𝑝 1 -𝑞𝑞 1 ) 2 + (𝑝𝑝 2 -𝑞𝑞 2 ) 2 + (𝑝𝑝 3 -𝑞𝑞 3 ) 2 ).

Detection of fibrillation waves

A wave was defined as a contiguous area in which all nodes have trans-membrane voltages above the excitation threshold of -60mV. The number of waves was calculated each millisecond of simulated time.

Dyssynchrony

Dyssynchrony was defined as the percentage of nodes that were excited (membrane potential above -60 mV), while the opposing segment of the atrial wall was not excited (membrane potential below -60mV).

In order to calculate dyssynchrony, model nodes were categorized either as epicardial or endocardial. Pairs of endocardial and epicardial nodes were created by finding the nearest epicardial node for each endocardial node.

Dyssynchrony was defined as the number of pairs with different state (one node activated and the other not activated) divided by the total number of pairs. Dyssynchrony was computed for each millisecond of simulated time.

Detection of breakthroughs

A breakthrough (BT) is a wave that appears in epicardial layer and cannot be related to the propagation of other waves in that layer. To detect BTs, wave sizes at each time step were calculated. If a wave smaller than a threshold (9 nodes in this study) appears within the epicardial layer and this wave has an overlap with a wave in endocardial layer, this was considered as a BT candidate. Each BT candidate was monitored and if its size increased in two folds, it was confirmed as a BT.

Results

Effect of fibrosis

Fibrillation patterns

In figure 2, an example of fibrillation patterns in both control and severe fibrosis are depicted. In the presence of fibrosis, fibrillation patterns in RA, LA anterior wall, and pulmonary veins area were more complex compared to control. Longitudinal dissociation, electric uncoupling of side-to-side connections between parallel-oriented atrial muscle bundles, within the epicardial plane increased and waves became narrower.

Dissociation and BT

Quantitative analysis of fibrillation patterns

The numbers of waves, phase singularities, and breakthroughs (BTs) were quantified at different degrees of fibrotic tissue. Increase in the "fibrotic" volume from zero (Control) to slight (50%), moderate (60% Fibrotic), and severe (70% Fibrotic) increased both the number of waves and the number of phase singularities. Along with the increase in fibrosis, the endo-epicardial electrical activity dyssynchrony and the incidence of BTs increased (see figure 4).

5.

Discussion

In this study we present a highly detailed 3D anatomical model of the human atria. This model was used to test the hypothesis that fibrosis in the epicardial layer would be sufficient to produce increased endo-epicardial dissociation of electrical activity and as a consequence an increase in transmural conduction. Our simulations demonstrate that epicardial fibrosis increased the number of waves and PSs, endo-epicardial dissociation of electrical activity and the BT rate.

Study limitations

There were some limitations regarding this study. (1) Human experimental and histological data were sparse and limited e.g. the histological data regarding fiber orientation in the RA of the human heart is very limited. For this reason, the fiber orientation in the RA thin wall was implemented using derived histological data from goats.

(2) Several studies have proposed AP heterogeneity as a main mechanism for AF initiation and perpetuation, but in this model only one myocyte model with the same characteristics was used in different region of atrium. Our model offers an excellent opportunity to investigate the presence of this property in the future. (3) Random focal discharges, which have been investigated in several AF computer models and can be one of the important AF mechanisms, are also not implemented in this model [START_REF] Gong | Mechanism underlying initiation of paroxysmal atrial flutter/atrial fibrillation by ectopic foci: a simulation study[END_REF].

Figure 1 .

 1 Figure 1. Anatomical model of the atria used for simulations. A) Posterior view. B) Saggital view. C) Fiber orientations in anterior view. D) RA free wall fiber orientations.

Figure 2 .

 2 Figure 2. A) Potential map in LA and RA (control). B) Potential map in LA and RA (severe fibrosis).

Figure 3

 3 Figure 3 illustrated examples of conduction patterns in control (3A and C) and severe fibrosis (3B and D).As illustrated in this figure in both cases (control and severe fibrosis) BTs occurred but incidence of BTs in control group was more repetitive and limited to only few locations. An increase in the degree of fibrosis enhanced endo-epicardial electrical dissociation which increased the chance for BT occurrence and therefore more BTs were observed at different locations. For example, in none of the control group simulations a BTs was observed in RA free wall near the PMs or in anterior side of LA but multiple breakthroughs occurred at those sites in the presence of fibrosis (figure3B and D).

Figure 3 .

 3 Figure 3. A) Incidence of BT in (black arrow) in control simulation. B) Incidence of BT (black arrow) in severe fibrosis simulation. C) Electrical activity association between PMs and right atrial epicardial layer. D) Increase in electrical activity dissociation in PMs and right atrial epicardial layer and incidence of BT.

Figure 4 .

 4 Figure 4. A) Number of waves, B) phase singularities, C) percentages of endo-epicardial electrical activity dissociation, and D) summation of all occurred BTs during simulation period in different degrees of fibrosis.

Computing in Cardiology 2016; 43:237-240.

Address for correspondence.

Ali Gharaviri Institute of Computational Science, Faculty of Informatics Via Giuseppe Buffi 13,6900 Lugano, Switzerland ali.gharaviri@usi.ch


R Chitic

R Bearee

F Giraud

B Lemaire-Semail

M Favre

P Tierce

J Jehanno

Industrial equipment for Powder transportation using piezoelectric "friction control" method

Keywords: powder transportation, ultrasonic vibration, friction control

This paper presents a new powder transportation system that uses a high frequency flexural stationary wave coupled with a low frequency horizontal displacement of a beam to produce the transport of the powder. The ultrasonic wave is produced with the help of piezoelectric cells glued under the beam and is used to decrease the friction coefficient between the powder and the beam surface.

Introduction

Powder transportation is sometimes a delicate stage of an industrial process, although it gives the final product no added value: in general, it is expected that the transported powder does not segregate nor suffers chemical alteration and that no loss of material is observed. Today many powder transport systems exist but they suffer from a significant number of defects to be corrected: conveyor belts or vibratory conveyor are often noisy mechanical systems, requiring a massive energy consumption and costly maintenance due to their moving parts. To cope with those problems, Ultrasonic Powder Transport Systems use ultrasonic vibrations to transport powder. For example, [START_REF] Mracek | A system for powder transport based on piezoelectrically excited ultrasonic progressive waves[END_REF] [START_REF] Kühne | Dosage of highly disperse powders by ultrasound agitated tube modules[END_REF] propagates an ultrasonic progressive wave into a specific pipe material with appropriate absorption properties. They obtained a simple device with no moving part which can transport powder with good performances. However, the pipe material needs to absorb the reflective wave, dissipating the electrical power into heat. In this paper, we use a high frequency flexural stationary wave coupled with a low frequency horizontal displacement. Our goal is to decrease the amplitude of horizontal displacement compared to classical technologies, while ultrasonic vibrations are imposed with a good efficiency because they are stationary. In the first section, we describe the friction reduction principle which is used to propel the powder, then a description of the test set-up is given, while we present in the final section the experimental results.

Friction control principle

A beam is submitted to a sinusoidal low frequency movement (LFM) in the direction of the excited transport.

Fig. 1. Illustration of the friction control principle

In order to create powder transportation, a dissymmetry in contact condition between the powder and the beam has to be introduced synchronized on the LFM. Such a dissymmetry is obtained by decreasing the friction coefficient. This is achieved if a high frequency vibration (HFV) at 20kHz is produced [3] [4]. Its amplitude is modulated by a signal named (MS), synchronized with the LFM. The synchronization principle consists in using HFV only when the friction coefficient has to be decreased. For transportation creation, the reaction force transmitted to the powder by the transportation system (for instance a beam) has to be only directed according to the desired direction of transportation. Figure 1 presents the friction control principle. Hence, when the LFM acceleration is opposite to the desired direction, the HFV is introduced to reduce the contact coefficient with the transportation device. When the HFV is used, one can assumed that the powder motion is still in the expected direction due to inertial effect.

Experimental set-up

Figure 2 depicts the tests bench, which will be used for the experimental investigations. A shaped beam is used to eliminate the possibility that the powder will fall off the plate [START_REF] Hashimoto | Transporting objects without contact using flexural traveling waves[END_REF]. A linear motor is used to create the LFM of the beam. In addition, we glued piezoelectric cells under the beam to create the HFV. One can note that the location of the piezoelectric cells is optimized using modal analysis to correspond to the maximum amplitude locations of the flexural vibrational mode excited by HFV.

Fig. 2. Tests set-up overview

The elastic supports of the beam needs to be positioned in precise locations determined using the modal analysis of the beam. They need to be mounted on the nodes position of the HF flexural wave, to avoid the perturbation of flexural deformation [START_REF] Szlabowicz | Contribution au dimensionnement et à la réalisation d'actionneur piézoélectrique à rotation de mode fort couple pour applications aéronautiques[END_REF]. The voltage supplied to the ceramics is synchronized with the beam's position, which is measured by a position sensor. Closed loop control is achieved by a computer embedding DSpace controller. The test bench was designed to be easily inclined for helping the powder transportation.

Results

A series of preliminary tests have been conducted. They consist in using only the high frequency flexural stationary wave in order to establish the influence of the ultrasonic vibration on the friction reduction between the powder and the beam. During the test we changed the quality surface of the beam by adding a layer of tape with different thickness h: tape #1 -h=0.05mm and tape #2 -h=0.12mm. Flow rate (g/s)

Fig. 3. Flow rate in function of HF amplitude for a corundum powder with a 15° inclination of the system and Tape #1 quality of the beam surface.

In the figure 3, we show the evolution of the flow rate as a function of the HFV amplitude. Such a curve is generic for all tested powder. It shows that if there is no vibration, no powder transportation is obtained, which is normal because there is no dissymmetry in contact conditions then.

On the same figure, we can see that at least 400nm p-p is required to produce powder transportation, and that 1000nm p-p produces good transport conditions. Figure 4 gives the evolution of the flow rate as a function of HFV amplitude for different qualities of the beam surface. These trials have been done with a 15° inclination of the system. As we can see, higher tape thickness produces lower flow rate.

-

Conclusion

This paper presents a new powder transportation system. It consists in a beam which is excited in a low frequency movement. A dissymmetry of the contact conditions according to the low frequency displacement is produced by creating flexural vibration of the beam. We found that ultrasonic vibrations with amplitude synchronized with the beam position could help to propagate the powder.

  frequency amplitude p-p (nm)

Fig. 6 .

 6 Fig. 6. Flow rates in function of low frequency displacement for different powders

  Now we present the results obtained using HFV synchronized with the LF displacement. In figure5we can see the curves for the flow rate in function of the LF displacement with different inclinations α of the system. The HFV amplitude is kept constant for the two inclinations but not equal: we have 830nm peak amplitude for α=15° and 530nm peak amplitude for α=0°. The flow rate is more important for a 15° inclination of the system, because the transversal component of the weight force is added to the force displacement produced by the LF movement of the beam.In the figure6, we plotted the curves of the flow rates for different types of powders in function of the horizontal displacement and for constant HFV amplitude. These tests have been done considering a 15° constant inclination of the system. As we can see the arrangement of the curves depend on the morphological characteristics of the powders (in ex. average diameter of the particles noted Φ).

		0.8					1.2
				Tape #1				α=15°α
	Flow rate (g/s)	0.2 0.4 0.6	Stainless steel Tape #2				Flow rate (g/s)	0.3 0.6 0.9	=0°
			0
			0	500	1000	1500	2000	0
		0.2					0	200	400	600	800	1000
				HF amplitude p-p (nm)			LF displacement p-p (µm)
	Fig. 4. Flow rate in function of HF amplitude for	Fig. 5. Flow rate in function of LF displacement for
	different beam surface using a corundum powder.	different inclination of the system using a corundum
								powder. Stainless steel quality of the beam surface.
				Flow rate in fonction of LF displacement, stainless steel surface,
			1.2				inclination α=15°0
				Corundum Φ50µm HFA=830nm
			1	Alumina Φ30µm HFA=820nm Cement Φ17µm HFA=810nm
	Flow rate (g/s)	0.4 0.6 0.8	Alumina Φ10µm HFA=850nm
			0.2
			0		200		400	600	800	1000
							LF displacement p-p (µm)


Francois Grinnaert

Jean-Yves Billard

Jean-Marc Laurens

KG max Curves Associated With Second Generation Intact Stability Criteria for Different Types of Ships

Keywords: pure loss of stability, parametric roll, KGmax curve, second generation intact stability criteria

Currently, second generation intact stability criteria are being developed and evaluated by the International Maritime Organization (IMO). In this paper, we briefly present levels 1 and 2 assessment methods for the criteria of pure loss of stability and parametric roll failure modes. Subsequently, we show the KGmax curves associated with these criteria. We compute these curves for five different types of ships and compare them with the curves embodied in the current regulations. The results show that the safety margin ensured by the first level-1 method of calculation for both pure loss of stability and parametric roll seems to be excessive in many cases. They also show that the KGmax given by the second level-1 method and by the level-2 method may be very similar. In some cases, the level-2 method can be more conservative than the second level-1 method, which is unanticipated by the future regulation. The KGmax curves associated with parametric roll confirm that the C11 container ship is vulnerable to this failure mode. The computation of the second check coefficient of parametric roll level 2 (C2) for all possible values of KG reveals the existence of both authorized and restricted areas on the surface formed by both the draft and KG, which may replace the classical KGmax curves. In consequence, it is not sufficient to check that C2 is lower than the maximum authorized value (RPR0) for a fixed ship's loading condition.

1 Introduction 1 Secondgeneration intact stability criteria are currently being developed and evaluated by the International Maritime Organization (IMO) [START_REF] Bassler | A review of available methods for application to second level vulnerability criteria[END_REF][START_REF] Francescutto | Current status of new generation intact stability criteria development[END_REF][START_REF] Wandji | Test application of second generation IMO intact stability criteria on a large sample of ships[END_REF]. These criteria were presented in detail at the 13 th International Ship Stability Workshop [START_REF] Umeda | Current status of second generation intact stability criteria development and some recent efforts[END_REF]. In this paper, we address this version, as amended in February 2015 and January 2016 by the Sub-Committee on Ship Design and Construction of the IMO (SDC 2/WP.4 and SDC 3/WP.5). The new criteria are structured in five failure modes: pure loss of stability in astern waves, parametric roll, broaching/surf-riding, dead ship condition, and excessive * acceleration. Three assessment levels are defined for each failure mode. The first level is based on a simplified deterministic approach of the phenomena and ensures high safety margins. A ship that does not pass this first level is designated as "unconventional" by the new regulation (Fig. 1). The second assessment level is based on probabilistic approaches to the phenomena. It requires more complex computations based on hydrostatic considerations with regard to static waves, and ensures medium safety margins. The third assessment level consists of numerical simulations of the ship's behavior on waves, as performed by specialized institutes. This last level ensures optimized safety margins. In order to comply with future rules, a ship must comply with current rules (IMO IS Code 2.2, 2009) and with at least one level of each failure mode (Fig. 1).

After a short presentation of the first and second level criteria of the pure loss of stability and parametric roll failure modes, in this paper, we show the curves of the maximum height of the center of gravity as a function of displacement (KGmax curves) with regard to the criteria of these two failure modes. We compute these curves for ships of five different types, chosen for their variety of behavior with regard to both failure modes: a naval ship, two container ships, a roll-on roll-off (Ro-Ro) vessel, and a tanker. KGmax curves allow a crew to assess their ship's stability with regard to a specific criterion or several criteria (envelope curve) without the need to performa any complex computation. KGmax curves are included in each ship's stability booklet. In this paper, the avoidance any assumption about the height of the center of gravity is permitted. This also enables a comparison of the results given by all levels and methods in order to determine the efficiency and relevance of the new criteria. We compare curves based on these new criteria with those based on the first-generation criteria. While not the main objective of this paper, the consideration of KGmax curves allows for the determination of the vulnerability of vessels to these two failure modes. Similar studies have been performed [START_REF] Krueger | Investigation of the 2nd generation of intact stability criteria for ships vulnerable to parametric rolling in following seas[END_REF][START_REF] Krueger | On the application of the 2 nd generation intact stability criteria to Ro-Pax and container vessels[END_REF] but were based on former versions of the new criteria (respectively SLF 54/19, 2012 andSDC 1/INF.8, 2013).

Presentation of ships

The main particulars of the five ships in this study are given in Table 1. The superstructures of four of them have been modelled in order to compute the weather criterion (IS Code 2.3, IMO, 2009), without taking into account their flooding points. The bilge keels of all the ships have also been modelled. The first ship is the David Taylor Model Basin hull number 5415, a naval ship presented by [START_REF] Moelgaard | PMM-tests with a model of a frigate class DDG-51[END_REF]. Included in her model are superstructure concepts inspired by those of the DDG51 Arleigh Burke.

The second ship is the C11 container ship, which is known to be vulnerable to parametric roll [START_REF] France | An investigation of head-sea parametric rolling and its influence on container lashing systems[END_REF].

The third ship is a 319 m container ship. A fatal accident due to an extreme roll motion (not a parametric roll) occurred on a ship of this class in 2008. One of the causes determined by the experts was insufficient roll damping [START_REF] Kaufmann | Fatal Accident on Board the CMV Chicago Express During Typhoon "Hugupit" on September 24 2008 off the Coast of Hong Kong[END_REF]. Numerical models of both container ships include the containers in full-load configuration.

The fourth ship is a Ro-Ro vessel presented by [START_REF] Garme | Model seakeeping experiments presented in the time-domain to facilitate validation of computational tool[END_REF]. Data for this ship have been provided by the KTH Royal Institute of Technology. Her superstructure has not been modelled. The KGmax curve associated with the first-generation criteria does not take weather criterion into account.

We chose the last ship from the DELFT ship database-a tanker that is wall-sided on 80% of her length-for her supposed non-vulnerability to both pure loss of stability and parametric roll. A simplified box-shaped superstructure was added to her model but has no influence on the KGmax curve associated with the current criteria.

We used the reference height of the center of gravity (KGref) to calculate the associated metacentric height in calm water (GMref) and natural roll period (T0) for each ship. The latter depends on the added mass coefficient and the radius of inertia coefficient (a and k, as defined in subsection 4.2. 1). The roll period is given for the full-load displacement, except for the C11 container ship, whose associated draft is 12.339m [START_REF] France | An investigation of head-sea parametric rolling and its influence on container lashing systems[END_REF]. We calculated the roll damping coefficient B44 as described in subsection 4.2.2. Its value is given for the reference height of the center of gravity (KGref), the full load draft (12.339 m for the C11 container ship), the service speed (Vs), and a roll amplitude of 10 degrees.

3 Pure loss of stability failure mode

Presentation of criteria

Here, we briefly present the first-and second-level criteria for the pure loss of stability failure mode. These criteria are applied to ships having a service speed Froude number larger than 0.24. Although neither container ship nor the tanker fulfill this condition, we computed their KGmax curves that are associated with the levels 1 and 2 criteria.

First level

To meet this criteria level, the minimum metacentric height in waves (GMmin) must be higher than 5 centimeters. The new regulation proposes two methods to calculate its value.

The first method considers a parallel water-plane area at the lowest draft dL:

full 0.25 ; 2 W L LS d d d d          (1)
with: d (m) draft corresponding to the loading condition in calm water (black line in Fig. 2); dfull (m) draft, full load; L (m) ship's length; SW 0.0334 wave steepness. In most cases, the lowest draft can be simplified as follows:

2 W L LS d d   (2)
The minimum metacentric height is calculated as follows:

min GM KB KG L I     (3)
with: KB (m) height of the vertical center of buoyancy in calm water; IL (m 4 ) moment of inertia of the water-plane area at lowest draft dL (gray area in Fig. 2);  (m 3 ) volume of displacement; KG (m) height of the vertical center of gravity. This method considers that the moment of inertia of the water-plane area in waves is similar to that of the parallel water-plane area at the lowest draft dL (gray area in Fig. 2).

In the second method, the minimum metacentric height is computed for the ship balanced in trim and sinkage on sinusoidal waves. The wave length is equal to the ship's length L and the height is h=SWL. The wave crest is centered at the longitudinal center of gravity and at each L/10 forward and aft.

Second Level

Determination of the second-level criterion consists of a probabilistic approach to the phenomenon, associated with a wave scattering table. For an unrestricted sailing area, the new regulation imposes that found in the IACS Recommendation No. 34 (2001), which corresponds with the Northern Atlantic A ship is found to be not vulnerable to pure loss of stability if both CR1 and CR2 values are lower than RPL0 = 0.06. CR1 and CR2 are computed as follows:

197 197 1 1 CR1 C1 CR2 C2 i i i i i i W W       (4)
with:

Wi weighting factor of the wave number i taken from Table 2 (number of occurrences divided by 100,000); C1i coefficient for the wave number i regarding the minimum angle of vanishing stability V.min (Fig. 3); C2i coefficient for the wave number i regarding the maximum angle of stable equilibrium S.max under action of a heeling lever RPL3 (Fig. 3). Waves are assumed to be sinusoidal. Their heights are taken from Table 2 (HS is twice the wave amplitude), and their lengths  are linked to their zero-crossing periods TZ (Table 2) by the infinite depth relation given by:

2 2π Z gT   (5)
For each wave, the heeling lever RPL3 is defined as follows:

2 3 8 s PL n H R d F   (6)
Then the minimum angle of vanishing stability (V.min) and the maximum angle of stable equilibrium (S.max) are computed for the 197 non-zero-weighted waves of the wave scatter diagram, and used these values to calculate C1i and C2i as follows:

min 1 1 if < C1 0 otherwise V . PL i R      (7) min 2 1 if < C2 0 otherwise S . PL i R      (8)
with:

RPL1 30 degrees; RPL2 25 degrees (15 degrees for passenger ships).

Since there are a large number of non-zero-weighted waves, the new regulation proposes a method to make the computation faster, based on Grim's effective wave height concept (1961). For each wave of Table 2 defined by its zero-crossing period (TZ) and significant height (HS), an effective height for the top 3% highest waves H3% eff is computed. It corresponds to an equivalent wave whose length  is equal to the ship's length L. The effective wave height is computed using the following formulae:

eff 3% 0 5.9725 H m  (9)   2 2 2 3 5 4 0 2 2 0.01 2 sin 2 exp d π 2 L L L L g g m A B L g                                        (10) 2 4 01 173 S A H T   (11) 4 01 691 B T   (12) 01 1.086 Z T T  (13) 2 π L g L   (14)
The minimum angle of vanishing stability (V.min) and the maximum angle of stable equilibrium (S.max) are calculated for 11 effective waves with heights from zero to the maximum effective wave height corresponding to the ship's length and the IACS wave scatter diagram (Table 2). For each effective wave, 10 GZ curves are computed, considering the ship balanced in trim and sinkage, and for the wave crest centered at the longitudinal center of gravity and at each /10 forward and aft. Only one GZ curve is needed for the first wave, which corresponds to calm water.

The values of V.min and S.max associated with any effective wave height are calculated using linear interpolation in the interval between two successive values of these 11 heights. Table 3 gives the effective wave heights (in meters) computed according to the Grim method for all the ships in this study, as associated with the waves of Table 2 having a zero-crossing period (TZ) equal to 10.5 seconds.

Results

We computed the KGmax curves associated with the pure loss of stability levels 1 and 2 criteria for the five ships in this study. Hydrostatic computations in calm water and static waves were made using Calcoque software, which is published in a handbook [START_REF] Grinnaert | Stabilité du navire -Théorie, réglementation, méthodes de calcul -Cours et exercices corrigés. Ellipses, collection Technosup[END_REF] and was presented at the 12 th International Conference on the Stability of Ships and Ocean Vehicles [START_REF] Grinnaert | Calcoque: a fully 3D ship hydrostatic solver[END_REF]. All computations were performed assuming a water density of 1.025 t/m 3 and zero trim (the sinkage and trim remain free while the ships are balanced on waves). The KGmax curves are shown in Fig. 4 to Fig. 8. The KGmax curve determined by the first level-1 method (parallel water-plane) is marked "PL Level 11st method" in the figures (square markers). The KGmax curve determined by the second level-1 method (ship balanced in trim and sinkage on a wave of the same length) is marked "PL Level 12nd method" (round markers). The KGmax curve determined by the level-2 criterion is marked "PL Level 2" (diamond markers). A ship may be considered to be vulnerable to pure loss of stability if her KGmax curve associated with the current regulation (civilian or military) is located above the highest curve associated with the pure loss of stability criteria (we formulated this consideration in this paper, for both pure loss of stability and parametric roll failure modes). From these figures, we can observe the following:

1) The two level-1 methods yield significantly different results for all vessels, except the tanker.

2) The KGmax value given by the second level-1 method is relatively close to that given by the level-2 method for the C11 container ship and the Ro-Ro vessel.

3) The level-2 method can yield more conservative results than the second level-1 method (by one point for the C11 container ship and more markedly for the naval ship and the Ro-Ro vessel). This is contrary to our understanding of the intent regarding future regulation.

4) The level-2 method allows for a negative GM for the tanker. This is due to the consideration of the GZ curve instead of GM in this criterion. Fig. 3 shows a GZ curve with a negative GM, which has a positive contribution to the criterion (C1i and C2i are not increased): V is larger than RPL1 (30 degrees) and S is lower than RPL2 (25 degrees).

5) According to these criteria, the naval ship and the tanker are not vulnerable to pure loss of stability, but both container ships are vulnerable, irrespective of the Froude number. No conclusions can be drawn with regard to the Ro-Ro vessel, since her superstructure has not been modelled and the weather criterion is not taken into account.

Parametric roll failure mode 4.1 Presentation of criteria

Next, we briefly present the first-and second-level criteria for the parametric roll failure mode. For more details please refer to SDC 2/WP.4 and SDC 3/WP.5.

First level

The first-level criteria of the parametric roll requires that the amplitude of the variation of the metacentric height in waves (GM) be less than RPR, whose value is between 0.17 (ship without a bilge keel) and 1.87 (ships with a large bilge keel). The new regulation proposes two methods for calculating GM.

In the first method, two parallel water-planes are considered. The first is a low draft (dL) water-plane, and its moment of inertia is IL. The second is a high draft (dH) water-plane, and its moment of inertia is IH, with:

full min 0.25 ; 2 W L LS d d d d          (15) min ; 2 W H LS d d D d          (
W L LS d d   (17) 2 W H LS d d   (18)
The amplitude of the variation of the metacentric height is calculated as follows: GM 2

H L I I     (19)
In the second method, 10 metacentric heights are computed for the ship balanced in trim and sinkage in sinusoidal waves. The wave length is equal to the ship's length L and the height is h = SWL. The wave crest is centered at the longitudinal center of gravity and at each L/10 forward and aft. We use the minimum and maximum values of GM to calculate GM, as follows:

max min GM GM GM 2    (20)

Second level

For the second level criterion, a probabilistic approach to the phenomenon is considered, by which a ship is found to be not vulnerable to parametric roll if the values of C1 and C2 are lower than RPR0 = 0.06.

The value of C1 is computed as follows, based on a table that lists 16 weighted waves (Table 4, calculated from the wave scatter diagram using the method described in SDC 3/WP.5):

16 1 C1 C i i i W     ( 21
)
where Wi is the wave weighting factor from Table 4, and Ci is equal to 0 if condition A or B is satisfied, and is equal to 1 otherwise.

i i i R    (22)
with: GMi (m) average value of 10 metacentric heights computed for the wave number i, with the wave crest centered at the longitudinal center of gravity and at each i/10 forward and aft;

GMi (m) one-half of the difference between the maximum and minimum values of the 10 metacentric heights. Condition B is as follows:

. PR i S V V  (23)
with: VS (m/s) ship's speed; VPR.i (m/s) reference speed corresponding to the parametric resonance conditions, calculated as follows:

. 0 2 2π

i i i PR i c GM V g T GM     (24)
with:

i (m) wave length (from Table 4); T0 (s) natural roll period of the ship in calm water; GMi (m) as defined above; GMc (m) metacentric height in calm water. The C2 coefficient is calculated from the maximum roll angle of the ship on each of the 197 non-zero-weighted waves of the wave scatter diagram (Table 2, or another table if the ship is sailing in a restricted area), considering seven speeds. For a given wave and a given speed, the maximum roll angle is the maximum absolute value of the solution (t) of the differential equation of the parametric roll, which can be established as follows: The C2 coefficient is calculated as follows:

  44 44 GZ , 0 J B W t          (
197 7 , 1 1 1 C2 .C2 7 i i j i j W     (26)
The weighting factor Wi is extracted from the wave scattering table (Table 2). The coefficient C2i,j is equal to 1

if the maximum roll angle of the (t) function solution of the differential equation for wave number i and the speed KjVS is higher than 25 degrees, and is 0 otherwise. The speed factors Kj are given in Table 5. Since there are a large number of non-zero-weighted waves in Table 2, the new regulation proposes the use of the [START_REF] Grim | Beitrag zu dem problem der Sicherheit des Schiffes im seegang[END_REF] effective wave height concept to make the computation faster. This methodology is similar to that used for pure loss of stability. However, the effective height is considered to be one of the top 1/3 highest waves, and is calculated as follows:

eff 1/3 0 4.0043 H m  (27)
The maximum roll angle max is computed by solving the differential equation Eq. ( 25) for 10 effective waves whose heights are a fraction (from 0.1 to 1) of the maximum effective wave height, corresponding to the ship's length and the wave scatter diagram (Table 2). Linear interpolation is used to calculate the maximum roll angle associated with any effective wave height. Table 6 gives the effective wave heights (in meters) computed using the Grim method for all ships and associated with the waves of Table 2 that have a zero-crossing period (TZ) equal to 10.5 seconds.

The future regulation proposes to calculate the maximum roll angle by analytically solving the differential equation. It also allows for a one-degree-of-freedom time domain simulation to solve the differential equation. Some solving conditions have recently been specified (initial roll angle, initial roll speed, time step, simulation duration, see SDC 3/WP.5). In this study, we used this second method, and specify below the methods used and the conditions for solving the differential equation and computing the value of KGmax.

Maximum roll angle and KGmax computation methods

We have already presented the differential equation to be solved as Eq. ( 25) above.

Moment of inertia

We calculate the moment of inertia J44 as follows: metacentric height in calm water. When the ship's roll period is known, we then choose the coefficients a and k. If it is not known, we set coefficient a to 0.1 and k to 0.4. The values of both a and k are given in Table 1 for all ships.

   2 44 Δ 1 J a kB   (

Roll damping

We calculate the roll damping coefficient B44 using the method proposed by [START_REF] Kawahara | A simple prediction formula of roll damping of conventional cargo ships on the basis of Ikeda's method and its limitation[END_REF], and take into account the ship's speed using the method proposed by [START_REF] Ikeda | Components of roll damping of ship at forward speed[END_REF]. The roll damping coefficient B44 depends on the roll amplitude (a in [START_REF] Kawahara | A simple prediction formula of roll damping of conventional cargo ships on the basis of Ikeda's method and its limitation[END_REF]. In this study, we pre-computed B44 for 11 roll amplitude values from 1 to 50 degrees and used linear interpolation using the value of the current roll amplitude during the solving of the differential equation. The roll damping coefficient B44 also depends on the distance between the center of gravity and the calm water surface (OG in [START_REF] Kawahara | A simple prediction formula of roll damping of conventional cargo ships on the basis of Ikeda's method and its limitation[END_REF] and on the ship's roll period (T0). Both are linked to the KG. Thus, pre-computation of B44 is required for all roll amplitudes in all iterations of the search for KGmax.

Righting arm

We calculate the righting arm using a simplified method proposed by [START_REF] Peters | Requirements for computational methods to be used for the IMO second generation intact stability criteria[END_REF] and the new regulation (SDC 3/WP.5, IMO, 2016), as follows:

        0 0 m o y GZ GZ GM GM ΔGM sin sin e t          (30)
with: GZ0() (m) righting arm in calm water associated with a user-fixed height of the center of gravity (KG0); GM0 (m) metacentric height in calm water associated with KG0. GMmoy (m) average value of the metacentric height in waves, associated with the considered KG, which may differ from KG0; GM (m) amplitude of the variation of the metacentric height in waves; e (rad/s) wave encounter frequency. [START_REF] Belenky | Development of second generation intact stability criteria. US Navy[END_REF] demonstrated that this simplified method yields acceptable GZ values up to an angle of 30 degrees. Since the C2 coefficient is increased if the maximum roll angle exceeds 25 degrees, we can used this method in this study and it speeds up the computation of the maximum roll angle and KGmax. We pre-computed the righting arm in calm water (GZ0) for heel angles from 0 to 50 degrees with a step of 5 degrees. During the simulation, we used linear interpolation to calculate its value for any roll angle . We computed GMmoy and GM before the simulation by computing the GM of the considered sinusoidal wave for 10 ship positions on the wave. GMmoy is the average value of the 10 metacentric heights. GM is half the difference between the maximum and minimum values. During the search for KGmax, we updated the value of GMmoy with the considered value of KG.

We used GMmoy to calculate the ship's natural roll period T0, which is needed to calculate the damping coefficient B44. Thus, if GMmoy is negative, simulation is not possible and the C2 coefficient value is forced to 1.

We calculated the wave encounter frequency e as follows:

1 j s e w w K V g            (31) with:
w (rad/s) wave frequency.

The positive values of the speed factors Kj (from Table 5) correspond to head seas, and negative values correspond to following seas.

Solving method

To solve the differential equation, we used the Runge-Kutta 4 th -order method.

Initial conditions

As proposed by [START_REF] Peters | Requirements for computational methods to be used for the IMO second generation intact stability criteria[END_REF] and the new regulation (SDC 3/WP.5), the initial roll angle is equal to 5 degrees, and the initial roll speed is null.

Simulation time

The time duration of the simulation is equal to 15 times the ship's natural roll period, as also proposed by [START_REF] Peters | Requirements for computational methods to be used for the IMO second generation intact stability criteria[END_REF] and the new regulation (SDC 3/WP.5). Computations performed with 10 and 20 roll periods give exactly the same value of KGmax, which proves that a time duration of 15 roll periods is long enough to attain a steady-state roll amplitude. The simulation is interrupted if the roll angle exceeds 50 degrees.

Time step

We set the time step to one 40 th of the ship's natural roll period, and validated this value by removing the damping and GM variation from the differential equation to simulate an undamped roll on calm water. The roll amplitude remained equal to the initial roll angle. This test fails if the Euler method, rather than the Runge-Kutta method, is used to solve the differential equation (the roll amplitude increases). [START_REF] Peters | Requirements for computational methods to be used for the IMO second generation intact stability criteria[END_REF] and SDC 3/WP.5 propose a time step equal to one 30 th of the ship's natural roll period.

Simulation implementation

We performed the hydrostatic computations using Calcoque software [START_REF] Grinnaert | Calcoque: a fully 3D ship hydrostatic solver[END_REF], with the water density equal to 1.025 t/m 3 . We also solved the differential equation and calculated the damping coefficient with this software to make a user-friendly computation of the C2 coefficient and the associated KGmax curves.

KGmax computation

We determined the KGmax associated with the C2 coefficient in two steps.

In the first step, we found the first value of KG for which the C2 coefficient is lower than RPR0 and a second value of KG for which C2 is larger than RPR0. To find this interval, we began with a low value of KG for which C2 is lower than RPR0, near zero. We called this KG value the "KG starting value," denoted as "KGstart". We then increased the value of KG by 10 centimeters in each iteration until C2 became higher than RPR0.

The second step is based on a classical dichotomy to find the highest value of KG that corresponds to a C2 value lower than RPR0, at the required precision of 0.1 mm. Since C2 is the sum of the weighting coefficients drawn from the IACS wave scatter diagram (number of occurrences of each wave in Table 2 divided by 100 000), it is not possible to let the value of C2 be equal to RPR0. Table 7 shows the values of KG and C2 while finding the KGmax for the naval ship in this study, at a draft of 6.125 m. The associated KGmax is 9.2236 m (in gray in).

Results

We computed the KGmax curves associated with the parametric roll levels 1 and 2 criteria for the five ships in this study, as shown in Fig. 9 to Fig. 13. The KGmax curve determined by the first level-1 method (parallel water-plane) is marked "PR Level 11st method" (square markers), and that by the second level-1 method (ship balanced in trim and sinkage on a wave with the same length) is marked "PR Level 12nd method" (round markers). The KGmax curve determined by the level-2 C1 criterion is marked "PR Level 2C1" (solid line with diamond markers). The KGmax curve determined by the level-2 C2 criterion is marked "PR Level 2-C2" (dashed line with diamond markers). Although coefficients C1 and C2 use the same criterion, they are considered as separate criteria in this paper. Based on our results, we can make following observations, some of which are similar to those observed for the pure loss of stability failure mode:

1) For the naval ship, the results of the two level-1 methods differ significantly.

2) The KGmax value determined by the second level-1 method is almost equal to that given by level-2 C1 for the naval ship and is close to that of the tanker.

3) For the Ro-Ro vessel, the level-2 C1 is somewhat more conservative than the second level-1 method.

4) Both first-level methods yield the same KGmax for the tanker at full load displacement and beyond. This result is as expected due to the cylindrical shape of the hull.

5) For the tanker, the criterion of the level-2 C1 is limited by the zero-GM at full load displacement and beyond (condition A, criterion C1).

6) The level-2 C2 curve of the tanker is fully coincident with that linked to KMT. This is due to the limit determined by the negative GMmoy (see subsection 4.2.3) and shows that the ship is not vulnerable to parametric roll.

7) Both curves given by the C1 and C2 criteria differ significantly for all ships (except for the tanker at full load displacement). The C2 criterion is always less conservative than C1.

8) The curve linked to the C2 criterion for the 319-m container ship has two peculiar jumps. A lesser jump occurs on the same curve for the C11 container ship (at full load displacement). This point is discussed in the next subsection. 9) According to these criteria, the naval ship and tanker are considered to be not vulnerable to parametric roll. The Ro-Ro vessel may also be considered to be not vulnerable because its highest KGmax curve (C2) has only one point below the curve given by the 1 st -generation criteria, which may be lower if the weather criterion are taken into account. According to these criteria, both container ships are considered to be vulnerable. This confirms what is already known for the C11 [START_REF] France | An investigation of head-sea parametric rolling and its influence on container lashing systems[END_REF].

Focus on C2 criterion

For any ship at any draft, the KGmax linked to the C2 criterion is defined as the highest value of KG for which the value of C2 is lower than RPR0 = 0.06. Thus, it is interesting to check the C2-versus-KG curve. These curves are shown in Fig. 14 to Fig. 18. For the three non-vulnerable ships (naval, Ro-Ro, tanker), we made the computations assuming a full load draft for an interval containing KGmax with a step of 2 centimeters. For both container ships, we made the computations assuming a larger interval of KG with a step of 1 centimeter, at a draft of 9 m (C11 container ship, Fig. 15) and 12 m (319 m container ship, Fig. 16). Fig. 18 shows the curve of the tanker, for which C2 is equal to 0 for all values of KG lower than KGmax (13.70 m) and is equal to 1 for all higher values. This shows that a parametric roll never occurs on this ship. The value of C2 is forced to 1 when the average value of GM on waves becomes negative (see subsection 4.2.3).

Fig. 14 shows the curve for the naval ship, where we observe a small interval of KG (centered approximatively at 8.70 m) in which the C2 coefficient is not null. This shows that a parametric roll occurs for some lightly weighted waves that correspond with the ship's natural roll period. For higher values of KG, C2 reduces to near zero and then rapidly increases to 1. While parametric roll occurs in these KG conditions, the average value of GM is near zero: the ship becomes statically unstable on waves. Fig. 17 shows the C2-versus-KG curve for the Ro-Ro vessel, from which we observe that the increasing part of the curve is longer than both the tanker and naval ship. We also observe that two values of KG larger than KGmax (12.57 m) give values of C2 that are lower than RPR0 (KG = 12.60 and 12.62 m, marked with * on Fig. 17).

Fig. 15 and Fig. 16 show the same curves, respectively, for the C11 and 319 m container ships. On both, we observe many peaks and relatively large intervals of KG larger than KGmax, for which the value of C2 is lower than RPR0, thus fulfilling the associated criterion. These intervals are colored in gray on the corresponding figures. This non-monotonically increasing configuration of the C2 curve makes the starting value of KG (KGstart in this paper, 15 m for both container ships) very important in the search for KGmax. The value of the increment used in this search (10 centimeters in this study) is also very important. Both parameters must be chosen such that they avoid overlooking a small zone of KG for which C2 is larger than RPR0. The computation of the C2-versus-KG curves for all drafts makes it possible to define authorized and restricted areas according to the C2 criterion. Fig. 19 and Fig. 20 show these zones for both container ships (authorized areas in white, restricted areas in gray). The lower envelopes of the restricted areas give the KGmax curves linked to C2 that are shown in Fig. 10 and Fig. 11. The jumps observed on these curves correspond to the passage from the upper to the lower restricted area (at drafts equal to 10 and 13 m for the 319m container ship). In this study, we computed the KGmax curves associated with the level-1 and level-2 criteria of pure loss of stability and parametric roll failure modes for five ships of different types. The results reveal that level-2 criteria can be more conservative than the second level-1 method for both failure modes. This configuration is not anticipated in the future regulation. The results also show that the discrepancy between the level-1 methods can be very large, especially for the pure loss of stability failure mode. We argue that the safety margins ensured by the first level-1 method (parallel water-plane) seem to be excessive and may conflict with the future excessive acceleration criteria. Thus, if equipped to make hydrostatic computations on waves, naval architects and shipyards should not use the first level-1 method. However, this method could have real added value if it ensures an acceptable safety margin for small vessels, which are likely to be designed by architects of modest means.

The computation results of the second check of the parametric roll level-2 criterion (C2) for all possible values of KG shows that the more vulnerable the ship is to this failure mode, the more her C2-versus-KG curve is characterized by peaks and a long interval of C2 going from 0 to 1. This makes the surface formed by both the draft and KG foliated with both authorized and restricted zones. Consequently, it is not sufficient to check that C2 is lower than RPR0 for a given ship's loading condition. It is necessary to check this condition for all lower values of KG. We recommend that the future regulation clearly specify the requirement for this additional check.

The KGmax curves associated with the parametric roll criteria confirm the well-known vulnerability of the C11 container ship to this failure mode.

Fig. 1

 1 Fig. 1 Future intact stability regulation (Bulian and Francescutto, 2011)

Fig. 2

 2 Fig. 2 Parallel water-plane at the lowest draft dL

Fig. 3

 3 Fig. 3 Righting arm curve on static wave with angle of stable equilibrium under heeling lever RPL3 (S) and angle of vanishing stability (V)

Fig. 4 Fig. 8

 48 Fig. 4 KGmax curves associated with the pure loss of stability criteria for the naval ship

Fig

  Fig. 9 KGmax curves associated with the parametric roll criteria for the naval ship

Fig. 12

 12 Fig. 12 KGmax curves associated with the parametric roll criteria for the Ro-Ro vessel

Fig. 14 FigFig

 14 Fig. 14 Curve C2 versus KG for the naval ship (d = 6.125 m)

Table 1 Main particulars of ships

 1

	Particular	Symbol Unit	Naval ship Container ship C11 Container ship 319 m Ro-Ro Tanker
	Length overall	LOA	m	152.9	275.8	335.5	147.9	236.5
	Length between perpendiculars LPP	m	142	262	319	135	227.5
	Breadth	B	m	19.06	40	42.8	24.2	32.2
	Draft, full load	dfull	m	6.15	12	13	5.50	11
	Freeboard, full load	f	m	3.85	12.45	11.60	12.50	7.00
	Displacement, full load		t	8 634	73 340	107 350	11 544 70 397
	Speed	Vs	kts	30	20	25	20	15
	Froude number	Fn	-	0.413	0.203	0.230	0.283	0.163
	Block coefficient	Cb	-	0.507	0.569	0.590	0.625	0.852
	Bilge keels length	Lbk	m	35.7	76.53	81	45	75
	Bilge keels breadth	Bbk	m	0.55	0.40	0.50	0.30	0.30
	Bilge keels projected area	Ak	m 2	34.009	58.023	57.276•	19.092 31.82
	Added mass coefficient (*)	a	-	0.1	0.1	0.1	0.094	0.1
	Radius of inertia coefficient (*)	k	-	0.4	0.429	0.4	0.41	0.4
	Reference height of G	KGref	m	7.998	18.239	18	12.215 12
	Metacentric height @ KGref	GMref	m	1.5	1.973	2.174	1.638	1.7
	Roll period @ KGref	T0	s	14.12	25.7	24.43	16.27	20.78
	Starting value of KG (*)	KGstart m	7.00	14.00	15.00	9.00	10.00
	Pure loss limit angle (*)	RPL2	-	25	25	25	25	25
	GM variation limit (*)	RPR	-	0.418	0.405	0.348	0.418	0.355
	Roll damping coefficient	B44	N•m•s/rad4.018e7	1.027e9	1.099e9	1.721e8 2.463e8
	Superstructures	-	-	Yes	Yes	Yes	No	Yes

Data marked with (*) are defined later in this paper.

Table 2 .

 2 This table lists 16 wave periods and 17 wave heights, and 197 waves have a non-zero number of occurrence. The new regulation allows for the use of another wave scattering table if the ship is sailing in a restricted area.

Table 2 Wave scatter diagram (from IACS, 2001)

 2

					Number of occurrences: 100000 / TZ(s)=average zero up-crossing wave period
	HS /m	3.5	4.5	5.5	6.5	7.5	8.5	9.5	10.5	11.5	12.5 13.5 14.5 15.5 16.5 17.5 18.5
	0.5	1.3 133.7 865.6	1186	634.2	186.3	36.9	5.6	0.7	0.1	0	0	0	0	0	0
	1.5	0	29.3	986	4976	7738	5569.7 2375.7 703.5	160.7	30.5	5.1	0.8	0.1	0	0	0
	2.5	0	2.2	197.5 2158.8	6230	7449.5 4860.4	2066	644.5	160.2 33.7 6.3	1.1	0.2	0	0
	3.5	0	0.2	34.9	695.5 3226.5	5675	5099.1	2838	1114.1 337.7 84.3 18.2 3.5	0.6	0.1	0
	4.5	0	0	6	196.1 1354.3 3288.5 3857.5 2685.5 1275.2 455.1 130.9 31.9 6.9	1.3	0.2	0
	5.5	0	0	1	51	498.4 1602.9 2372.7 2008.3	1126	463.6 150.9 41	9.7	2.1	0.4	0.1
	6.5	0	0	0.2	12.6	167	690.3 1257.9 1268.6 825.9	386.8 140.8 42.2 10.9 2.5	0.5	0.1
	7.5	0	0	0	3	52.1	270.1	594.4	703.2	524.9 276.7 111.7 36.7 10.2 2.5	0.6	0.1
	8.5	0	0	0	0.7	15.4	97.9	255.9	350.6	296.9	174.6 77.6 27.7 8.4	2.2	0.5	0.1
	9.5	0	0	0	0.2	4.3	33.2	101.9	159.9	152.2	99.2 48.3 18.7 6.1	1.7	0.4	0.1
	10.5	0	0	0	0	1.2	10.7	37.9	67.5	71.7	51.5 27.3 11.4	4	1.2	0.3	0.1
	11.5	0	0	0	0	0.3	3.3	13.3	26.6	31.4	24.7 14.2 6.4	2.4	0.7	0.2	0.1
	12.5	0	0	0	0	0.1	1	4.4	9.9	12.8	11	6.8	3.3	1.3	0.4	0.1	0
	13.5	0	0	0	0	0	0.3	1.4	3.5	5	4.6	3.1	1.6	0.7	0.2	0.1	0
	14.5	0	0	0	0	0	0.1	0.4	1.2	1.8	1.8	1.3	0.7	0.3	0.1	0	0
	15.5	0	0	0	0	0	0	0.1	0.4	0.6	0.7	0.5	0.3	0.1	0.1	0	0
	16.5	0	0	0	0	0	0	0	0.1	0.2	0.2	0.2	0.1	0.1	0	0	0

Table 3 Effective height of the top 3% highest waves (m) for a zero-crossing period of 10.5 seconds

 3

	HS / m	Naval	Container	Container	Ro-Ro	Tanker
		L=142m	L=262m	L=319m	L=135m	L=227.5m
	0.5	0.430	0.539	0.538	0.417	0.526
	1.5	1.289	1.618	1.614	1.251	1.579
	2.5	2.149	2.697	2.691	2.085	2.631
	3.5	3.008	3.776	3.767	2.919	3.683
	4.5	3.868	4.855	4.843	3.753	4.736
	5.5	4.727	5.934	5.919	4.587	5.788
	6.5	5.587	7.013	6.996	5.421	6.841
	7.5	6.446	8.092	8.072	6.255	7.893
	8.5	7.306	9.171	9.148	7.089	8.945
	9.5	8.165	10.250	10.225	7.923	9.998
	10.5	9.025	11.329	11.301	8.757	11.050
	11.5	9.884	12.408	12.377	9.591	12.103
	12.5	10.744	13.487	13.453	10.425	13.155
	13.5	11.603	14.566	14.530	11.258	14.207
	14.5	12.463	15.645	15.60	12.092	15.260
	15.5	13.322	16.724	16.682	12.926	16.312
	16.5	14.182	17.803	17.758	13.760	17.365

Table 4 Wave cases for C1 coefficient of parametric roll level 2

 4

	Wave	Weight Wi	Wave length	Wave height
	number i		i / m	Hi / m
	1	0.000013	22.574	0.350
	2	0.001654	37.316	0.495
	3	0.020912	55.743	0.857
	4	0.092799	77.857	1.295
	5	0.199218	103.655	1.732
	6	0.248788	133.139	2.205
	7	0.208699	166.309	2.697
	8	0.128984	203.164	3.176
	9	0.062446	243.705	3.625
	10	0.024790	287.931	4.040
	11	0.008367	335.843	4.421
	12	0.002473	387.440	4.769
	13	0.000658	442.723	5.097
	14	0.000158	501.691	5.370
	15	0.000034	564.345	5.621
	16	0.000007	630.684	5.950
	Condition A is as follows:

Table 5 Speed factors Kj

 5

	Number j Speed factor Kj
	1	-1
	2	-0.866
	3	-0.5
	4	0
	5	0.5
	6	0.866
	7	1

Table 6 Effective height of the top 1/3 highest waves (m) for a zero-crossing period equal to 10.5 seconds

 6

	HS /m	Naval	Container	Container	Ro-Ro	Tanker
		L=142m	L=262m	L=319m	L=135m	L=227.5m
	0.5	0.288	0.362		0.361		0.280	0.353
	1.5	0.865	1.085		1.082		0.839	1.059
	2.5	1.442	1.809		1.804		1.399	1.764
	3.5	2.018	2.532		2.525		1.958	2.470
	4.5	2.595	3.256		3.247		2.518	3.176
	5.5	3.171	3.979		3.968		3.077	3.882
	6.5	3.748	4.703		4.690		3.637	4.587
	7.5	4.325	5.426		5.411		4.196	5.293
	8.5	4.901	6.149		6.133		4.756	5.999
	9.5	5.478	6.873		6.854		5.315	6.705
	10.5	6.054	7.596		7.576		5.875	7.410
	11.5	6.631	8.320		8.297		6.434	8.116
	12.5	7.208	9.043		9.019		6.993	8.822
	13.5	7.784	9.767		9.740		7.553	9.528
	14.5	8.361	10.490		10.462	8.112	10.234
	15.5	8.937	11.214		11.183	8.672	10.939
	16.5	9.514	11.937		11.905	9.231	11.645
	The ship's natural roll period in calm water (T0) is linked
	to her moment of inertia, as follows:
		0 T		2π	44 ΔGM J g		2π k B	1 GM a g 	(29)
	with:
	g	9.81 m/s 2 gravitational acceleration;
	GM (m)

Table 7 Finding the value of KGmax

 7

	First step: finding the interval
	Iteration	KG (m)	C2	C2<RPR0
	1	KGstart = 7	0	Y
	2	7.1	0	Y
	3	7.2	0	Y
	4	7.3	0	Y
	5	7.4	0	Y
	6	7.5	0	Y
	7	7.6	0	Y
	8	7.7	0	Y
	9	7.8	2.86E-07	Y
	10	7.9	1.57E-06	Y
	11	8.0	4.29E-07	Y
	12	8.1	0	Y
	13	8.2	0	Y
	14	8.3	0	Y
	15	8.4	0	Y
	16	8.5	0	Y
	17	8.6	2.49E-05	Y
	18	8.7	0.0293483	Y
	19	8.8	0.017636	Y
	20	8.9	0.00266057	Y
	21	9.0	0.000234857	Y
	22	9.1	3.13E-05	Y
	23	9.2	0.0219234	Y
	24	9.3	0.154703	N
	Second step: dichotomy
	Iteration	KG (m)	C2	C2<RPR0
	25	9.25	0.0760626	N
	26	9.225	0.0749527	N
	27	9.2125	0.0371166	Y
	28	9.21875	0.0449244	Y
	29	9.22188	0.0462103	Y
	30	9.22344	0.0464114	Y
	31	9.22422	0.07403	N
	32	9.22383	0.0740291	N
	33	9.22363	0.0502479	Y
	34	9.22373	0.0740291	N

Acknowledgements

The authors would like to acknowledge Clève Wandji (Bureau Veritas), Philippe Corrignan (Bureau Veritas), and Jean-François Leguen (DGA/TH) for their valuable discussions during the preparation of this paper.