Myline Cottance
email: m.cottance@esiee.fr

Sebastien Nazeer
email: s.nazeer@esiee.fr

Rousseau Lionel

Gaelle Bazin Lissorgues
email: g.lissorgues@esiee.fr

Alexandre Bongrain
email: alexandre.bongrain@cea.fr

Raphael Kiran
email: raphael.kiran@cea.fr

Emmanuel Scorsone
email: emmanuel.scorsone@cea.fr

Philippe Bergonzo
email: philippe.bergonzo@cea.fr

Amel Bendali
email: amel.bendali@inserm.fr

Serge Picaud
email: serge.picaud@inserm.fr

Sébastien Nazeer

Lionel Rousseau
email: l.rousseau@esiee.fr

Gaëlle Lissorgues

Raphaël Kiran

Sébastien Joucla
email: sebastien.joucla@u-bordeaux1.fr

Blaise Yvert
email: blaise.yvert@u-bordeaux1.fr

Diamond micro-electrode arrays (MEAs): A new route for in-vitro applications

Keywords: microelectrode arrays, electrical stimulation, recording, neural network, diamond electrode

I NTRODUCTION Micro-Electrode Arrays (MEAs) were developed since the early 1970s and play an important role in the comprehension of the mechanisms of electrogenic cell biological activity. They are widely used in various domains especially in neuroscience as they correspond to noninvasive techniques that allow both long-term stimulation and recording of bio-electrical signals of cultured cells and tissues [START_REF] Morin | Investigating Neural Activity with Planar Microelectrode Arrays: Achievements and New Perspectives[END_REF]. For instance, MEAs can help to understand some currently misunderstood phenomena such as extracellular electrical stimulation mechanisms that excite nervous cells. In clinical research, extracellular electrical stimulation of the central nervous system is used in therapeutic aims to restore different sensory or motor functions (retinal [START_REF] Meyer | Retina implant -a bioMEMS challenge[END_REF] and cochlear implants [START_REF] Wilson | Cochlear implants: Current designs and future possibilities[END_REF], lesioned spinal cord stimulation [START_REF] Grabow | Spinal Cord Stimulation for Complex Regional Pain Syndrome: An evidence-based Medicine Review of the literature[END_REF], deep brain stimulation [START_REF] Voon | Deep Brain Stimulation: Neuropsychological and Neuropsychiatric Issues[END_REF]). In fundamental research, this type of stimulation is used as a tool to investigate network of neurons to highlight their organization, connectivity, and also to study their modification in response to external stimulation. MEAs also enable extracellular electrical recordings which allow gaining of information about interactions between electrogenic cells, sampling distribution of electrophysiological behavior and/or monitoring changes of electrical activity which are suitable for drug screening and discovery applications [START_REF] Stett | Biological application of microelectrode arrays in drug discovery and basic research[END_REF].

Here, we report on the fabrication of diamond planar MEAs of two matrix geometries, 8×8 and 4×15, designed for the stimulation and recording of respectively rat retina and mouse embryonic spinal cord. The material electrode chosen in this study is the boron-doped nanocrystalline diamond (BNCD) as it exhibits many valued advantages due to its physical and chemical properties [START_REF] Fraft | Doped Diamond: A Compact Review on a New, Versatile Electrode Material[END_REF] compared with common platinum and gold electrodes. Besides its 3V electrochemical water potential, BNCD also display biocompatible properties which render this material ideal for biological in-vitro applications. Moreover, the novelty of the work is related to an original process development since diamond growth is not obviously compatible with standard silicon micro-fabrication steps due to temperature constraints. Two designs of MEAs depicted in Fig. 1, were fabricated in order to suit the geometry of studied tissues for a better recording and stimulation of cells.

II. M ETHODS

A. Designs of MEAs

The first design consists of circular electrodes of 80 µm 100 µm diameter of unique size and 300×900 µm pitch arranged in a 4×15 matrix. It is adapted to the geometry of mouse embryonic hindbrain-spinal cord preparation with an active area of 11.34 mm².

The second design is made of circular electrodes of 34 µm diameter of unique size and 100 µm pitch arranged in an 8×8 matrix and it is adapted to the geometry of rat retina preparation with an active area of 0.49 mm². We have investigated previously a process where we deposited the metal layer for contacting the diamond before growing diamond on top of the metal. Because of high temperature required for the diamond growth (above 850°C), the metal exhibits a blistered and cracked surface resulting in poor adhesion to the substrate. As the steps order was mainly governed by diamond growth temperature, a new fabrication process of diamond MEAs is developed (Fig. 2) where diamond is grown before metal deposition. First, detonation diamond nanoparticles were spin coated on a 4 inches silicon or glass wafer using the process described in [START_REF] Scorsone | Enhanced control of diamond nanoparticle seeding using a polymer matrix[END_REF]. An aluminum hard mask was then sputter deposited to define electrode shape by lithography. Diamond nanoparticles outside these protected areas were etched away using a short oxygen reactive ion etch (RIE) plasma. Next, the metal hard mask was removed chemically to reveal the diamond particles patterns which thus were exposed to a gas mixture of methane, hydrogen and trimethylboron in a MPECVD (Microwave Plasma Enhanced Chemical Vapor Deposition) diamond growth reactor (Seki technotron corp. AX6500X). The grown electrodes were individually contacted by the sputtered deposition of Ti (50 nm)/ Pt (150 nm) metal tracks using the lift-off process and Nlof 2020 as a sacrificial photoresist layer. A 600 nm silicon nitride (Si 3 N 4 ) layer was deposited on the substrate by plasma-enhanced chemical vapor deposition in order to insulate the metal tracks from the electrolyte solution. Then, a sulfur hexafluoride (SF 6 ) RIE was performed to etch locally the silicon nitride to open contact pads and diamond electrodes, which also defines the electrode size. Finally, the photoresist used to selectively etch the nitride layer was removed and the MEAs were diced using a diamond saw, being ready to be used for biology. The grown diamond electrodes, characterized using a UV-Vis Spectroscopic Ellipsometer (Horiba Jobin-Yvon, UVISEL), present a thickness of approximately 300 nm over a diameter of 80 µm for the 4×15 MEA. After the passivation step which defines the active area, electrodes exhibit a smooth surface aspect and metal annular contacts and tracks are well defined as shown in Fig. 3. The morphology and grain size were observed using scanning electron microscopy (SEM). The typical morphology of a boron-doped diamond layer is shown in Fig. 4. The diamond crystals are highly faceted, with an average grain size of 100 nm and the obtained diamond films are uniform and continuous. The measured boron concentration determined with a secondary ion mass spectrometry (SIMS) depth profiling using a Cameca IMS 4f ion microprobe was approximately 2×10 21 at.cm -3 .

B. Fabrication process of MEAs

III. R ESULTS

A. Morphological characterization

B. Electrochemical characterization

Cyclic voltammetry (CV) provides information on the potential window of the electrode which is the voltage range between which the material electrode does not induce water oxidation nor reduction. This window is important for the stability of cells and tissues facing the electrodes to prevent them to be exposed to ionic species nor hostile pH. 0.1 V.s -1 . This experiment demonstrates that diamond electrodes have a wider window potential (3.5 V) than platinum and gold electrodes (2.5 V) reported in literature [START_REF] Kraft | Doped Diamond: A Compact Review on a New , Versatile Electrode Material[END_REF][START_REF] Nebel | Diamond for bio-sensor applications[END_REF].

Fig. 6 presents the CV of the same electrode in a ferri/ferrocyanide solution where the peak to peak separation of the oxidation and reduction potential was probed at 63 mV at a 0.1 V.s -1 scan rate. This experiment demonstrates the extreme reactivity of the diamond surface.

This wide potential window of diamond can be explained by the inertness of its surface for adsorption (hydrogen termination) [START_REF] Alexander | The influence of surface termination on the electrochemical properties of boron doped diamond[END_REF]. Electrochemical impedance spectroscopy (EIS) was also performed to study the impedance of the electrode. Measurements were made at 1 kHz (sine wave) thus corresponding to characteristic frequency of neural biological activity. They are performed using a NanoZ device (Fig. 7) from Multichannel Systems (Reutlingen, Germany) with an Ag/AgCl electrode in phosphate buffered saline (PBS) aqueous solution of pH 7.2. The equipment gives the magnitude and phase of the impedance value for each electrode.

We found at 1 kHz that impedance values were within 200-300 kΩ for 80-µm planar electrodes. C. Functional Tests Two sets of experiments are currently being performed with the fabricated MEAs: the first one is the extracellular signal recording and electrical stimulation of an embryonic mouse hindbrain-spinal cord and the second one is the recording of spontaneous extracellular signals of a retina.

Preliminary results obtained on a whole embryonic hindbrain-spinal cord (as in [START_REF] Yvert | Artificial CSF motion ensures rhythmic activity in the developing CNS ex vivo : A mechanical source of rhythmogenesis?[END_REF]) show that BNCD MEAs could successfully detect spontaneous activity composed of spike bursts with and amplitude in the range of a few µV and lasting for a few seconds (see Fig. 8).

IV. C ONCLUSION

We have presented a novel microfabrication process of diamond planar MEAs, assessing that the combination of standard silicon micro-fabrication processes with diamond growth is possible to fabricate functional devices for neural recordings.

The manufactured MEAs were characterized with cyclic voltammetry demonstrating that MEAs were successfully fabricated. The wider potential window of diamond compared to other material electrode is an important attractive feature for sensor applications [START_REF] Nebel | Diamond for bio-sensor applications[END_REF].

In-vitro electrical recordings of spontaneous activity of an embryonic mouse hindbrain-spinal cord with our 4×15 MEA and of a retina rat with our fabricated 8×8 MEA were successfully accomplished. Extracellular electrical neural stimulation using MEAs is the purpose of current development.

New materials or combination of materials will be also studied to improve our MEAs properties. Moreover, as our MEAs and diamond material exhibit promising results, they also led us to head towards in-vivo applications such as working on the development of retinal implants with similar diamond electrodes.

Fig. 1 .

 1 Fig. 1. On the left, electrodes arranged in a 4×15 matrix and on the right electrodes arranged in an 8×8 matrix.

Fig. 2 .

 2 Fig. 2. Process flow of diamond MEAs fabrication process

Fig. 3 .

 3 Fig. 3. Optical microscope images of a 80 µm diamond electrode of a 4×15 MEA before passivation (on the left), and a 34 µm diamond electrode of a 8×8 MEA after passivation (on the right).

Fig. 4 .

 4 Fig. 4. SEM image of BNCD film

Fig. 5

 5 presents the result of CV characteristics as probed on our manufactured diamond MEA electrode (80-µm diameter) in a 0.5 M aqueous LiClO 4 solution scanned at

Fig. 5 .Fig. 6 .

 56 Fig. 5. CV of diamond electrodes on a 0.5 M aqueous LiClO4 solution scanned at 0.1 V.s -1 .

Fig. 7 .

 7 Fig. 7. A 8×8 MEA fixed on NanoZ device to measure the magnitude and phase of the electrode impedance at 1 kHz.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Recording of a spike burst with a 4×15 matrix diamond MEA.Extracellular recordings of spontaneous activity of an ex vivo rat retina on an 8×8 MEA are shown in Fig.8. Pieces of retina were placed on the array with the ganglion cells (spiking neurons in the retina) facing the electrodes. Spontaneous activities of a retina were recorded on several electrodes of which two recordings of two different electrodes can be seen on Fig.9(A) and an enlargement of the recording exhibits on Fig.9(B), a typical shape of a recorded action potential.

ACKNOWLEDGMENT

The authors would like to thank the French National Research Agency for their financial support through the MEDINAS project as well as the European Community's Seventh Framework Programme (FP7) under grant agreement n° 280433 (NEUROCARE project).

BIOGRAPHY

Myline Cottance got her bachelor's degree in organic chemistry in 2010 from the Science Faculty of Orsay and her master's degree in chemical engineering in 2011 from the University of Pierre and Marie Curie. She currently conducts PhD studies at ESIEE-Paris working on the fabrication of MEAs for retinal implants and neuroelectronic interfaces.


M. D T H C Le

' Ottavio

P Vidal

O Polit

A new robust quadrilateral four-node variable kinematics plate element for composite structures

that define the QC4 interpolation for

Introduction

The increasing use of composite laminates and sandwich structures in engineering applications drives the need for appropriate analysis and design tools with dedicated computational models. Based on geometric considerations, composite panels are conveniently modeled as two-dimensional plate/shell structures. However, complicating effects -such as anisotropy, heterogeneity and transverse shear compliance -call for plate/shell models that go beyond the so-called classical models, i.e., those relying on i ( ) Kirchhoff-Love assumptions and neglecting transverse deformation (Classical Laminate Plate Theory, CLPT), ii ( ) or on Reissner-Mindlin assumptions and retaining a merely constant transverse shear deformation through the thickness (First-order Shear Deformation Theory, FSDT). Numerous review papers devoted to high-order plate/shell models witness the scientific progress in this specific topic [START_REF] Noor | Assessment of shear deformation theories for multilayered composite plates[END_REF][START_REF] Reddy | Theories and computational models for composite laminates[END_REF][START_REF] Noor | Computational models for sandwich panels and shells[END_REF][START_REF] Carrera | Theories and finite elements for multilayered, anisotropic, composite plates and shells[END_REF][START_REF] Sayyad | On the free vibration analysis of laminated composite and sandwich plates: a review of recent literature with some numerical results[END_REF].

A useful classification discerns Equivalent Single Layer (ESL) and Layer-Wise (LW) models [START_REF] Reddy | Mechanics of Laminated Composite Plates and Shells: theory and Analysis[END_REF]: In the former model class, the number of unknowns is independent of the number of layers constituting the composite plate, while in the latter one, the number of unknowns increases with the number of constituting layers. CLPT and FSDT evidently pertain to the ESL models. ESL models that enhance the kinematics for the transverse shear deformation while still discarding the transverse normal deformation, are referred to as High-order Shear Deformation Theories (HSDT). Contrary to FSDT, no numerical shear correction factors are required in HSDT thanks to an at least parabolic transverse shear distribution that may also exactly verify the stress boundary conditions at the laminate's top and bottom surfaces. Seminal examples are the third-order theory by Reddy [START_REF] Reddy | A simple higher-order theory for laminated composite plates[END_REF] and the sinus-model by Touratier [START_REF] Touratier | An efficient standard plate theory[END_REF]. Conventional ESL models employ a single approximation for the displacement field through the entire laminate's thickness; the resulting transverse strain field is continuous across the stacked layers, which contradicts the equilibrium conditions between adjacent plies with different material properties [START_REF] Carrera | C 0 z-Requirements -models for the two dimensional analysis of multilayered structures[END_REF]. Zig-Zag models are special ESL models with a piece-wise continuous displacement field that allows to fulfill the interlaminar continuity conditions of the transverse shear stresses, see, e.g., [START_REF] Vidal | A refined sinus plate finite element for laminated and sandwich structures under mechanical and thermomechanical loads[END_REF][START_REF] Tessler | Refined zigzag theory for homogeneous, laminated composite, and sandwich beams derived from Reissner's mixed variational principle[END_REF]. The review paper by Carrera [START_REF] Carrera | Historical review of zig-zag theories for multilayered plates and shells[END_REF] offers a comprehensive discussion about the various approaches to Zig-Zag models. In presence of highly compliant core layers, or for accurately resolving local stress gradients whose wavelength is comparable to the plate thickness, it is necessary to further refine the model upon including the transverse normal deformation, see, e.g., the early contribution by Lo et al [START_REF] Lo | A high-order theory of plate deformation -Part 1: homogeneous plates -Part 2: laminated plates[END_REF] and the more recent papers [START_REF] Barut | C 0 -continuous triangular plate element for laminated composite and sandwich plates using the {2,2}-Refined Zigzag Theory[END_REF][START_REF] Vidal | Assessment of the refined sinus plate finite element: free edge effect and Meyer-Piening sandwich test[END_REF]. Retaining the full three-dimensional constitutive law is particularly useful for coupled problems, such as those related to thermo-mechanics [START_REF] Polit | Robust C 0 high-order plate finite element for thin to very thick structures: mechanical and thermo-mechanical analysis[END_REF].

In order to reduce the computational effort, it has been proposed to limit the use of expensive high-order models to those local regions

MARK

Variable kinematics plate model suffer transverse shear locking also, and the first employed countermeasures consisted in a selectively reduced quadrature [START_REF] Robbins | Modelling of thick composites using a layerwise laminate theory[END_REF][START_REF] Carrera | Classical and advanced multilayered plate elements based upon PVD and RMVT. Part 1: derivation of finite element matrices. Part 2: numerical implementations[END_REF][START_REF] Ottavio | Considerations on higher-order finite elements for multilayered plates based on a unified formulation[END_REF]. Rectangular four-and nine-nodes plate elements have been been subsequently implemented upon extending the MITC approach to CUF-based high-order kinematics [START_REF] Carrera | MITC technique extended to variable kinematic multilayered plate elements[END_REF]. A fournode quadrilateral plate element has been proposed by Kulikov and Plotnikova [42] by resorting to a hybrid-mixed ANS approach in conjunction with a variable kinematics approach formulated in terms of Sampling Surfaces (SaS). In these works, all transverse shear strain terms issued from the high-order kinematics are constrained according to the adopted MITC or ANS approach. However, since high-order shear deformation terms depend on the plate thickness and will vanish in thin-plate limit, the locking behavior is produced by the first-order Reissner-Mindlin kinematics only. As a matter of fact, the convergence rates of CUF elements do not depend on the polynomial order N defining the plate kinematics [START_REF] Ottavio | Considerations on higher-order finite elements for multilayered plates based on a unified formulation[END_REF].

Refined two-dimensional structural models including thickness change require special care for preventing thickness or Poisson locking, i.e., an overly stiff bending response that is produced if the transverse normal strain is not allowed to vary along the thickness direction [START_REF] Sze | Three-dimensional continuum finite element models for plate/shell analysis[END_REF]. Two class of remedies have been proposed for correctly resolving the Poisson coupling along the thickness: to retain an at least linearly varying transverse normal strain, either within an Enhanced Assumed Strain approach pioneered by Büchter et al [START_REF] Büchter | Three-dimensional extension of nonlinear shell formulation based on the enhanced assumed strain concept[END_REF], see also [START_REF] Naceur | On the modeling and design of composite multilayered structures using solid-shell finite element model[END_REF], or directly in the model kinematics [START_REF] Polit | Robust C 0 high-order plate finite element for thin to very thick structures: mechanical and thermo-mechanical analysis[END_REF]; or to modify the constitutive law for the bending contribution by referring to a generalized plane stress condition, as proposed, e.g., in [START_REF] Kulikov | Non-linear geometrically exact assumed stressstrain four-node solid-shell element with high coarse-mesh accuracy[END_REF].

Based on the outlined background, the present paper proposes a four-node quadrilateral element for variable kinematics displacementbased CUF plate models. A special transverse shear locking correction is formulated by referring to the field consistency paradigm and applied only to the constant, thickness-independent part of the transverse shear strain. For this, the method first proposed by Polit et al [START_REF] Polit | A new eight-node quadrilateral shear-bending plate finite element[END_REF] for FSDT, and subsequently extended to a refined kinematics [START_REF] Polit | Robust C 0 high-order plate finite element for thin to very thick structures: mechanical and thermo-mechanical analysis[END_REF], is here further extended to high-order plate models with arbitrary kinematics. The resulting plate FE is implemented as User Element in the commercial package Abaqus, along with dedicated Python plug-ins for generating the model within the graphical interface of Abaqus. Following the recommendations expressed by MacNeal and Harder [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF], extensive numerical tests are reported for demonstrating the absence of spurious mechanisms and assessing the element's accuracy for extreme thickness ratios as well as distorted meshes. The paper is organized as follows: the CUF-based variable kinematics approach is recalled in Section 2 and the new QC4 FE approximation is presented in Section 3. The numerical results are discussed in Section 4, where a comprehensive investigation is proposed that concerns the rank of the stiffness matrix, the robustness of the element with respect to lengthto-thickness ratio and mesh distortion, as well as the accuracy of the predicted displacements and stresses for homogeneous and composite plates. Finally, Section 5 summarizes the main conclusions and proposes an outlook towards further studies. . hosting the stress gradients of interest, while employing classical loworder models for the remaining large portions of the structure. Several approaches have been followed for coupling heterogeneous kinematics and their respective Finite Element (FE) approximations, see, e.g., the transition element proposed by Feng and Hoa [START_REF] Feng | Partial hybrid finite elements for composite laminates[END_REF] within their hybridmixed composite plate element formulations, the overlapping mesh approach based on Arlequin method [START_REF] Hu | Multi-scale modelling of sandwich structures using the Arlequin method -Part I: linear modelling[END_REF], or the direct interface coupling proposed by Wenzel et al [START_REF] Wenzel | Coupling of heterogeneous kinematics and finite element approximations applied to composite beam structures[END_REF] on the basis of an extended variational formulation.

A very flexible manner for implementing such substructuring methods is to resort to so-called variable kinematics models and corresponding finite elements, as first pioneered by Reddy [START_REF] Reddy | An evaluation of equivalent-single-layer and layerwise theories of composite laminates[END_REF] and more recently systematically developed by Carrera and co-workers thanks to a dedicated Unified Formulation [START_REF] Biscani | Variable kinematic plate elements coupled via Arlequin method[END_REF][START_REF] Carrera | Use of Lagrange multipliers to combine 1D variable kinematic finite elements[END_REF]. By virtue of an extensive index notation, Carrera's Unified Formulation (CUF) allows to implement a large number of finite plate/shell elements within a unique software, whose inputs provided at runtime define the actual plate/shell model to be used in the analysis [START_REF] Carrera | Finite Element Analysis of Structures through Unified Formulation[END_REF]. Within CUF, classical displacement-based plate elements are defined through the choice of i() the description at laminate level (ESL or LW), (ii) whether a Zig-Zag function is to be used or not, and (iii) the order N of the polynomial expansion assumed for the displacement field. The original version of CUF as presented by [START_REF] Carrera | Classical and advanced multilayered plate elements based upon PVD and RMVT. Part 1: derivation of finite element matrices. Part 2: numerical implementations[END_REF] has been implemented as Abaqus User Element with a dedicated Python-based preprocessing tool [START_REF] Ottavio | Damage and its evolution in fiber-composite materials: simulation and non-destructive evaluation[END_REF]; this package has been successfully employed in, e.g., [START_REF] D'ottavio | Assessment of plate theories for free-edge effects[END_REF][START_REF] Wenzel | Assessment of free-edge singularities in composite laminates using higher-order plate elements[END_REF] and forms the basis of the developments proposed in this paper.

From the computational point of view, the development of robust Finite Elements (FE) is required to cope with the adopted twodimensional plate/shell model, see, e.g., the discussion by MacNeal [START_REF]Perspective on finite elements for shell analysis[END_REF] about the FE technology employed for Kirchhoff-Love and Reissner-Mindlin shell models. A general, highly predictive plate FE should not rely on numerical tuning coefficients, should have only six rigid body modes, for spurious zero-energy modes could be particularly detrimental, e.g., in non-linear analyses, and should be free from numerical pathologies that could degrade the solution's accuracy in case of distorted elements or extreme thickness ratios. The most characteristic example for this latter issue is transverse shear locking, a spurious over-constraint that dramatically underestimates the bending deformation of a thin, shear-deformable plate element. Several techniques have been devised for correcting the transverse shear locking pathology affecting FSDT-based plate/shell elements, most of which can be stated from hybrid-mixed approaches [START_REF] Pian | State-of-the-art development of hybrid/mixed finite element method[END_REF]. The most widespread techniques are reduced integration methods, which, however, require a dedicated stabilization for preventing spurious zeroenergy modes [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF], or so-called B-bar methods [START_REF] Hughes | The Finite Element Method[END_REF], in which a specific constraint is used for the transverse shear strain field. Different approaches have been followed for constructing this modified strain field, such as Kirchhoff mode [START_REF] Hughes | Finite elements based upon Mindlin plate theory with particular reference to the four node bilinear isoparametric element[END_REF], Assumed Natural Strain (ANS) [START_REF]Derivation of element stiffness matrices by assumed strain distributions[END_REF][START_REF] Park | The ANS shell elements: earlier developments and recent improvements[END_REF], Mixed Interpolation of Tensorial Components (MITC) [START_REF] Bathe | A four-node plate bending element based on Mindlin/ Reissner plate theory and a mixed interpolation[END_REF], the field-consistency paradigm [START_REF] Somashekar | A field-consistent, four-noded, laminated, anisotropic plate/shell element[END_REF], Discrete Shear elements [START_REF] Batoz | A discrete shear triangular nine d.o.f. element for the analysis of thick to very thin plates[END_REF] or Discrete Shear Gap (DSG) [START_REF] Bletzinger | A unified approach for shear locking free triangular and rectangular shell finite elements[END_REF].

Variable kinematics plate model

Let us consider a multilayered plate occupying the domain

V Ω x = × {-≤ ≤ } e e 2 3 2
in a Cartesian coordinate system x y z x ( , , ) = ( ) i , see Fig. 1. Unless otherwise stated, Latin indices range in {1, 2, 3}, Greek indices range in {1, 2} and tensorial repeated index convention is employed. Ω is the reference surface of arbitrary shape lying in the x x ( , ) -

1 2
plane located for convenience at z=0. The plate has constant thickness e, which is composed of k N = 1, 2, … L ortho- tropic, elastic and perfectly bonded layers, each with a thickness e k

( )
and with an orientation of the material orthotropy axes defined by the rotation angle θ k ( ) about the thickness direction x z ≡ 3 .

The weak form of the boundary value problem

The Principle of Virtual Displacement (PVD) is employed to obtain the weak form for a displacement-based plate approximation. The volume boundary is split as

V Γ Γ ∂ = ⋃ u σ ,
∫ ∫ δu σ u x δu t x - ϵ ( ) ( ) d + d = 0 V ij l ij m i Γ i i α σ ( 1
)
where t i are the surface loads at Γ σ , a bar denoting prescribed values; δu ϵ ( ) ij l is the compatible virtual strain tensor and σ u ( ) ij m the stress tensor defined by means of the linear elastic constitutive law in terms of the actual strains u ϵ ( ) rs m .

Variable kinematics assumptions

Carrera's Unified Formulation (CUF) is a technique that permits to implement a large number of bi-dimensional models in a unified manner by means of an extensive use of a compact index notation [START_REF] Carrera | Classical and advanced multilayered plate elements based upon PVD and RMVT. Part 1: derivation of finite element matrices. Part 2: numerical implementations[END_REF][START_REF] Carrera | Finite Element Analysis of Structures through Unified Formulation[END_REF]. Within the displacement-based approach, all plate theories are defined in CUF by specifying i ( ) whether the displacement field is described in Equivalent Single Layer (ESL) or Layer-Wise (LW) manner, and ii ( ) the order N of the polynomial expansion used for approximating the behavior along the thickness direction. The twodimensional variable kinematics plate model is generally formulated upon separating the in-plane variables x α from the thickness direction z, along which the displacement field is a priori postulated by known functions F(z):

u x z F z u x ( , ) = ( ) ( ), i α τ i α τ ( 2
)
where τ N = 0, 1,…, is the summation index and the order of expansion N is a free parameter of the formulation. In this work N can range from 1 to 4, in agreement with the classical CUF implementation [START_REF] Carrera | Classical and advanced multilayered plate elements based upon PVD and RMVT. Part 1: derivation of finite element matrices. Part 2: numerical implementations[END_REF].

In order to deal with both ESL and LW descriptions within a unique notation, it is convenient to refer to a layer-specific thickness coordinate z

z z ∈ { , } k b k t k ( ) ( )
that ranges between the z-coordinates of the bottom (subscript b) and top (subscript t) planes delimiting the k th layer, see Fig. 1. Eq. (2) can thus be formally re-written for each layer as

u x z F z u x F z u x F z u x ( , ) = ( ) ( ) + ( ) ( ) + ( ) ( ) i k α k t k i k α b k i k α r k i k α ( ) ( ) ( ) ( ) t b r (3)
with τ t b r = , , and r N = 2, … . The displacement field u i for the whole multilayered stack is then defined through an opportune assembly procedure of the layer-specific contributions u i k ( ) , which depends on the ESL or LW description.

In an ESL approach, the thickness functions are defined as Taylortype expansion and only one variable u iτ is used for the whole multilayer, i.e., the layer index k ( ) in Eq. ( 3) may be dropped off and the following thickness functions are used:

F F z F z r N = 1, = , = ( = 2, … -1) b t N r r (4)
The ESL description can be enhanced by including the Zig-Zag function F ZZ (z) proposed by Murakami in order to allow slope discontinuities at layers' interfaces [START_REF] Murakami | Laminated composite plate theory with improved in-plane response[END_REF]. In this case, the Zig-Zag function replaces the highest expansion order and the following functions are used:

F F F z F z r N = 1, = ( ), = ( = 2, … -1) b t Z Z r r ( 5
)
where Murakami's ZigZag Function (MZZF) is defined as

⎛ ⎝ ⎜ ⎞ ⎠ ⎟ F z ζ z ζ z z z z z z ( ) = (-1) ( ) with ( ) = 2 - - + 2 ZZ k k k t k b k t k b k ( ) ( ) ( ) ( ) (6)
Note that F ZZ (z) is expressed in terms of the non-dimensional layer-

specific coordinate ζ -1 ≤ ≤ +1 k
and it provides a linear piecewise function of bi-unit amplitude across the thickness of each layer k. More details about the use of MZZF in variable kinematics PVD-based models can be found in the paper by Demasi [START_REF] Demasi | Partially Zig-Zag advanced higher order shear deformation theories based on the Generalized Unified Formulation[END_REF].

The assumptions for a LW description are formulated in each layer k as in Eq. ( 3), where the thickness functions are defined by linear combinations of Legendre polynomials P ζ ( ) r k as follows:

F ζ P ζ P ζ F ζ P ζ P ζ ( ) = ( ) + ( ) 2 ; ( )= ( ) -( ) 2 ; t k k k b k k k 0 1 0 1 (7)
F ζ P ζ P ζ r N ( ) = ( ) - ( ) ( = 2, … ) r k r k r k -2
where ζ k is the non-dimensional coordinate introduced in Eq. ( 6). The Legendre polynomials of degree 0 and

1 are P ζ ( ) = 1 k 0 and P ζ ζ ( ) = k k 1
, respectively; higher-order polynomials are defined according to the following recursive formula:

P ζ n ζ P ζ nP ζ n ( ) = (2 + 1) ( ) - ( ) + 1 n k k n k n k +1 -1 (8)
which leads to the following expressions for the polynomials employed if N=4:

P ζ ζ P ζ ζ ζ P ζ ζ ζ ( ) = 3 -1 2 ; ( )= 5 -3 2 ; ( )= 35 8 - 15 4 + 3 8 k k k k k k k k 2 2 3 3 4 4 2 (9)
It is finally emphasized that the chosen thickness functions for a LW model satisfy the following properties are the corresponding linear Lagrange interpolation functions.

ζ F F F = 1: = 1, = 0, = 0 k t b r (10a) ζ F F F = -1: = 0, = 1, = 0 k t b r (10b) Therefore, u i k ( )

The stress and strain fields

The contributions to the strain and stress fields in each layer k are identified with respect to the bending (b), transverse normal (n) and transverse shear (s) deformation of the plate:

ϵ ϵ ϵ σ σ σ σ σ σ σ σ σ = [ϵ ϵ ϵ ]; = ϵ ; = [ϵ ϵ ] = [ ]; = ; = [ ] b k k k k n k k s k k k b k k k k n k k s k k k ( ) 1
( ) 2
( ) 6
( ) ( ) 3
( ) ( ) 5
( ) 4
( ) ( ) 1
( ) 2
( ) 6
( ) ( ) 3
( ) ( ) 5
( ) 4
( ) (11)
where Voigt's contracted notation for the symmetric strain and stress tensors has been invoked. The generic layer k is assumed to have a monoclinic material symmetry in the plate's reference frame x z ( , ) α ; thus the linear elastic constitutive law reads in matrix form as follows

⎡ ⎣ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ σ σ σ C C 0 C C 0 0 0 C ϵ ϵ ϵ = ∼ ∼ ∼ ∼ ∼ b k n k s k bb k bn k bn k nn k ss k b k n k s k ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) T ( 12
) u i u u i j i j ϵ = for = 1, 2, 3; ϵ = + for , = 1, 2, 3 and ≠ i i i j i j , (9--) , , i j i ( 13
)
Recalling the separation of the in-plane variables from the thickness direction, employed for expressing the assumed displacement field as in Eq. ( 3), the bending, transverse normal and transverse shear components of the strains defined by Eq. ( 13) are recast in the following matrix notation

ϵ Fb V x z z x ( , ) = ( ) ( ); b α k τ k τ k α ( ) (14a) ϵ Fn V x z z x ( , ) = ( ) ( ); n α k τ k τ k α ( ) (14b) ϵ Fs V x z z x ( , ) = ( ) ( ) s α k τ k τ k α ( ) (14c) where V x ( ) τ k α ( )
is the generalized strain vector of each layer defined as

V x u u u u u u u u u ( ) = [ ⋮ ⋮ ] τ k α k k k k k k k k k T ( ) 1 ( ) 1 , ( ) 1 , ( ) 2 ( ) 2 , ( ) 2 , ( ) 3 ( ) 3 , ( ) 3 ,
( )

τ τ τ τ τ τ τ τ τ 1 2 1 2 1 2 (15)
The explicit expressions for the matrices Fb Fn , τ τ and Fs τ can be found in Eq. (A-1).

In order to introduce the transverse shear locking correction proposed in the next section, the transverse shear strain field given in Eq. ( 14c) is split into the classical z -constant contribution γ 0 of standard FSDT, and a contribution γ h that depends on the thickness coordinate z and is related to high-order terms:

ϵ γ γ F s V F s V x z x x z x z x ( , ) = ( ) + ( , ) = ( ) + ( ) ( ) h s α k α h α k τ τ k α τ k τ k α 0 0 ( ) ( ) (16)
Note that Fs τ 0 is a matrix containing only constant values for all zdependency is contained in the matrix Fs Fs Fs z z ( ) = ( ) -

h τ k τ k τ 0 .
The matrices Fs τ 0 for ESL and LW models of order N are reported explicitly in Eq. (A-2) and Eq. (A-3), respectively.

The two-dimensional weak form of the boundary value problem

The displacement field defined in Eq. ( 2) is inserted into the PVD Eq. ( 1) yielding the following expression

∫ ∫ F δu σ F u x F δu t x - ϵ ( ) ( ) d + ( ) d = 0 V ij τ l ij s m i Γ τ i i α τ s σ τ ( 17
)
where the summation indexes τ and s are used for the expansions of the virtual and the actual displacements, respectively, and

F F z = ( ) τ τ
are the thickness functions evaluated at the z -coordinate at which the external tractions t are applied. Taking into account the separation of the thickness direction from the in-plane variables, recalling the strain definitions in Eqs. ( 14) and expressing the layer-specific stresses through the constitutive law Eq. ( 12), the virtual internal work can be expressed in the following matrix form

⎧ ⎨ ⎪ ⎪ ⎩ ⎪ ⎪ ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎫ ⎬ ⎪ ⎪ ⎭ ⎪ ⎪ ∑ ∫ ∫ ∫ V Fb Fn Fs C C 0 C C 0 0 0 C Fb Fn Fs V F δu σ F u x z x ϵ ( ) ( ) d = d d ∼ ∼ ∼ ∼ ∼ V ij τ l ij s m i Ω k N τ k e τ τ τ T bb k bn k bn k nn k ss k s s s k s k α =1 ( ) ( ) ( ) ( ) ( ) ( ) ( ) τ s L T k T ( ) (18)
The layer integrals along the thickness z k are explicitly carried out according to the following notation

∫ Z Z Z Z Z Fb C Fb Fb C Fn Fn C Fb Fn C Fn Fs C Fs z { , , , , } = { , , , , }d ∼
∼ ∼ ∼ ∼ ∼ ∼ ∼ ∼ ∼ bb k nb k bn k nn k ss k e τ T bb k s τ T bn k s τ T nb k s τ T nn k s τ T ss k s k ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) τs τs τs τs τs k ( ) (19)
Furthermore, the special treatment of the transverse shear strain given in Eq. ( 16) is introduced, which yields to the following subdivision

Z Z Z Z Z = + + + ∼ ∼ ∼ ∼ ∼ ss k ss k ss h k ss h k ss hh k ( ) 00( ) 0 ( ) 0( ) ( )
τs τs τs τs τs [START_REF] Reddy | An evaluation of equivalent-single-layer and layerwise theories of composite laminates[END_REF] This permits to identify those terms that involve the z -constant transverse shear strains responsible of the locking phenomenon. The two-dimensional matrix expression for the virtual internal work is finally written as [START_REF] Biscani | Variable kinematic plate elements coupled via Arlequin method[END_REF] The terms in Eq. ( 17) related to external tractions are developed in a similar manner; details can be found elsewhere [START_REF] Carrera | Classical and advanced multilayered plate elements based upon PVD and RMVT. Part 1: derivation of finite element matrices. Part 2: numerical implementations[END_REF][START_REF] Carrera | Finite Element Analysis of Structures through Unified Formulation[END_REF] and are here omitted for brevity.

⎪ ⎪ ⎧ ⎨ ⎩ ∫ ∫ ∑ V Z Z Z Z Z Z Z Z V F δu σ F u x x ϵ ( ) ( ) d = [ + + + + + + + ] }d ∼ ∼ ∼ ∼ ∼ ∼ ∼ ∼ V ij τ l ij s m i Ω k N τ k bb k bn k nb k nn k ss hh k ss h k ss h k ss k s k α =1 ( ) ( ) ( ) ( ) ( ) ( ) 0 ( ) 0

Finite element approximations

The solution of the two-dimensional boundary value problem expressed in the previous section is sought in weak form according to the finite element (FE) method. This section presents the four-node quadrilateral FE approximation and details out the new QC4 approximation for avoiding transverse shear locking problems and minimizing the convergence rate loss for distorted meshes.

Approximation for the geometry

The physical domain x x ( , ) 1 2 is mapped onto the reference, bi-unit square element, defined by the natural non-dimensional coordinates ξ η ( , ), through to the bilinear Lagrange interpolation functions N ξ η ( , ):

x N ξ η X i = ( , ) ( = 1, …4) α i α i (22)
where X αi are the xα coordinates of the four corner nodes, see Fig. 2.

Isoparametric interpolation

According to the isoparametric approach, the same interpolation defined by Eq. ( 22) is employed for approximating the unknown functions and their in-plane derivatives contained in V x ( )

τ k α ( ) : V B q i = ( = 1, …4) τ k τ i τ i k ( ) ( ) (23)
where

q u u u = [ ] τ i k k k k i T ( ) 1
( ) 2
( ) 3
( ) τ τ τ
is the DOF vector of the node i related to the layer k and the expansion order index τ. The 9×3 matrix B τ i Fig. 2. Four-node quadrangular element in the physical frame x x , 1 2 and natural frame ξ η , .

where superscript T is the transposition operator. The constitutive law is obviously defined for each layer k for it depends on the layer's orthotropic elastic properties C P k Q () (PQ,∈{1,6)}and on the orientation angle θ (k) . Within the small strain and small displacement approximation, the linearized strain components are defined as

∫ K B Z Z Z Z Z B x = [ + + + + ] d ∼ ∼ ∼ ∼ ∼ ij τs k Ω τ i T bb k nb k bn k nn k ss k s j α ( ) ( ) ( ) ( ) ( ) ( )
τs τs τs τs τs [START_REF] Carrera | Classical and advanced multilayered plate elements based upon PVD and RMVT. Part 1: derivation of finite element matrices. Part 2: numerical implementations[END_REF] Explicit expressions of this model-invariant kernel can be found in several works [START_REF] Carrera | Classical and advanced multilayered plate elements based upon PVD and RMVT. Part 1: derivation of finite element matrices. Part 2: numerical implementations[END_REF][START_REF] Carrera | Finite Element Analysis of Structures through Unified Formulation[END_REF]. These references report also the expansion procedure over the summation indexes τ and s as well as the assembly of all layers k, which are here omitted for brevity.

For the present four-node element, a 2×2 Gauss quadrature scheme is employed for an exact numerical evaluation of the domain integral (the 4 Gauss points are represented by the black triangles in Fig. 2). In order to contrast the shear locking phenomenon, the transverse shear stiffness contribution related to Z ∼ ss k ( ) τs has been reduced-integrated [START_REF] Carrera | Classical and advanced multilayered plate elements based upon PVD and RMVT. Part 1: derivation of finite element matrices. Part 2: numerical implementations[END_REF] or treated within a MITC approach [START_REF] Carrera | MITC technique extended to variable kinematic multilayered plate elements[END_REF]. In this work, however, a new interpolation for the z -constant transverse shear strain is presented, which affects only the contributions Z Z , ∼ ∼ ss ss h 00 0 and Z ∼ ss h 0 .

The new QC4 interpolation

The isoparametric interpolation for the transverse shear field leads to a locking phenomenon because of the incompatibility of the polynomial spaces defined by the sum of u ατ and the in-plane derivative u 3 , τ α (α = 1, 2 for γ 13 0 and γ 23 0 , respectively) [START_REF] Polit | Robust C 0 high-order plate finite element for thin to very thick structures: mechanical and thermo-mechanical analysis[END_REF]. The locking pathology is associated to the z -constant part only, for higher-order contributions depend on the plate thickness and vanish naturally in the thin plate limit. Thus a new, field-compatible interpolation for the four-node element is constructed for the z -constant part γ 0 : this is an extension to arbitrary variable kinematics plate models of the approach proposed in [START_REF] Polit | A new eight-node quadrilateral shear-bending plate finite element[END_REF] for FSDT plate elements and in [START_REF] Polit | Robust C 0 high-order plate finite element for thin to very thick structures: mechanical and thermo-mechanical analysis[END_REF] to a refined plate element. The new interpolation is constructed as follows:

• In order to enhance the element's robustness for distorted shapes, the zconstant part of transverse shear strain components is first written in the element's natural coordinate system ξ η ( , ) as:

⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ Fs U γ ξ η γ ξ η ξ η ( , ) ( , ) = (, ) ξ η k τ τ k 0 0 ( ) 0 ( ) (25)
where

⎡ ⎣ ⎢ ⎤ ⎦ ⎥ U ξ η u u u u u u u u u ( , ) = ⋮ ⋮ τ k ξ k ξ k ξ k η k η k η k k k k T ( ) ( ) , ( ) , ( ) ( ) , ( ) , ( ) 3 ( ) 3 , ( ) 3 , ( )
τ τ ξ τ η τ τ ξ τ η τ τ ξ τ η ( 26
)
is the projection onto the reduced natural coordinates ξ η , of the generalized strain vector V τ k ( ) .

• The field-compatible approximation shall be constructed upon enhancing the polynomial space of the transverse deflection u k 3 ( ) τ so that its derivative matches the bi-linear space of the in-plane displacements u ξτ and u η τ . For this, a quadratic serendipity inter- polation is assumed for u k

A 3 , τ ξ , u ( ) B 3 , τ η , u ( ) C 3 , τ ξ and u ( ) D 3 , τ η
, at the mid-side points, see Fig. 2. The supplementary DOFs are subsequently expressed in terms of the DOFs at the four corner nodes by imposing the tangential component of the transverse shear strain to be constant along the element's edges, viz.

γ ξ η γ ξ η ( , = ± 1) = constant; ( = ± 1, ) = constant ξ η 0 0 (27)
A new interpolation for the transverse deflection u 3τ is thus obtained through these four equations. The resulting four-node FE shall be denoted QC4 due to the quadratic and constant approximations employed for the transverse deflection and the tangential transverse shear strain, respectively.

• The new field-compatible interpolation for γ γ , ξ η 0 0 is defined by the polynomial basis obtained from the intersection sets of monomial terms in ξ and η:

γ u u η γ u u ξ ( ) = ( ) + ( ) = {1, } ( ) = ( ) + ( ) = {1, } ξ ξ η η 0 3 , 0 3 , τ τ ξ τ τ η (28)
• According to the polynomial basis spanning the transverse shear strain components, only two points are required for defining the resulting linear variations along the reduced natural coordinates.

Choosing these points to be the four mid-side points A D -(see Fig. 2), one has

γ ξ η C ξ η γ ξ η C ξ η γ ξ η γ ξ η C ξ η γ ξ η C ξ η γ ξ η ( , ) = ( , ) ( , ) + ( , ) ( , ) ( , ) = ( , ) ( , ) + ( , ) ( , ) ξ ξ ξ A A ξ ξ C C η η η D D η η D D 0 1 0 2 0 0 1 0 2 0 (29a)
with the following interpolating functions:

Cξ ξ η η Cξ ξ η η Cη ξ η ξ Cη ξ η ξ ( , ) = (1 -) 2 ; ( , )= (1 + ) 2 ( , ) = (1 -) 2 ; ( , )= (1 + ) 2 1 2 1 2 (29b)
• The physical transverse shear strains are finally deduced from the transverse shear strains in the reference domain as

⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ J γ ξ η γ ξ η γ ξ η γ ξ η ( , ) ( , ) = ( , ) ( , )
where the J is the 2×2 Jacobian matrix evaluated at the 4 Gauss points used for integrating the stiffness matrix.

The expression for the QC4 approximation for the z -constant transverse shear strain field can be finally written as

γ J Fs B q = τ i τ i k 0 -1 0 () (31)
where B i is the matrix containing the modified QC4 interpolation functions and their derivatives with respect to the physical coordinates x α . The explicit expression of the four B i matrices is given in Eq. (A-4).

Numerical results

Several numerical benchmark problems are considered for displaying the accuracy and robustness of the proposed QC4 CUF-based plate element. The classical CUF acronyms are used for naming the various plate models: the polynomial order N is appended to a string that identifies whether the model is ESL (EDN), Zig-Zag (EDZN) or LW (LDN).

A first set of problems involving a simple homogeneous and isotropic plate is used for assessing the element's performance, i.e., the properties of the stiffness matrix via an eigenvalue analysis, the convergence rate for thin and thick plates, and the sensitivity with respect to distorted element shapes; these latter tests are conducted on a square plate with distorted mesh, a skew plate and a circular plate. Previous analyses have shown that all CUF-based elements have the same sensitivity to shear locking irrespective of the order of expansion N [START_REF] Ottavio | Considerations on higher-order finite elements for multilayered plates based on a unified formulation[END_REF]. The case studies considered in this first problem set are hence computed with the ED2 model only, an ESL model that retains the three-dimensional constitutive law without Poisson locking.

A final set of problems is subsequently proposed for demonstrating the accuracy of the predicted displacements and stresses for homogeneous and composite, thin and thick plates. In this latter case study, several CUF-based elements are used for highlighting the variable kinematics framework.

The influence of the boundary conditions and type of loading are contains the shape functions i Nξ(,η) of Eq. ( 22) and their derivatives with respect to the physical coordinates x α . Introducing the FE approximation Eq. ( 23) into Eq. ( 21) yields the definition of the 3×3 fundamental nucleus of the classical CUF-based plate elements also addressed; Table 1 lists the acronyms used for denoting the various configurations. For all tests except the skew plate, symmetry is exploited to reduce the model to only a quarter plate.

Present QC4 results are compared against solutions obtained with the following isoparametric approaches: ISO full-integrated isoparametric element(2×2 Gauss points for the whole stiffness matrix) ISO-SI isoparametric element with selective integration(1 Gauss point for transverse shear stiffness contributions)

All elements are implemented as user subroutines into the commercial ABAQUS software. A dedicated pre-processing tool allows to prepare the FE model within the ABAQUS/CAE graphical interface.

The properties of the stiffness matrix

The eigenvalues of the stiffness matrix of a square element are analyzed for a thin and a thick plate, according to the following configuration geometry square element a a × (a=1), thickness e = 10 n -(n = 1, 3) boundary conditions (SA) material properties isotropic with E=10.92 and ν = 0.3

The results of the eigenvalue analysis are reported for thick and thin plates in Figs. 3 and4, respectively. The proposed graphics allow to recognize at a glance the eigenvalues associated to rigid-body modes as well as the gap between them and those associated to deformation modes. It is apparent that the proposed QC4 element and the fullintegrated isoparametric element have the correct number of rigidbody modes, i.e., 6, regardless of the element's slenderness S=a/e. On the contrary, the selective reduced integration scheme entails 3 spurious zero-energy modes, which indicates the possibility of an unstable behavior of both, thick and thin elements.

The transverse shear locking

A numerical test is carried out to assess the sensitivity of the proposed QC4 element to the transverse shear locking. The test is described as follows:

geometry square plate a a × (a=1), thickness e = 10 n -with n ∈ {0, 4} boundary conditions (SA) or (CL) on all sides loading (P) or (C) material properties isotropic with E=10.92 and ν = 0.3 mesh regular with N = 2, 4, 8, 16, 32 (see Fig. Four configurations are considered, defined by the combination of boundary conditions and loading. The concentrated force is applied only to thin plates with S ≥ 10 2 in order to avoid the singularity at the loaded node; thick plates (S=10) are considered only for a uniform pressure loading. In view of an assessment of the sensitivity to the transverse shear locking, the results are summarized in two different ways:

• Constant mesh, varying slenderness ratio Three regular meshes with N = 4, 8, 16 have been considered for the quarter plate (see Fig. 5) and, for each mesh, the results are given for varying length-to-thickness ratios S=a/e: Figs. 6 and7 report the results for plates subjected to a uniform pressure load and a concentrated force, respectively. For both loading types, the fullintegrated isoparametric element is shown to suffer a very strong locking as the plate becomes thin (S ≥ 10 2 ). On the contrary, both the QC4 and the ISO-SI elements are free from transverse shear locking and provide the reference Kirchhoff-Love values with good accuracy: a N=8 mesh is sufficient to obtain accurate results for all configurations.

• Convergence analysis for a thin or thick plate

The convergence curves of the transverse displacement at the plate's center with respect to the mesh density are shown in log-log scale in Fig. 8 and9 for various length-to-thickness ratios (S = 10, 10 , 10 2 4 ). The previous comments are confirmed: the convergence rate of the full-integrated ISO element is shown to strongly degrade as the length-to-thickness ratio S increases, a behavior that clearly indicates the presence of the transverse shear locking pathology; on the contrary, the convergence rate of the QC4 and ISO-SI elements remains sensibly the same irrespective of the length-to-thickness ratio. The monotone convergence of the QC4 element can be appreciated. As expected, the accuracy of the results depends on the test case for a given mesh.

The distortion tests

In this section, the sensitivity of the present FE to the mesh distortion is illustrated on three test cases widely used in open literature, namely the square, skew and circular plate.

The square plate test

This standard test is classically used in order to investigate the mesh sensitivity in plate bending problems. The data are given as follows:

geometry square plate a a × with a = 100 and thickness e = 0.1 boundary conditions (SA) or (CL) on all sides loading (P) or (C) material properties isotropic with E = 10.92 and ν = 0.3 mesh N = 2 for the quarter plate, with the distortion parameter s ∈ {-12, -8, -4, 0, 4, 8, 12} (see Fig. 10) results transverse displacement U u a a = ( /2, /2, 0) 3 3 at plate's center reference value transverse displacement U 3 (0) for the regular mesh (s=0)

In this test, the distorted meshes are characterized by the parameter s defining the coordinates of the mid-node of the quarter plate, which is located in the undistorted mesh (s=0) at

X X a = = /4 1 2
. The parameter s may be positive or negative, as illustrated in Fig. 10, and defines the coordinates of the mid-node as a s a s ( /4 + , /4 + ): for the most distorted meshes (s = ± 12) the mid-node is hence located at ( ± 37, ± 37). Note that it is not usual in open literature to consider positive and negative values for the parameter s.

The transverse displacement U 3 at the center node is normalized with respect to the value U 3 0 obtained with the regular, undistorted mesh (s=0). Fig. 11 shows the results obtained by the QC4, ISO and ISO-SI elements with ED2 kinematics for the four configurations defined by combining the two boundary and loading conditions. As expected, the largest errors occur for the most distorted meshes. The maximum absolute errors occurring for the four configurations are summarized in Table 2. The results indicate that QC4 is the most robust element with respect to mesh distortion for all boundary conditions, with a maximum error that does not exceed 20%. The highest distortion sensitivity is displayed by the full-integrated ISO element.

Razzaque's skew plate test

This benchmark originally proposed by Razzaque [START_REF] Razzaque | Program for triangular bending elements with derivative smoothing[END_REF] is illustrated in Fig. 12 and defined by the following data: geometry skewed square plate a a × with a = 100 and thickness e = 0.1, skew angle ψ = 60°b oundary conditions (SA) on straight edges, oblique edges are free loading (P) material properties isotropic with E = 10 3 and ν = 0.31 mesh N = 2, 4, 8, 16, 32, 64 skew elements for whole plate results non-dimensional transverse deflection at plate's center:

U u D q a = (0, 0, 0) 10 Table 3 reports the non-dimensional transverse displacement at the center of the skew plate, obtained by the QC4, ISO-SI and ISO elements with ED2 kinematics for different mesh densities N. The table also reports the percentage errors with respect to Razzaque's reference solution, which allows to identify the convergence rate of the three elements: the ISO-SI and QC4 elements show the same convergence rate and recover the reference solution with approximately 1% error with a N=16 mesh. The ISO element yields very poor results as a consequence of two causes, namely the severe transverse shear locking affecting the response of the thin plate (S a e = / = 10 3 ) and the high sensitivity to non-rectangular element shapes. This test thus confirms the robustness of the present QC4 element.

The circular plate test

A last mesh distortion sensitivity test is considered, which concerns the bending of a clamped circular isotropic plate subjected to a uniform pressure load according to the following data: geometry circular plate of radius R = 5 and thickness e = 0.1 boundary conditions (CL) on the external perimeter loading (P) material properties isotropic with E = 1.7472 10 7 and ν = 0.3 mesh N = 3, 12, 28, 60 elements for a quarter plate (see Fig. 13) results

deflection U u = (0, 0, 0) 3 3 at the center of the plate reference value Kirchhoff-Love theory U = 0.61147 10

3 ref -6 [53]
The evolution of the ratio between the central deflection and the Kirchhoff-Love solution with respect to the mesh density is reported in Fig. 14. The convergence is compared for the QC4, ISO and ISO-SI elements with ED2 kinematics. As previously discussed, the ISO element suffers the distorted shapes and for the present case it shows errors exceeding 90% even with the more refined mesh. It can be noticed that the QC4 element has a good convergence rate and its accuracy is very satisfactory. In particular, the error becomes less than 3% for N ≥ 30, whereas the converged result of the ISO-SI FE is still affected with an error of about 10%.

Mechanical analysis of isotropic plates

The focus is here set on the convergence of displacement and stresses predicted by the proposed QC4 element for a plate bending problem. The benchmark problem is defined as follows:

geometry square plate a a × and length-to-thickness ratio

) plates under a pressure load. The convergence behavior for displacements and stresses of the QC4 element is reported in Fig. 15, where both very thick (S=5) and thin (S=100) plates are considered. The present results are obtained with an ED2 kinematics and the FE error is defined with respect to the closed-form Navier-type solutions obtained with the same ED2 model. The convergence rate is shown to be independent of the length-to-thickness ratio and is very good for both, displacements and stresses: note that all discretization errors are less than 1% if a mesh with N ≥ 8 elements is used. Fig. 16 displays the distributions along the z -axis of the normalized transverse stresses σ z σ a z q S ( ) = (0, /2, )/( ) [START_REF] Lo | A high-order theory of plate deformation -Part 1: homogeneous plates -Part 2: laminated plates[END_REF] 13 0 and σ z σ a a z q ( ) = ( /2, /2, )/ 33 33 0 , computed directly from the constitutive equations, when the plate thickness is subdivided into N = 1, 3, 5 L numerical layers with an LD4 kinematics. A FE mesh with N=8 has been used for the results in Fig. 16. It is shown that, for the present case of isotropic plate, the refinement of the response across the thickness by means of numerical layers has no effect and the converged solution is reached with N L =1. For a thick plate (S=5), the LD4 model is capable to accurately recover the natural boundary conditions at the outer surfaces of the plate, see Fig. 16 (left). In case of a thin plate (S=100), however, the transverse normal stress distribution is highly inaccurate, see Fig. 16 (right). Fig. 17 shows that, for a thin plate, it is necessary to refine the in-plane FE mesh in order to enhance the accuracy of the σ 33 distribution across the thickness, up to the fulfillment of the natural boundary conditions at the outer surfaces of the plate, see also [START_REF] Polit | Robust C 0 high-order plate finite element for thin to very thick structures: mechanical and thermo-mechanical analysis[END_REF]. reference value 3D exact elasticity results given by Pagano [START_REF] Pagano | Exact solutions for rectangular bidirectional composites and sandwich plates[END_REF] Three ESL models (ED2, ED4 and EDZ4) and two LW models (LD2 and LD4) are considered in order to assess different plate theories. Present QC4 results, identified by the suffix -FEM, are compared against the exact solution proposed by Pagano [START_REF] Pagano | Exact solutions for rectangular bidirectional composites and sandwich plates[END_REF] and the Navier closed-form solutions obtained with the corresponding plate model. These latter results are identified by the suffix -anlt.

Tables 4 and5 list the non-dimensional displacement and stress values, respectively, for all considered models, along with the percentage error with respect to the reference solution. Three different lengthto-thickness ratios have been considered (S = 4, 10, 100). While simple ESL models are sufficiently accurate for the global response of thin plates, high-order LW models are necessary for thick plates and for an accurate prediction of the local transverse stresses. The results show that the N=8 mesh provides converged FEM results with respect to the corresponding closed-form Navier solution, thus confirming the performances for the isotropic plate. Moreover, the same behavior is found for all considered plate models, which means that the QC4 approximation could be successfully extended to arbitrary kinematics including Zig-Zag and LW models.

It is also interesting to illustrate the through-thickness distribution of the response. Figs. 18, 19 and 20 display the distributions of nondimensional displacements, in-plane stresses and transverse stresses, respectively, obtained for the moderately thick plate S ( = 10) and with the high-order ED4 kinematics. All stresses are directly evaluated from the compatible strains by means of the constitutive relation. The distributions of the displacements (Fig. 18) and the in-plane stresses (Fig. 19) agree well with the reference elasticity solution. However, note that the Zig-Zag shape of the u 1 displacement is not captured by the considered ED4 model. Larger discrepancies are obviously obtained in Fig. 20 for the distributions of the transverse shear stresses because of their interlaminar discontinuity, which is typical of any PVD-based formulation. Thanks to the variable kinematics approach, the proposed FE is nonetheless particularly efficient and accurate in the computation of displacements and stresses for multilayered plates.

Conclusion

This paper has introduced a new locking-free four-node quadrilateral finite element (FE) for displacement-based variable kinematics plate models expressed through Carrera's Unified Formulation (CUF). The remedy to the shear locking pathology pertains to the B-bar method and consists in a new QC4 approximation that is constructed for the zconstant term of the transverse shear strain only. The robustness and accuracy of the proposed element has been thoroughly assessed by referring to a number of tests that are recommended whenever new FE applications are proposed: classical eigenvalue counts, convergence behavior for thin and thick plates under various boundary and loading conditions, as well as three case studies involving distorted meshes. The numerical results confirmed the superiority of the proposed FE in comparison to classical isoparametric approaches with full or reduced integrations, i.e., it has a correct rank, is free of transverse shear locking and is less sensitive to distorted element shapes, showing a high convergence rate for both displacements and stresses. The proposed variable kinematics FE, implemented as a user subroutine into Abaqus, provides a robust tool for the analysis of composite laminates, for which different models may be used for adapting the computational cost in case of thin or thick plates or whether a global or local response is required.

Further studies shall concern the optimization of the computational

Table 4

Simply-supported (0°/90°/0°) square plate under bi-sinusoidal pressure load: Nondimensional displacements for length-to-thickness ratios S = 4, 10, 100.

S Model ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ u - e 1 2 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ u - e 2 2 u (0) 3 4
Ref. [ Simply-supported (0°/90°/0°) square plate under bi-sinusoidal pressure load: Non-dimensional stresses for length-to-thickness ratios S = 4, 10, 100.

S Model

⎛ ⎝ ⎜ ⎞ ⎠ ⎟ σ e 11 2 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ σ - - e 22 6 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ σ - e 12 2
σ (0) cost for FE models of composite structures: on the one hand, the FE implementation of a more general variable kinematics modeling approach will be considered [START_REF] Ottavio | A Sublaminate Generalized Unified Formulation for the analysis of composite structures and its application to sandwich plates bending[END_REF]; the variable kinematics models shall, on the other hand, be effectively employed within a global-local approach that limits the use of expensive, highly accurate and quasi-3D models to small model portions.

Appendix A. Appendix

A.1. Matrices with thickness functions

The arrays defining the through-thickness behavior of the strain field in Eqs. ( 14 The thickness functions describing the z -constant transverse shear strain γ 0 in Eq. ( 16) have the following expressions depending on the model:

• ESL model: the following matrices are used for all expansion orders N ≥ 1:

⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎤ ⎦ ⎥

Fs

Fs

= 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 = 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 • LW model: with reference to the notation in Eq. ( 3), the approximation is defined as Fs Fs Fs Fs

Fig. 1 .

 1 Fig. 1. Coordinates and notation used for the description of the composite plate.

b

  are the physical displacement components at the top and bottom of the k th layer, respectively, and F ζ ( ) t k and F ζ ( ) b k

Fig. 4 .

 4 Fig. 4. Eigenvalues of the stiffness matrix for an ED2 model and thin plate (a e / = 10 3 ).

Fig. 5 .

 5 Fig. 5. Regular meshes for a quarter of the plate.

Fig. 6 .

 6 Fig.6. Transverse displacement with respect to slenderness ratio for thick (S=10) to very thin plates (S = 10 5 ) under a pressure load.

Fig. 7 .

 7 Fig. 7. Transverse displacement with respect to slenderness ratio for thin plates (S ≥ 10 2 ) under a concentrated load.

Fig. 8 .

 8 Fig. 8. Convergence of the transverse displacement for thick (S=10) and thin (S = 10 , 10

Fig. 9 .Fig. 10 .

 910 Fig. 9. Convergence of the transverse displacement for thin plates (S = 10 , 10 2 4 ) under a

  , 5, 10, 100 boundary conditions (SA) loading bi-sinusoidal pressure on the top surface: with E=73 GPa and ν = 0.34 results non-dimensional response defined as reference value closed-form Navier-type solution (ED2 model)

Fig. 11 .

 11 Fig. 11. Variation of the normalized central deflection with respect to the parameter s for square plates with various loading and boundary conditions (ED2 model).

4. 5 .

 5 Mechanical analysis of laminated composite platesThis final case study involves simply-supported cross-ply laminates subjected to a bi-sinusoidal pressure load, for which Pagano provided exact elasticity results[START_REF] Pagano | Exact solutions for rectangular bidirectional composites and sandwich plates[END_REF]. The following configuration is considered: geometry square plate a a × and length-to-thickness ratio S three-layered cross-ply plate (0°/90°/0°); each la- yer has thickness e e transverse directions, respectively) mesh N=8 for the quarter plate results displacements and stresses are made non-dimensional ac-

Fig. 13 .

 13 Fig. 13. FE meshes with N = 3, 12, 28 elements of a quarter of the circular plate.

Fig. 14 .

 14 Fig. 14. Convergence study of the central deflection for the circular plate.

Fig. 15 .

 15 Fig.[START_REF] Vidal | Assessment of the refined sinus plate finite element: free edge effect and Meyer-Piening sandwich test[END_REF]. Simply-supported square isotropic plate under bi-sinusoidal pressure load: convergence of displacements and stresses with respect to corresponding closed-form solutions (ED2 model) for length-to-thickness ratios S=5 (left) and S=100 (right).

Fig. 16 .Fig. 17 .

 1617 Fig.[START_REF] Polit | Robust C 0 high-order plate finite element for thin to very thick structures: mechanical and thermo-mechanical analysis[END_REF]. Influence of the number of numerical layers N L on the through-thickness distribution of σ 13 and σ 33 for a thick (S=5, left) and a thin (S=100, right) plate (isotropic plate loaded at top surface, LD4, N=8).

Fig. 18 .

 18 Fig.18. Through-thickness distribution of u 1 (left), u 2 (middle) and u 3 (right) for the (0°/90°/0°) square plate (S=10).

Fig. 19 .

 19 Fig.[START_REF] Wenzel | Coupling of heterogeneous kinematics and finite element approximations applied to composite beam structures[END_REF]. Through-thickness distribution of σ 11 (left), σ 22 (middle) and σ 12 (right) for the (0°/90°/0°) square plate (S=10).

Fig. 20 .

 20 Fig.[START_REF] Reddy | An evaluation of equivalent-single-layer and layerwise theories of composite laminates[END_REF]. Through-thickness distribution of σ 13 (left), σ 23 (middle) and σ 33 (right) for the (0°/90°/0°) square plate (S=10).

Table 1

 1 Acronyms for boundary conditions and type of loading.

	Boundary conditions	Loading
	(SA)	Simply supported	(P)	Uniform load q 0 at top surface
	(CL)	Clamped	(C)	Concentrated force P at plate center

Fig. 3. Eigenvalues of the stiffness matrix for an ED2 model and thick plate (a e / = 10).

Table 2

 2 Maximum absolute error of the central deflection and corresponding svalue for various boundary and loading conditions.

	Config.	QC4		ISO-SI		ISO
		error [%]	s	error [%]	s	error [%]	s
	SA-P	8.6	-12	14.2	-12	27.0	-12
	SA-C	8.9	+12	16.5	-12	17.9	-12
	CL-P	11.2	-12	19.7	-12	32.8	-12
	CL-C	18.9	+12	18.9	+12	30.8	+12

Fig. 12. Razzaque's skew plate test (mesh N=4).

Table 3 Non

 3

	S	N=2	N=4	N=8	N=16	N=32	N=64
	ISO	0.0000	0.0000	0.0002	0.0007	0.0029	0.0114
		(-100.0%) (-100.0%) (-100.0%) (-99.9%) (-99.6%) (-98.6%)
	ISO-SI 0.3269	0.6616	0.7547	0.7843	0.7908	0.7927
		(-58.9%)	(-16.7%)	(-5.0%)	(-1.3%)	(-0.5%)	(-0.2%)
	QC4	0.3727	0.6609	0.7533	0.7838	0.7908	0.7927
		(-53.1%)	(-16.8%)	(-5.2%)	(-1.3%)	(-0.5%)	(-0.2%)

-dimensional deflection at the center point of Razzaque's skew plate (values in parentheses are the percentage errors with respect to the reference solution).

Table 5

 5

A.2. QC4 interpolation

The non-zero terms of the 9×3 matrices B i ( = 1 -4)

(3) 3 (A-4)

where J αβ i ( ) is the Jacobian at node i.


Marieke Huisman
email: m.huisman@utwente.nl

Rosemary Monahan
email: rosemary.monahan@nuim.ie

Peter Müller
email: peter.mueller@inf.ethz.ch

Andrei Paskevich
email: andrei.paskevich@lri.fr

Gidon Ernst
email: gidon.ernst@unimelb.edu.au

Gidon Ernst Verifythis

VerifyThis 2018: A Program Verification Competition

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

Introduction

VerifyThis 2018 took place on April 14 and 15, 2018 in Thessaloniki, Greece, as a two-day verification competition in the European Joint Conferences on Theory and Practice of Software (ETAPS 2018). It was the sixth edition in the VerifyThis series after the competitions held at FoVeOOS 2011, FM 2012, Dagstuhl Seminar 14171 (April 2014), ETAPS 2015-2017.

The aims of the competition were:

to bring together those interested in formal verification, and to provide an engaging, hands-on, and fun opportunity for discussion; to evaluate the usability of logic-based program verification tools in a controlled experiment that could be easily repeated by others.

Typical challenges in the VerifyThis competitions are small but intricate algorithms given in pseudo-code with an informal specification in natural language. Participants have to formalise the requirements, implement a solution, and formally verify the implementation for adherence to the specification. There are no restrictions on the programming language and verification technology used. The time frame to solve each challenge is quite short (90 minutes) so that anyone can easily repeat the experiment. The verification challenges are available from the VerifyThis website http://www.pm.inf.ethz. ch/research/verifythis.html.

The correctness properties which the challenges present are typically expressive and focus on the inputoutput behaviour of programs. To tackle them to the full extent, some human guidance within a verification tool is usually required. At the same time, considering partial properties or simplified problems, if this suits the pragmatics of the tool, is encouraged. The competition welcomes participation of automatic tools as combining complementary strengths of different kinds of tools is a development that VerifyThis would like to advance.

Submissions are judged for correctness, completeness, and elegance. The focus includes the usability of the tools, their facilities for formalizing the properties and providing helpful output.

VerifyThis 2018 consisted of three increasingly difficult verification challenges, selected to showcase various aspects of software verification. Before the competition, an open call for challenge submissions was made. Overall, we received six proposals, of which we selected one for the competition, and some were excluded because published solutions are available. The challenges (presented later) provided reference implementations at different levels of detail.

Eleven teams participated (Table 1) in this edition of the competition. Teams of up to two people were allowed and physical presence on site was required. In addition, one non-competing solution was developed off site simultaneously. We particularly encouraged participation of:

student teams (including PhD students); non-developer teams using a third-party tool; several teams using the same tool.

Teams using different tools for different challenges (or even for the same challenge) were also welcome. Student participants could apply for travel grants, sponsored by Amazon Web Services, to ease their participation.

The website of the 2018 instalment of VerifyThis can be found at http://www.pm.inf.ethz.ch/research/ verifythis/Archive/2018.html. More background information on the competition format and the rationale behind it can be found in [START_REF] Huisman | On the organisation of program verification competitions[END_REF]. Reports from previous competitions of similar nature can be found in [START_REF] Klebanov | The 1st Verified Software Competition: Experience report[END_REF][START_REF] Bormer | The COST IC0701 verification competition 2011[END_REF][START_REF] Filliâtre | The 2nd Verified Software Competition: Experience report[END_REF][START_REF] Huisman | A program verification competition[END_REF][START_REF] Huisman | A program verification competition[END_REF] and in the special issue of the International Journal on Software Tools for Technology Transfer (STTT) on the VerifyThis competition 2012 (see [START_REF] Huisman | VerifyThis[END_REF] for the introduction).

Invited Tutorial

We started the competition day with an interactive invited tutorial by Alexander Summers on the Viper tool [START_REF] Müller | Viper: A verification infrastructure for permission-based reasoning[END_REF]. Viper is a verification tool for permission-based reasoning about sequential and concurrent programs. Its primary purpose is to establish an infrastructure that can be relied on by front ends for higher-level programming languages to discharge their verification tasks (an example is VerCors [START_REF] Blom | The VerCors tool for verification of concurrent programs[END_REF], which was used in the competition by teams 7 and 8). A distinguishing feature of the tool is its support for the "magic wand" connective and quantified permissions.

Additionally, Viper provides a user-friendly surface language and integration into the Visual Studio Code editor. It can therefore be used stand-alone and was used in this mode in the competition (by team 10). During the tutorial, the participants had the opportunity to learn about Viper's methodology and to try out Viper on example verification problems. 1

Post-mortem Sessions

Two concurrent post-mortem sessions were held the day after the competition, where participants explained their solutions and answered questions of the judges.

During one session, the judges asked the teams questions in order to better understand and appraise their solutions. These sessions provided an excellent opportunity to explore the strengths and weaknesses of the solutions presented by each team, while acquiring more detailed knowledge of the verification tools used. In parallel, all other participants presented their solutions, leading to lively discussion and exchange about tool developments. These presentations were also attended by some nonparticipants.

http://www.pm.inf.ethz.ch/research/ verifythis/Participation.html

Judging Criteria

Limiting the duration of each challenge assists the judging and comparison of each solution. However, this task is still quite subjective and hence, difficult. Discussion of the solution with the judges typically results in a ranking of solutions for each challenge.

Criteria that were used for judging were:

-Correctness: is the formalisation of the properties adequate and fully supported by proofs? Were any bugs found in the code? -Completeness: are all tasks solved, and are all required aspects covered? Are any assumptions made? Is termination verified? -Readability: can the submission be understood easily, possibly even without a demo? -Effort and time distribution: what is the relation between time expended on implementing the program vs. specifying properties vs. proving? -Automation: how much manual interaction is required, and for what aspects? Does the solution make use of information from libraries? -Novelty: does the submission apply novel techniques?

What special features of the tool are used?

First Challenge: Gap Buffer

A gap buffer is a data structure for the implementation of text editors, which can efficiently move the cursor, as well as add and delete characters. 2 A gap buffer represents the editor's content as a character array a of size n, which has a gap of unused entries a[l], . . . , a[r -1], with respect to two indices l ≤ r. The data it represents is composed as

a[0], . . . , a[l -1], a[r], ..., a[n -1]
. The current cursor position is at the left index l, and if we type a character, it is written to a[l] and l is increased. When the gap becomes empty, the array is enlarged and the data from r is shifted to the right.

The gap buffer provides four operations (shown in Fig. 1, given to the participants as part of the description). Procedures left() and right() move the cursor by one character; insert() places a character at the beginning of the gap a[l]; delete() removes the character at a[l] from the range of text.

Verification Tasks

The intended behavior of the buffer should be specified in terms of a contiguous representation of the editor content. We suggested to use strings, functional arrays, sequences, or lists. The task was to verify that the gap buffer implementation satisfies this specification, and that every access to the array is within bounds. We encouraged to verify insert() with a precondition l < r first and to ignore the call to grow(). An extended task was to lift this restriction and to specify and verify the procedure grow(). We explicitly asked the participants to do this modularly, not by referring to the implementation of grow() in the proof of insert(). We encouraged to make use of built-in functionality for copying array ranges (such as System.arraycopy() in Java) instead of the loops in grow().

The participants had 60 minutes to implement, specify, and verify the challenge. It was intended to be straightforward to solve, not requiring advanced features of verification tools or familiarity with specific theories. The challenge exercises tool support for arrays, abstraction capabilities in the contracts/specifications of procedures, and modular verification.

Comments on Solutions

We received 10 complete solutions for part 1, compared to 1 partial and 4 complete solutions for part 2.

We observed that not all teams specified a contiguous representation, i.e., the fact that there is a gap was reflected in the specification which is arguably not an elegant approach when it would come to reasoning about a client (i.e., an editor) of this interface. Users of Dafny, Isabelle, KIV, KeY, and VerCors could benefit from their tool support for abstract sequences. Among those, we find Dafny's syntax for array slices and implicit promotion of arrays to mathematical sequences the most elegant (submitted by Matthias Ulbrich after the competition),

0 <= l <= r <= a.Length && s == a[..l] + a[r..]
where s: seq<char> is the abstract representation.

The solutions by the KIV and Isabelle teams expressed the specification independently from the code and used refinement to prove the correspondence.

The Frama-C team expressed the specification of delete() as two disjoint behaviors, nicely corresponding to the case distinction made in the implementation. Another remarkable aspect of their solution is the use of the realloc() library function in grow(), which is not only idiomatic in C but also entirely avoids reasoning about the first loop.

Among the few errors made in the preconditions and representation invariant, a common issue was to mistakenly specify l < r in the invariant instead of l ≤ r, which precludes moving the cursor to the rightmost position.

Separation Logic, as used by the Coq/CFML team, is very suitable to reason modularly about the array modifications. Their model is a representation based on Huet's "zipper" [START_REF] Huet | The zipper[END_REF], which is well-suited to specify the contract of grow() but exposes the gap in the buffer. The authors of that solution considered another abstraction to lists but did not attempt it due to the limited time.

With the exception of CIVL, all tools support modular reasoning about procedures (pre-/postconditions). In the latter case, a manual encoding was done.

Second Challenge: Colored Tiles

This challenge is based on Project Euler problem #114.

Alice and Bob are decorating their kitchen, and they want to add a single row of fifty tiles on the edge of the kitchen counter. Tiles can be either red or black, and for aesthetic reasons, Alice and Bob insist that red tiles come by blocks of at least three consecutive tiles. Before starting, they wish to know how many ways there are of doing this. The algorithm is as follows: The participants had to verify that at the end, count[50] contains the correct number. Since the algorithm works by enumerating the valid colorings, a good solution is expected to provide a nice specification of a valid coloring and to prove that each coloring counted by the algorithm is valid, no coloring is counted twice, and no valid coloring is missed.

Comments on Solutions

While no on-site participant completed the challenge, a number of submissions provided a full specification and partially finished proof. The approach most used for specification consisted in describing explicitly the set of all valid colorings for a given number of tiles and linking the values computed by the algorithm to the cardinality of this set. In the Isabelle solution, quite elegantly, the set was defined directly by set comprehension. In other solutions, the set of colorings was computed in a recursive or iterative fashion, following the structure of the algorithm. In the Dafny submission, the sought value was computed by counting valid entries in the set of all possible colorings.

The suggested definitions of a valid coloring turned out to be quite diverse. Some submissions defined it in a pointwise fashion, by stating that each read tile must have two red adjacent tiles or two red tiles to its left or to its right. The others defined it recursively, by describing valid extensions of shorter valid colorings. Most submissions described a coloring as a sequence of colors (Booleans), while the Viper solution used the run-length representation.

Having access to rich mathematical library, as was the case for Isabelle, considerably simplified specification and proof. In particular, the second-order sum operator allowed for succinct and convincing specification.

Worthy of note, the CIVL submission produced an automated proof of unicity: for the three categories of valid colorings (starts with a black tile, starts with n ≥ 3 red tiles followed by a black one, contains only red tiles), no coloring belong to two categories at once, and in the second category, no coloring has two different positions of the first black tile.

Third Challenge: Array-Based Queue Lock

Array-Based Queuing Lock (ABQL) is a variation of the Ticket Lock algorithm with a bounded number of concurrent threads and improved scalability due to better cache behaviour.

We assume that there are N threads and we allocate a shared Boolean array pass[] of length N. We also allocate a shared integer value next. In practice, next is an unsigned bounded integer that wraps to 0 on overflow, and we assume that the maximal value of next is of the form kN -1. Finally, we assume at our disposal an atomic fetch _ and _ add instruction, such that fetch _ and _ add(next,1) increments the value of next by 1 and returns the original value of next.

The elements of pass[] are spinlocks, assigned individually to each thread in the waiting queue. Initially, each element of pass[] is set to false, except pass[0] which is set to true, allowing the first coming thread to acquire the lock. Variable next contains the number of the first available place in the waiting queue and is initialized to 0.

Here is an implementation of ABQL in pseudocode:

procedure procedure procedure abql _ init() procedure procedure procedure abql _ release(my _ ticket) for for for i = 1 to to to N - Each thread that acquires the lock must eventually release it by calling abql _ release(my _ ticket), where my _ ticket is the return value of the earlier call of abql _ acquire(). We assume that no thread tries to reacquire the lock while already holding it, neither it attempts to release the lock which it does not possess. Notice that the first assignment in abql _ release() can be moved at the end of abql _ acquire().

The verification tasks were as follows:

1. Safety of ABQL under the given assumptions. Specifically, participants had to prove that no two threads can hold the lock at any given time. 2. Fairness, that is, the threads acquire the lock in order of request. 3. Liveness under a fair scheduler: each thread requesting the lock will eventually acquire it.

The participants had the liberty to adapt the implementation and specification of the concurrent setting as best suited for the verification tool of their choice. In particular, solutions with a fixed value of N were acceptable. It was expected, however, that the general idea of the algorithm and the non-deterministic behaviour of the scheduler were to be preserved.

Comments on Solutions

We received two solutions specifying and proving safety of ABQL in the general case (VerCors and Why3) and one solution proving safety and fairness for a bounded number of threads and a fixed maximal value of next (CIVL). Why3 submission also provided an unproved specification for fairness using a ghost queue.

Two main approaches were used for this challenge: direct implementation in tools with support for concurrency (Vercors, CIVL, Frama-C) or encoding the whole N-thread environment as a transition system. In the first case, the lock-free nature of the algorithm created an additional difficulty; both Vercors submissions added extra mutual exclusion primitives to the provided implementation. In the second case, the participants had to identify the relevant sequence points in the code; the number of states per thread ranged from 3 (KeY, Isabelle) to 4 (Viper) to 7 (Why3).

We would like to single out the specification framework provided by the VerCors tool, which allowed for natural and unencumbered thread-modular annotations.

Results and Closing Remarks

Awarded Prizes

We awarded five prizes. The announcements were made during the conference banquet and each winning team received a certificate and a prize in form of an Amazon voucher.

-Overall Best Team: Peter Lammich and Simon Wimmer (Isabelle, TU Munich) -Student Gold Medal: Raphaël Rieu-Helft (Why3, Inria) -Student Team Silver Medal: Wytse Oortwijn and Mohsen Safari (VerCors, Twente) -Distinguished Tool Feature: Jörg Pfähler and Stefan Bodenmüller (Augsburg) for their program refinement method in KIV -Best Challenge Submitted: Jean-Christophe Filliâtre, Colored Tiles

Lessons Learned

The two major challenges in organizing a competition like this are coming up with good riddles and judging the submissions.

For the first task, we would recommend to start searching as early as possible. Not many people will answer the call for problems and not all suggestions will be suitable for a 90 minute challenge. It may well happen that the organizers will have to invent a problem or two by themselves. Recruiting a knowledgeable and altruistic colleague or community member (who would have to forgo participation in the competition that year) for brainstorming sessions may also be of great help.

The organizers should remember that with more than a few actively developed tools for program verification on the scene, quite a lot of nice specification and verification problems (especially ones that can be handled in a limited time) have already been solved and published somewhere: in scientific papers, in online tutorials or simply somewhere on the tool's web page. Having rejected two promising ideas this year, we would advise our successors to spend some time looking for any previous work on the proposed challenges.

As for the second task, it is a good idea to prepare a long list of checkboxes for each problem beforehand and revise it constantly during the post-mortem interviews. For each aspect of each verification task, there will be some typical and some unique mistakes or omissions. Each of them may be only partially proved or not proved at all. Having a common evaluation grid, being ready to extend it, and to update it retroactively for the previously discussed submissions will greatly simplify the final systematization of the results. This evaluation grid must also provide space for qualitative observations: highlights or curious aspects of a given solution or tool.

Fig. 1 .

 1 Fig. 1. Gap Buffer implementation in pseudo-code.

  var var var count[51] // count[i] is the number of valid rows of size i count[0] := 1 // [] count[1] := 1 // [B] -cannot have a single red tile count[2] := 1 // [BB] -cannot have one or two red tiles count[3] := 2 // [BBB] or [RRR] for for for n = 4 to to to 50 do do do count[n] := count[n-1] // either the row starts with a black tile for for for k = 3 to to to n-1 do do do // or it starts with a block of k red tiles count[n] := count[n] + count[n-k-1] // followed by a black one end end end-for for for count[n] := count[n]+1 // or the entire row is red end end end-for for for

  1 do do do pass[my _ ticket] := false pass[i] := false pass[(my _ ticket + 1) mod mod mod N] := true end end end-for for for end end end-procedure procedure procedure pass[0] := true next := 0 end end end-procedure procedure procedure function function function abql _ acquire() var var var my _ ticket := fetch _ and _ add(next,1) mod mod mod N while while while not not not pass[my _ ticket] do do do end end end-while while while return return return my _ ticket end end end-function function function

Table 1 .

 1 Teams participating in VerifyThis 2018 (alphabetically by tool).

	# Team members	Tool	Team attributes
	1 Stephen Siegel	CIVL	[19] developer
	2 Armaël Guéneau, Cyprien Mangin	Coq+CFML	[4] student, developer
	3 Lionel Blatter, Jean-Christophe Léchenet Frama-C	[14] student
	4 Peter Lammich, Simon Wimmer	Isabelle Refinement Framework	[16] developer
	5 Stefan Bodenmüller, Jörg Pfähler	KIV	[5] student, developer
	6 Matthias Ulbrich	KeY, Dafny	[1, 17] developer (of KeY)
	7 Wytse Oortwijn, Mohsen Safari	VerCors	[2] student, developer
	8 Sebastian Joosten, Marieke Huisman	VerCors	[2] developer
	9 Jafar Hamin	Verifast	[13] student
	10 Marco Eilers, Alexander Summers	Viper	[18] developer
	11 Raphael Rieu-Helft	Why3	[6] student, developer
	Martin Clochard	Why3	[6] student, developer, off site
		procedure procedure procedure insert(x: char char char)	procedure procedure procedure grow()
	procedure procedure procedure left()	if if if l == r then then then
	if if if l != 0 then then then	// see extended task
	l := l -1	grow()
	r := r -1	end end end-if if if
	a[r] := a[l]	a[l] := x
	end end end-if if if	l := l + 1
	end end end-procedure procedure procedure
		procedure procedure procedure delete()
	procedure procedure procedure right()	if if if l != 0 then then then
	// similar to left()	l := l -1
	// but pay attention to the	end end end-if if if
	// order of statements	end end end-procedure procedure procedure
	end end end-procedure procedure procedure

var var var b := new new new char char char[a.length + K] // b[0..l] := a[0..l] for for for

Available from

See e.g. http://scienceblogs.com/goodmath/2009/02/18/ gap-buffers-or-why-bother-with-1

Acknowledgments

We thank Amazon Web Service for generous funding of travel grants and prizes. We thank Tej Chajed, Stefan Ciobaca, Ernie Cohen, Tim Denvir, Jean-Christophe Filliâtre, and Dennis Furey who submitted challenge proposals for the competition. The friendly easyconferences staff helped us with the organization of the competition dinner. Finally, we thank all participants for sitting through the-metaphorically and literally-hot sessions and for sharing and discussing their solutions.