Mariem Sbayou

Youssef Bouanan

Grégory Zacharewicz

Bruno Vallespir

BPMN Coordination And DEVS Network Architecture For Healthcare Organizations

Keywords: Business Process Management, Business Process, Healthcare Providers, Multi-Agent, Patient Pathways, Performance, Resources, Social Networks

HAL is

INTRODUCTION

Healthcare processes or pathways are typically described informally by text description or by semi-formal languages [START_REF] Eshuis | State of the Art Clinical Pathway Definition: Gap Analysis[END_REF]. A process is generally defined as a sequence of events that uses inputs to produce outputs. Besides, a business process is an activity or set of activities that will accomplish a specific organizational goal. Also, it is considered as a sequence of performed steps that drive information to produce goods and/or provides services. The business process needs to be managed and controlled; it attests the need to use a business process management methodology. According to [START_REF] Van Der Aalst | Business process management: A survey[END_REF] business process management (BPM) includes methods, techniques, and tools to support the design, enactment, management, and analysis of operational business processes.

It can be considered as a formalization of classical Workflow Management (WFM) systems and approaches. BPM solutions enable enterprise to measure and standardize processes and also provide reusable processes that can be networked. In this context, several languages are used for BPM; the most used one nowadays is Business Process Modeling Notation (BPMN). BPMN has been proposed by the Business Process Modelling Initiative (BPMI) and is currently maintained by the Object Management Group (OMG 2003) that provides this standard for Information Technology (IT) and business actors, it is an increasingly significant standard for process modeling and has received high attention and uptake in BPM practice [START_REF] Recker | BPMN Modeling-who, Where, How and Why[END_REF]. Besides, it is an established standard for business process modelling in industry and economy that is frequently supported by a computer program which enables a quite easy graphical description of complex processes [START_REF] Scheuerlein | New Methods for Clinical Pathways-business Process Modeling Notation (BPMN) and Tangible Business Process Modeling (t. BPM)[END_REF].

Furthermore, the analysis of performance of BPMN is crucial, as it helps Business analysts to predict whether the goal can be achieved or not. For this reason, simulation has been identified as a key technique for BP performance analysis [START_REF] Antonacci | A BPMN-Based Automated Approach for the Analysis of Healthcare Processes[END_REF]. It can be important for business processes as it helps and supports the decision-making process, reduces cost that can occurs in the case wrong decisions haven't been anticipated and as final objective provides a good quality of services.

As argued in [START_REF] Antonacci | A BPMN-Based Automated Approach for the Analysis of Healthcare Processes[END_REF], the use of BPM in healthcare sector is becoming a key enabler for the improvement of healthcare processes, since the healthcare environment is becoming more dynamic and volatile, and follows more complex processes, the combination of agent based model with business process may be efficient for resolving simulation's limitations in BP in terms of resource allocation, in the other hand, managing the availability of resources in healthcare sector is a very challenging problem with very little research attention. For this reason, the use of business process modeling combined with the adoption of simulation-based analysis provides a cost effective, accurate, and rapid way to evaluate alternatives before committing the required effort and resources [START_REF] Tumay | Business Process Simulation[END_REF][START_REF] Nakatumba | Business Process Simulation: How to Get It Right[END_REF].

The points we will explore in this paper are the concept of coupling processes with simulated models and also to propose a viable general approach based on modeling and simulation to address the problem of managing the resources allocation by applying it for healthcare case, so, the general purpose of this paper is to present a new architecture for modeling and simulation that aims to couple tasks followed by actors in the process with agents behaviors described with DEVS (Discrete Event System Specification). Overall, the reminder of this paper is structured as follows first related work and background is presented, then we have the contribution section when we first expose the problem statement, then we give a general overview of the proposed architecture and the methodology, after we will apply the architecture to ensure its feasibility.

RELATED WORK AND BACKGROUND

BPMN is a graphical notation for drawing business processes (OMG, 2003). It provides a standard notation that is easily understandable by all stakeholders; and also bridges the communication gap that frequently occurs between business process design and implementation. Nowadays BPMN becomes widely used by different organizations, as it is simple to learn, but yet powerful enough to depict the potential complexities of business processes. One of the main complexities of BP is the limitation of its simulation, which is due to the following reasons: lack of simulation know-how of BP analysts, costs and difficulties in retrieving and analyzing the data required for simulation model parameterization, large semantic gap between the business process model and the simulation model and finally, the use of models that may be (partially) incorrect or may not be at the right level of abstraction [START_REF] Bocciarelli | A Language for Enabling Model-Driven Analysis of Business Processes[END_REF]. In this context, the use of BPMN models in healthcare sector helps to facilitate the management of complex hospital BPs [START_REF] Antonacci | A BPMN-Based Automated Approach for the Analysis of Healthcare Processes[END_REF]. Accordingly, a lot of researchers in the literature, professionals and health administrators have given various definitions of patient pathways. In [START_REF] Benabdejlil | Modélisation Des Processus de Soins: Vers Une Implantation de Nouveaux Services à Valeur Ajoutée[END_REF], authors define a patient pathway as a description of process's elements of the disease management by the actors of the care network; it corresponds to a trajectory composed of steps that vary not only according to the actors of the network who can interact but also according to the place of residence and the medical history of the patient.

In the other hand, Multilayer Network has become a vital tool for understanding different relationships and their interactions in a complex system; it is often composed of many single-layer networks with nodes linked within edges. Each layer is sharing the same set of nodes, and edges are divided into two categories: intra-layer networks within each layer and inter-layer edges between vertices of different layers [START_REF] Murata | Comparison of Inter-Layer Couplings of Multilayer Networks[END_REF]. A graph (i.e., a single-layer network) is a tuple G = (V, E), where V is the set of nodes and E V ×V is the set of edges that connect pairs of nodes (Bollobas, 2013). If there is an edge between a pair of nodes, then those nodes are adjacent to each other. This edge is incident to each of the two nodes, and two edges that are incident to the same node are also said to be "incident" to each other [START_REF] Kivelä | Multilayer Networks[END_REF]. In the field of network science, two types of dynamic are studied: dynamics on networks and dynamics of networks. Dynamics on networks indicates that the activity (or some other property of interest) of nodes changes as a function of time and dynamics of networks indicates that the edges of the network themselves change either in their existence/absence or in their strength [START_REF] Bassett | Network Neuroscience[END_REF]. Multilayer Networks are useful for complex systems as they allow better understanding system's structure, in addition, many real-world systems do not operate in isolation.

In addition, simulation proved its value in the manufacturing sector and has been used to evaluate process problems in healthcare as well (Harrell and Price, 2000). It has been used for over a decade by health services for improving patient care. In this respect, there are a lot of examples in the literature which provide simulation tools for healthcare organization, in [START_REF] Laskowski | Agent-Based Simulation of Emergency Departments with Patient Diversion[END_REF], authors have developed a tool for managing emergency department, by planning capacity resources utilization and staff capacity. Moreover, authors in [START_REF] Norouzzadeh | Modeling to Optimize Healthcare Delivery in an Outpatient Clinic[END_REF] show how modelling and simulation of internal medicine practice process can help on decision making. The results of the simulation, based principally on patient waiting time were able to give an idea about improving resources utilization. Another example is given by Gűnal and Pidd in [START_REF] Gehlot | From Data to Processes??? Use of Modeling and Simulation in Healthcare[END_REF] which describes a model of the process flow of patients, that represents the multitasking behavior of medical staff (doctors and nurses), the only issue is that the ignorance of other possible factors such as: doctors' interactions with patients, other medical staff and their working environment lead to incomplete consideration of the problem (Jain et al., 2011).

Moreover, popular approaches include decision analysis, Markov process, mathematical modeling, systems dynamics and discrete event simulation [START_REF] Fone | Systematic Review of the Use and Value of Computer Simulation Modelling in Population Health and Health Care Delivery[END_REF] [START_REF] Kanagarajah | An Exploration into the Uses of Agent-Based Modeling to Improve Quality of Healthcare[END_REF]. The limitation of these approaches is that they ignore the effect of naturalistic human decision-making and behaviors on the performance of healthcare processes. In order to overcome such limitation, the use of Agent Based Modeling (ABM) offers complementary perspectives to model the process of health care domain [START_REF] Wainer | Discrete-Event Modeling and Simulation: A Practitioner's Approach[END_REF]. Agent based modeling (ABM) is formed by a set of autonomous agents that interact with their environment and other agents through a set of internal rules to achieve their objectives [START_REF] Onggo | Running Agent-Based Models on a Discrete-Event Simulator[END_REF] [START_REF] Grundspenkis | An Overview of the Agent Based Systems for the Business Process Management[END_REF].

In this context, there are some studies In the literature using BPMN model with DEVS model, these studies are based on meta-model approach (OMG 2003), which is one of the most used transformation techniques that includes the mapping of BPMN concepts to DEVS concepts. Based on the proposed approach of BPMN to DEVS (Cetinkaya, Verbraeck, and Seck, 2012) which presented a Model Driven Development (MDD) framework for modeling and Simulation (MDD4MS), and where a set of transformation rules were defined: some BPMN concepts (Pools, Lane, SubProcess) were mapped to DEVS coupled component while task, event (start, end and intermediate) and Gateway were mapped to DEVS atomic component [START_REF] Bazoun | SLMToolBox: Enterprise service process modelling and simulation by coupling DEVS and services workflow[END_REF]. This proposal doesn't cover the intervention of different resources like (human resources, devices and/or Software services) which may affect the execution of the task. To overcome such limitations, authors in (D'Ambrogio and [START_REF] Bouanan | Simulating Information Diffusion in a Multidimensional Social Network Using the DEVS Formalism (WIP)[END_REF] proposed new approach, by introducing reliability analysis that takes into consideration unexpected failures of the resources that execute the process tasks (unavailability of a resource allocated to task).

Discrete Event system Specification (DEVS) [START_REF] Zeigler | Theory of Modeling and Simulation: Integrating Discrete Event and Continuous Complex Dynamic Systems[END_REF] is formalism for modelling Discrete Events Systems. The hierarchical and modular structure of DEVS allows defining multiple models that are coupled to work together in a single and model by connecting their input and output through messages [START_REF] Wainer | Discrete-Event Modeling and Simulation: A Practitioner's Approach[END_REF]. In the same way, the resulting model can also be coupled with other models defining multiple layers in the hierarchical structure. In DEVS, atomic models define the behaviour of the system, and coupled models describe the structure of the system. The DEVS formalism has several advantages. The reason behind using DEVS is that it is based on dynamical systems theory and provides well defined concepts for coupling components, hierarchical and modular model construction, and an object-oriented substrate supporting repository reuse. Modular construction is one of the most important characteristics of DEVS because it allows the modeler to design and construct each model independently for optimal efficiency. As long as models adhere to certain protocols, they can interact which each other [START_REF] Pérez | Modeling and Simulation of Nuclear Medicine Patient Ser vice Management in DEVS[END_REF]. Accordingly, we use VLE (Virtual Laboratory Environment) for the implementation.

VLE is an open source software and API under GPL which supports multimodelling and simulation by implementing the DEVS abstract simulator [START_REF] Quesnel | The Virtual Laboratory Environment-An Operational Framework for Multi-Modelling, Simulation and Analysis of Complex Dynamical Systems[END_REF]. It integrates specific models developed in most popular programming languages into one single multi-model. It also proposes several simulators for particular formalisms; for example, cellular automata, ordinary differential equations (ODE), differential equations, various finite state automata (Moore, Mealy, Petri-Nets, etc.) and so on. VLE can be used to model, simulate, analyse and visualize dynamics of complex systems. Its main features are: multimodeling abilities (coupling heterogeneous models), a general formal basis for modelling dynamic systems and an associated operational semantic, a modular and hierarchical representation of the structure of coupled models with associated coupling and coordination algorithms, distributed simulations, a component-based development for the acceptance of new visualization tools, storage formats and experimental frame design tools [START_REF] Bouanan | Simulating Information Diffusion in a Multidimensional Social Network Using the DEVS Formalism (WIP)[END_REF]. The purpose of this paper is to present architecture for modeling and simulation in healthcare domain, the proposed architecture aims to overcome BPMN limitations and drawbacks, so for this we propose new approach which includes agent -based technologies. In the coming sections, we give an overview about the architecture, and then we study its feasibility in the case of emergency department (ED).

CONTRIBUTION

After giving an overview about the issues covered in this paper and before tackling the use case, it is required to introduce the general problem statement of our contribution. In this section, we will provide the problem statement of our study and then explain the proposed architecture and methodology.

Problem Statement

The use of workflow in healthcare domain is thriving, but since healthcare environment is evolving and complex to manage, it faces some drawbacks including: lack of coordination, limited flexibility during process enactment [START_REF] Bolcer | Advanced Workflow Management Technologies[END_REF], inability to cope with dynamic changes in resource levels and task availability, limited ability to predict changes, due to external events, in both the volume and composition of work, lack of performance, scalability and reliability as well [START_REF] Pang | Implementation of an Agent-Based Business Process[END_REF].

In the other hand, healthcare systems are complex to manage, and they face a lot of problems, poor coordination across healthcare levels is considered to be one of the main obstacles to attaining effective healthcare in many healthcare systems around the world, leading to difficulties in access to care, poor technical quality, discontinuity of care and inefficiencies in the use of resources (Organization and Control, 2008). Patient coordination problem can be insured by a representation of patient pathway where we can clearly see different interactions between resources. In this context, patient pathway has several definitions in the literature, an overall definition is given by ("ARS," 2016), that considers the patient pathway as the comprehensive and continuous care of patients closer to their place of life. It is a question of coordinating the interventions of the actors of the prevention and sanitary. The organization on the pathway helps to respond to the growing development of chronic diseases and situations of loss of autonomy that impact the patient's daily life. Patient pathways are based on the intervention of actors: of the health care system (doctors, pharmacists, hospitals, etc.), medical and social services and institutions (for the elderly or disabled), local authorities (town halls), other state services and social protection agencies. It's about organizing a coordinated care: throughout the patient's illnesses a territory and in the proximity of his home.

Overall, a multilayer representation may be helpful for having a meaningful representation of the system under study. Also, for meeting the individual needs of different roles of the patient pathway, BPMN can be considered the most efficient one, since it provides a good understanding of healthcare process [START_REF] Sang | BPMN Security Extensions for Healthcare Process[END_REF]. Besides, Patient pathway can be represented by BPMN as it offers an understandable representation of the system and its interactions, the only issue is that this representation is not dynamic, and as healthcare sector is known as dynamic system that evolve rapidly we propose a new architecture that uses BPMN elements for orchestrating the process where each element has its own behavior described by DEVS model. In order to well describe the architecture, we give a detailed description of the proposed approach and study its feasibility with a small number of data in the case of French healthcare system.

DEVS and BPMN for Healthcare Approach

In this section we will explain the used approach for modeling and simulation which is applied for healthcare sector, our main objective is to propose a solution for combining BPMN with dynamic models in order to manage availability of resources. For this, we choose a region for the study, specify participants and collect relevant data, then we draw the BPMN of the chosen case study for ensuring coordination and specifying exchanged information between resources. Interactions between entities are divided according to a multilayer network described in [START_REF] Sbayou | DEVS modelling and simulation for healthcare process application for hospital emergency department[END_REF].

General Overview of the Architecture

In Figure 1, we give an overview about the proposed architecture, our aim is to be able to use exchanged information between participants in the BPMN diagram for launching the DEVS network process, for this we extract the XML file of the BPMN model, and we read each of its component, then according to the message flow contained in the xml file, we create an action in the network. In the other hand, participants in the pathway are described as node in the multilayer network and as agents in DEVS model [START_REF] Sbayou | DEVS modelling and simulation for healthcare process application for hospital emergency department[END_REF]. Besides, healthcare resources' attributes are taken from database that contains information collected from yellow pages, some population's characteristics are generated according to specification given by the National Institute of Statistics and Economic Studies [START_REF] Ferret | Habitants En Nouvelle-Aquitaine Au 1er Janvier[END_REF], and geographical positions of patient are generated randomly. The agent-based simulation uses a MySQL database and an XML file generated from BPMN process for inputs. DEVS based simulation for modeling agent's behaviors and R-Studio tool for visualizing results. This approach allows us to study the dynamic part of BPMN model by integrating its elements in the global model that describes resources' behavior, the global DEVS model is a couple model of xml atomic model and the model that represent integration between participants in the BPMN. We draw an overview about the proposed architecture in Figure 1.

METHODOLOGY

As healthcare processes are complex to manage and simulate, we chose a category of model that can be used by non-specialist of M&S, and that is divided into 3 blocks.

The proposed modeling system is described in Figure 2, and it is divided as below: • Pre-Simulation: We define a conceptual model of the studied issue using BPMN. The BPMN process contains pools where each pool describes a participant in the workflow. Also, the structure of the organization of the service is described under a multiplex social network, as described in [START_REF] Sbayou | DEVS modelling and simulation for healthcare process application for hospital emergency department[END_REF]); • Simulation: The process starts by transforming the BPMN model to a formal model. Each agent is specified as an atomic model. BPMN model is then specified as a coupled DEVS model. DEVS is used for simulation which is run using VLE; • Post-Simulation: We use R to process the results. R is a programming language for statistical computing and graphics. RStudio tool is used to visualize the simulation results and to conduct analysis.

Operation of the XML Manager

As mentioned previously, we generate the xml file of the BPMN model for using it on orchestrating connection between patient and healthcare resources, and also it serves as input for simulation, the BPMN is used to ensure coordination between resources at first, and then an.xml file is generated where we will extract tasks that would be able to decide which resources to choose. For this, we create a library of DEVS atomic model, this model will read the xml file where each information is described by tags, the atomic model of the XML controller is describe in (Figure 3), the model has two states, "Init" and "Phase_0", at first, the model send an out information to other model according to the healthcare resources' profession, the output information are responsible for checking the available resources and establish connection between patient and available resources who are going to participate in his pathway.

CASE STUDY

In this section we will study the feasibility of our approach in the case of emergency department for a region located in southwest of France, this application is for proving the connection between actors' tasks that are contained in the BPMN process, and the system's behavior described by DEVS.

Healthcare System in France

The French healthcare system covers both public and private hospitals, doctors and other medical specialists who provide care services to French resident. It is accessible for all residents, independently of their age, income or status. French resident has to register a General Practitioner "Medecin Traitant", in order to ensure full eligibility to reimbursement of health costs. While following this process, the General Practitioner (GP or G), becomes the principal route which follows the patient care pathways "parcours de soin". In some cases, the patient can directly go to specialists without passing through GP. These specialists are: gynecologists, ophthalmologists, dentists and psychiatrists (for patients under 26 years old). In emergency cases, there is no appointment needed, in addition, emergency departments (ED) are the most complex system in healthcare sector, they offer care for life-threatening and other emergency cases over 24 hours daily, 7 days a week, accordingly, they usually require a lot of resources, at any time. The shortage of resources may lead to the recurrent problems such as overcrowding, congestion and also patient's dissatisfaction. For this reason, we propose to study the feasibility of our system for emergency departments and then apply it for other health cases.

Data Collection

Accordingly, we have retrieved geographical information about required resources in ED. As the number of existent resources is very important, we have located about 96 resources in a small sub region of "Nouvelle Aquitaine." Generally speaking, we extracted information about: Hospitals (H), Radiologists(R) and General Practitioner (GP and R not needed in the case of ED), we also generated population according to id, age, latitude and longitude parameters. The Figure 4 shows different existent resources, the orange bubbles represent patients with emergency needs, green babbles localize general practitioners, yellow ones are hospitals in the chosen region, and finally red bubbles represent radiologists in the chosen region. The database gives the positioning of the different resources that will be used at simulation time to compute their reachability. Rules of proximity threshold have been defined to be used at simulation time when a resource is reached by the patient himself or through the healthcare social network by the GP or the ED, in the ED, the process doesn't require nor GP nor R, and in emergency cases the patient is supported by ED without any appointment or recommendation.

Patient Pathway in ED

For our case study, we take the example of the University Hospital Center (CHU of Bordeaux).

Considering the CHU of Bordeaux; we have collected some statistical information about the patients as well as their care in the emergency department. The hospital is visited daily by about 149 patients, the completeness of reception is 100%, the sex ratio representing the number of men and number of women at birth during a given age group is 1,17. 1% of patients visiting the hospital are under one year old, 14% are patients under 18 years of age, patients aged 75 years and over account for 16%, and 7% of patients are not from the region. The average duration of process for each patient is about 4 hours according to ("ORU," 2017). In order to well understand steps followed by the patient in ED, we drew a BPMN diagram that describes roles of each participant, we can see that we have 2 pools in Figure 5, one for patient and the second one is for emergency department, in the case of ED in CHU of Bordeaux, the principal resources are the secretary, nurses and emergency physicians, hence, once the health issue occurs, the patient selects the GP or ED according to his state, for emergency cases the patient visits the hospital without any appointment, we can see that we have in Figure 4, colored tasks (orange, red and green), this tasks are described by tags information in the xml file, the orange tasks are called decisional tasks, like the task named "select hospital" in the patient pool, this task is responsible for launching the connection between the patient and the hospital, in the second pool named hospital, we have the secretary who is responsible of registration and providing patient's information to nurses, the available nurse starts the process by the triage of patient, the queuing system follows the first in first out (FIFO) option, red tasks in hospital pool are responsible of changing the state of resources, so, once the nurse is in red tasks he will be busy and not able to work with another patient, the nurse has to select the available physician that defines the priority of patients: resuscitation, treatment that lead to hospitalization, or examination which may be followed by hospitalization or redirection to home, in the physician lane we have green task named redirect home that means that resources are free and can be used by the next patient.

RESULTS AND DISCUSSION

In order to apply the architecture, we take the example of 10 patients who arrived to hospitals. First we extract the nearest hospital (proximity rule) to each one, when the hospital receives patients, and after registration, the process of checking of available resources is launched according to availability that some patients are only registered because there are not enough health resources which may take them in charge. We can also observe that that administrative agent is connected to all coming patients that are still waiting for their turn in order to join the network. This study ensures the feasibility of connecting BPMN and DEVS model for sample simulation of the allocation of available resources, the main idea is to test and verify the proposed architecture, then to apply the architecture for a big system with a lot of exchanged information and a complex queuing system.

CONCLUSION AND PERSPECTIVES

This paper presents a new architecture based on business process and multidimensional network for modeling in healthcare sector. It uses BPMN to define the patient pathway process connected with an ABM network. It defines and run the search for resources and their allocation based on a multi-layer different healthcare stakeholders' networks exploration. The objective is to study the impact of dynamic allocation of participant actors in BPMN healthcare pathways. Such an approach can be proposed to study a territory in terms of sufficient or insufficient resources available in a specific area. The case study, located in Bordeaux city, has been possible thanks to data coming from yellow pages that contains information about healthcare resources with their geographic location. The presented case study aims to manage resources and shows their impact on the performance of the process. As future work, we will include progressively new parameters as new resources (equipment), as well as more demographic information of patients and staff, budget constraints in order to manage healthcare system. We also propose to observe the impact of giving new geographical positions for healthcare resources, as well as to propose the implementation of new services, and run different scenarios to see how the system could be more performant. Finally, we want to apply our architecture in other sectors in order to make it functional for different application domain.

Figure 1 .

 1 Figure 1. General overview of the architecture

Figure 3 .

 3 Figure 3. Basic atomic model of XML controller

Figure 4 .

 4 Figure 4. Gironde department healthcare resources positioning

Figure 5 .

 5 Figure 5. Patient pathway in emergency cases


J.-J Marigo
email: marigo@lpmtm.univ-paris13.fr

L Truskinovsky

Initiation and propagation of fracture in the models of Griffith and Barenblatt

Keywords: fracture mechanics, variational methods, gamma-convergence, pull-out test, Barenblatt model PACS: 62.20.Mk, 68.35.Md, 81.40.Jj, 83.50.Lh

In the setting of the simplest debonding problem we give a systematic comparison of the fracture models due to Griffith and Barenblatt. We prove that the Griffith model represents an asymptotic Γ -limit of the Barenblatt model, when the ratio of the external and internal lengths increases indefinitely. We then illustrate the character of convergence by solving explicitly two sample problems with "initially rigid" and "initially elastic" cohesive energies. The geometrical simplicity of the setting allows us to study the small parameter dependence of both global and local minimizers of the total energy.

Introduction

While the well-known theories of fracture due to Griffith [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF] and Barenblatt [2] take into account cohesive forces in different ways, both have the same basic structure. The theory of Griffith assumes that the surface energy density is equal to a constant that does not depend on the crack opening, while the theory of Barenblatt postulates that the surface energy density is a monotonic function of the displacement discontinuity, approaching a constant only at some sufficiently large crack opening. For consistency of the two approaches, the constant appearing in the Griffith model and representing the toughness of the material must be equal to the limiting cohesive energy in the Barenblatt model. When the two models are consistent and when Barenblatt's cohesive forces act over sufficiently short range, a formal asymptotic analysis of [START_REF] Willis | A comparison of the fracture criteria of Griffith and Barenblatt[END_REF] showed that the stress fields near the tip of the crack predicted by the two models are asymptotically equivalent. In the present paper we explore the detailed character of this asymptotic convergence in the most elementary shear-lag setting of the fiber pull-off problem. The simplicity of the problem opens the way to studying the limiting behavior of both global and local minimizers without any a priori specification of the fracture geometry. An equivalent mathematical treatment of the general fracture problem presents a considerable challenge because of the complex technical issues associated with the variation of the generic discontinuity sets.

It is well known that the experimental characterization of cohesion parameters is usually based on either fiber pulling or pushing tests making the pull-out problem an important point of reference in the analysis of fiber-matrix interactions [START_REF] Marshall | Analysis of fibre debonding and sliding experiments in brittle matrix composites[END_REF][START_REF] Charalambides | Debonding properties of residually stressed brittle-matrix composites[END_REF]. The shear-lag approximation and its dynamic extensions were derived from the three-dimensional theory in [START_REF] Marshall | Analysis of fibre debonding and sliding experiments in brittle matrix composites[END_REF][START_REF] Geymonat | Stress distribution in anisotropic elastic composite beams[END_REF][START_REF] Hutchinson | Models of fibre debonding and pullout in brittle composites with friction[END_REF][START_REF] Müller | Static and dynamic pull out of an elastic rod from a rigid wall[END_REF][START_REF] Kerans | Theoretical analysis of the fibre pullout and pushout tests[END_REF][START_REF] Bechel | Application of debond length measurements to examine the mechanics of fiber pushout[END_REF][START_REF] Bilteryst | Amorçage de la décohésion dans l'essai d'arrachement[END_REF][START_REF] Bilteryst | An energy based analysis of the pull-out problem[END_REF][START_REF] Cox | Dynamic effects in the pull-out mechanism during dynamic loading of a bridged crack[END_REF]. The main conclusion of these and other investigations is that the shear-lag model can be used for sufficiently long and thin fibers and that the solutions are accurate except in the vicinity of the crack tip. Although the shear-lag approximation may misrepresent some fine features of the elastic fields, the error is known to decrease with the slenderness parameter in a controlled fashion. In the present paper we compare the Griffith and Barenblatt theories by explicitly solving the shear-lag problem for a composite bar in a hard device.

The main focus of our analysis is the influence of the structure of the cohesive energy on the geometry of the discontinuity set. The deficiencies of the simplest Griffith model in predicting the crack geometry are well known. For instance, since the ultimate strength of the material in this theory is infinite, it fails to describe the infinitesimal fracture (see e.g. [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF]). The problem disappears in the Barenblatt formulation, in which the ultimate strength is finite (see e.g. [START_REF] Del Piero | Macro-and micro-cracking in one-dimensional elasticity[END_REF]); to achieve the same result in the Griffith theory, one needs to introduce an incremental energy minimization (see e.g. [START_REF] Bilteryst | An energy based analysis of the pull-out problem[END_REF][START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF]). It remains unclear, however, how close are the predictions of the two theories concerning the maximally attained stress and the fine structure of the resulting fracture surface. Of particular interest is the exact domain of applicability of the much simpler and therefore more attractive Griffith approximation.

We show rigorously that the Griffith model furnishes the exact variational limit for a generic Barenblatt model when either the length of the bar tends to infinity or the internal length scale, originating from the cohesive law, tends to zero. We then predict the onset of debonding and compute the complete history of damage evolution until the full rupture. In addition to confirming asymptotic results by [START_REF] Willis | A comparison of the fracture criteria of Griffith and Barenblatt[END_REF], we obtain a detailed account of the character of convergence for both the energy and the geometry of the crack. Thus, we show that the "good" quantities that are continuous in the limit include the energy and the critical stress, but not the geometrical parameters of the damage zone such as, for instance, the length of the debonding. The Barenblatt theory is shown to play a major role when the crack is about to appear; in the corresponding range of loadings the convergence of the two theories is not uniform.

Contrary to most other investigations of fracture focused on the global minimum of the energy only, we explicitly study the local minimizers. The only other comparison of the local minimizers in the theories of Griffith and Barenblatt has been previously performed for a bar in tension (see e.g. [START_REF] Del Piero | Macro-and micro-cracking in one-dimensional elasticity[END_REF][START_REF] Del Piero | A one-dimensional model for localized and distributed fracture[END_REF][START_REF] Braides | Variational formulation of softening phenomena in fracture mechanics: the onedimensional case[END_REF][START_REF] Charlotte | Revisiting brittle fracture as an energy minimization problem: comparison of Griffith and Barenblatt surface energy models[END_REF]). The analysis of the associated one-dimensional problem shows that, while Griffith's global minimizer asymptotically approximates Barenblatt's global minimizer, the two theories generate rather different predictions concerning the structure of the local minima of the energy. Our pull-out problem for the fiber in a hard device is both simpler and more complex than the problem for a bar in tension. It is more complex because the displacement field in both the matrix and the fiber must be considered. It is more simple because, as we show, the bifurcation from the uncracked state is always supercritical so that the only local minimizer is the global one. The special structure of the bifurcation diagram in the pull-out problem is the result of the unavoidable stress concentration around the boundary of the domain of application of the kinematic boundary condition which creates a permanent nucleation site.

The paper is organized as follows. In Sect. 2 we introduce the shear-lag energy functional and set the formal variational problem by defining the appropriate space of competitors. We then reduce the vector problem to the scalar one and prove that both local and global energy minimizers can have at most one crack. In Sect. 3 the problem is nondimensionalized and the set of dimensional parameters is reduced to two main nondimensional criteria. We then prove the main theorem establishing the Griffith theory as an exact Γ -limit of the Barenblatt theory when one of the two nondimensional parameters tends to zero. In Sect. 4 we explicitly solve the local and global minimization problem in the Griffith theory and reduce the corresponding solution to quadratures in the general Barenblatt theory. In Sect. 5 we illustrate the general solution of the Barenblatt problem by explicitly solving two special cases: with "initially rigid" cohesive response (Dugdale model) and with "initially elastic" cohesive response (piece-wise quadratic model with zero ultimate strength). Finally in Sect. 6 we present a sketch of the theory in the case of a soft device. The paper ends with a discussion of the parameter sensitivity and of the size effect. All technical proofs of the mathematical statements are collected in the Appendix.

Main symbols

R

Radius of the cross-section R f

Radius of the fiber S m Area of the matrix cross-section S f Area of the fiber cross-section

Preliminaries

Consider a circular cylindrical composite bar of length L which includes a fiber of radius R f centered in the matrix of inner radius R f and outer radius R. Both the fiber and the matrix are isotropic, linearly elastic materials with Young moduli E f and E m , respectively, and with the same Poisson ratio ν. Both the matrix and the fiber are assumed to have infinite strength so that the fracture may only occur along the common interface. Before loading, the fiber is assumed to be glued to the matrix and the interface is assumed to be breakable. We adopt the following boundary conditions:

-the lateral part of the boundary (r = R) is free; -at one end (z = L) the matrix is free, while the fiber is loaded in a hard device generating a given longitudinal displacement; -at the other end (z = 0) both the matrix and the fiber are fixed.

In order to compute the equilibrium response we adopt the energy minimization point of view. Specifically, we require that, at each value of the loading parameter, the displacement field in the composite bar must deliver at least a local minimum to the corresponding energy functional. To give the precise definition of the minimization problem, we have to define -the set of admissible displacement fields; -the expression of the energy associated with these displacement fields; -the norm in the set of admissible displacements allowing one to evaluate the closeness between different displacement fields.

To simplify the analysis we consider a limiting case of sufficiently slender bars and approximate the 3D problem by a 1D problem (see [START_REF] Geymonat | Stress distribution in anisotropic elastic composite beams[END_REF] for more details). The longitudinal displacement fields inside the fiber, u f , and inside the matrix, u m , are assumed to depend only on the longitudinal coordinate z. It is convenient to take as unknowns the pair u = (u f , u d ), where u d = u f -u m . The difference of displacements u d vanishes if the bond remains unbroken and accordingly we can define the two sets

D(u) = {z ∈ [0, L] : u d (z) = 0}, C(u) = {z ∈ [0, L] : u d (z) = 0}. (1)
The displacements u f and u d must satisfy the following kinematic boundary conditions:

u f (L) = U, u f (0) = u d (0) = 0, (2)
where the prescribed displacement U of the end L of the fiber plays the role of the loading parameter. The total energy associated with the admissible displacement fields u = (u f , u d ) contains two terms: the elastic energy E e (u) and the cohesive (or surface) energy E s (u). The elastic energy can be written as

E e (u) = 1 2 L 0 E f S f u f (z) 2 + E m S m (u f (z) -u d (z)) 2 dz, (3)
where

S m = π(R 2 -R 2 f
) and S f = πR 2 f represent the area of the matrix and fiber cross-sections. For the cohesive energy we assume

E s (u) = L 0 2πR f κ(|u d (z)|)dz, (4)
where the presence of the absolute value means that the energy cost of sliding is independent of its sign. We consider the following two main choices for κ:

• Griffith model: κ is discontinuous at zero and is constant elsewhere

κ(δ) = 0 , δ = 0 γ , δ = 0 (5) • Barenblatt model: κ is a continuous function of δ such that κ(0) = 0, κ(δ) > 0 when δ > 0, κ(∞) = γ. ( 6
)
In both cases γ represents the (tangential) toughness of the interface. For the Barenblatt model σ c = κ (0) represents the ultimate shear stress, which may either be finite ("initially rigid" cohesive response) or zero ("initially elastic" cohesive response); in the Griffith model the ultimate strength is equal to infinity. We remark that the model with finite ultimate strength appears as the natural limit of a lattice model with Lennard-Jonestype interactions (see e.g. [START_REF] Truskinovsky | Fracture as a phase transition[END_REF]), while the model with zero ultimate strength can be linked to the presence of pre-existing fracture surfaces (see e.g. [START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF]). Now notice that the total energy of the bar,

E(u) = E e (u) + E s (u), (7)
is finite provided that u f and u d belong to the space W 1,2 (0, L); the set of admissible displacements then consists of pairs of functions belonging to this functional space and satisfying the boundary conditions [START_REF] Barenblatt | The mathematical theory of equilibrium cracks in brittle fracture[END_REF]. This set will be denoted by U. The precise statement of the local minimization problem reads

find u ∈ U such that ∃h > 0, ∀v : v -u ≤ h, E(u) ≤ E(v). ( 8
)
To complete the setting of the mathematical problem we choose the norm of W 1,2 (0, L) as defining the notion of a (strong) local minimum and denote it by • . We emphasize that this assumption is of an entirely physical nature. Weaker norms, considered for instance in [START_REF] Del Piero | Macro-and micro-cracking in one-dimensional elasticity[END_REF], would mean a broader set of local minimizers but would not affect the choice of the global minimum of the energy (and our main Theorem in Sect. 3.2).

Reduction to a scalar problem

We begin by eliminating u f and formulating the problem in terms of u d alone. Since E is a strictly convex function of u f , we can minimize it out by using the following Euler-Lagrange equation:

(E f S f + E m S m )u f (z) -E m S m u d (z) = const = F, (9)
meaning that the normal force is constant along the entire bar. By using the boundary condition u f (0) = u d (0) = 0 we can integrate (9) to give

u f (z) = A m A u d (z) + F A z. ( 10
)
Here we introduced the following notation:

A f = E f S f , A m = E m S m , A = A f + A m . ( 11
)
From the boundary condition u f (L) = U , we obtain

F L = A U -A m u d (L). (12)
By inserting these relations in the expression of the energy, we can rewrite [START_REF] Hutchinson | Models of fibre debonding and pullout in brittle composites with friction[END_REF] as

Ê(u d ) = A m A f 2A L 0 u d (z) 2 dz + A 2L U - A m A u d (L) 2 + 2πR f L 0 κ(|u d (z)|)dz. (13)
Here u d belongs to W 1,2 (0, L) and must satisfy the boundary condition u d (0) = 0. The energy Ê is positive and hence bounded from below. Since it is also weakly lower semi-continuous, there always exists a global minimizer of the semi-linear minimization problem for [START_REF] Cox | Dynamic effects in the pull-out mechanism during dynamic loading of a bridged crack[END_REF] and therefore of the original problem (8).

Monotonicity of the minimizers

Proposition 1. Assume that the surface energy function κ is an increasing function. Then the local minimizer of the energy (13) u d is a monotonic function of z, increasing when U ≥ 0, and decreasing when U ≤ 0.

Proposition 2. Let u be a local minimizer of the energy and C(u) be the set of points where the fiber and the matrix are bonded. If the surface energy κ is an increasing function, then

C(u) is a closed interval of the form [0, ], with 0 ≤ ≤ L.
Proof. The proof of Proposition 1 is given in the Appendix, whereas Proposition 2 is a consequence of the monotonicity of u d and of the boundary condition u d (0) = 0.

3 Relationship between the models of Griffith and Barenblatt

Dimensionless formulation

For the minimization of Ê it is convenient to introduce dimensionless quantities. We define

z = Lz, u d (z) = Lū d (z), U = A m A L Ū, F = A m L L F , a = A m A f , ( 14
)
where

L = 4πR f γA A f A m L = 4γA R f E f A m L ( 15
)
is the main internal length scale of the Griffith model. A nondegenerate Barenblatt model contains another internal length scale δ c = γ/κ (0), which can be used to form a dimensionless parameter

= δ c L . ( 16
)
In terms of the nondimensional variables the cohesive energy in the Barenblatt model can be written as

κ(u d ) = γ κ( ūd ). (17)
Now, if we normalize the total energy

2 Ê(u d ) = A f A m A L 2 L Ē(ū d ), (18)
we obtain

Ē(ū d ) = 1 0 ū d (z) 2 dz + 1 0 κ ūd (z) dz + a ūd (1) -Ū 2 . ( 19
)
The relation between the dimensionless force F and the dimensionless total displacement takes the form

F = Ū -ūd (1). (20)
As we see, the general Barenblatt problem contains two dimensionless parameters, namely a, which represents the relative stiffness of the matrix with respect to the fiber, and , which measures the ratio of internal (cohesive) and external (geometrical) length scales.

Limiting behavior as → 0

Since in the Griffith theory the nondimensional parameter is identically equal to zero, it is instructive to study how the structure of the minimizers changes as tends to 0. One of the goals is then to show that the longer the Barenblatt composite bar is, the closer its response is to that of a composite bar with the Griffith cohesive energy.

In what follows the main parameters of the Barenblatt theory will be explicitly marked by the subscript , while all other unnecessary subscripts will be dropped. For instance, we assume that ū delivers a local minimum to the energy functional Ē ,

Ē (v) = 1 0 v (z) 2 dz + a(v(1) -Ū ) 2 + 1 0 κ |v(z)| dz, ( 21
)
among all v from W 1,2 (0, 1) such that v(0) = 0. A local minimizer in the Griffith model can then be denoted by ū0 . By using Proposition 2 we can rewrite the corresponding energy functional Ē0 in the form

Ē0 (v) = 1 0 v (z) 2 dz + a(v(1) -Ū ) 2 + 1 -length {z ∈ [0, 1] : v(z) = 0} . ( 22
)
The relationship between the two theories is established by the following theorem:

Theorem 1. Suppose that a > 0, Ū ∈ R, > 0, and κ is a bounded positive function which satisfies κ(0) = 0 and lim δ→∞ κ(δ) = 1. Then the global minimizer of Ē converges (strongly in the sense of the W 1,2 (0, 1)-norm) to the global minimizer of Ē0 when goes to 0.

Proof. The proof is given in the Appendix.

We emphasize that this theorem is rather general and remains true even if κ is not monotonic or not continuous. The only limitation is that it ensures the convergence of the global minimizer and does not say anything about the local minimizers. As we show in the next section, in the case of a general Barenblatt surface energy, the unique local minimizer is also the global one, and therefore our Theorem ensures that the response of the Barenblatt tends to the response of the Griffith bar unconditionally.

Exact solution of the minimization problem

As a part of the proof of Proposition 1 we have shown that, when U = 0, u = (0, 0) is the unique local minimizer of E. Also, by symmetry, if u is a local minimizer corresponding to U , then -u is a local minimizer corresponding to -U . Therefore in what follows we consider only the case U > 0.

Griffith model

Owing to its simplicity, the Griffith problem can be solved in closed form. Recall that for a local minimizer

ū the set C(ū) = {z ∈ [0, 1] : ūd (z) = 0} is an interval of the form [0, (ū)] with 0 ≤ (ū) ≤ 1.
In the subspace of admissible displacements vd in W 1,2 (0, 1) such that vd (0) = 0 and C(v) = [0, ], the energy functional of the Griffith theory reduces to

Ē( , vd ) = 1 v d (z) 2 dz + (1 -) + a Ū -vd (1) 2 . ( 23
)
Since this functional is strictly convex with respect to vd for fixed , its unique (global) minimizer is given by

ūd (z) =      0 , 0 ≤ z ≤ a(z -) Ū 1 + a -a , ≤ z ≤ 1 . ( 24
)
To determine the value of we insert ( 24) into ( 23) to obtain Since Ě is strictly convex with respect to , there exists one and only one local minimizer of this one-dimensional problem which is also the global one. The minimum is reached at 0 when Ě (0) ≥ 0, at 1 when Ě (1) ≤ 0, and at such that Ě ( ) = 0 otherwise. Therefore, we obtain three different phases of behavior:

Ě( ) = Ē( , ūd ) = a Ū 2 1 + a -a + 1 -. ( 25
=          1 , a Ū ≤ 1 (phase I) 1 + 1 a -Ū , 1 ≤ a Ū ≤ 1 + a (phase II) 0 , 1 + a ≤ a Ū (phase III) . ( 26
)
The global response of the composite bar can also be illustrated by the following relation between the force F and the displacement Ū :

F =          Ū , 0 ≤ Ū ≤ 1/a (phase I) 1/a , 1/a ≤ Ū ≤ 1 + 1/a (phase II) Ū 1 + a , 1 + 1/a ≤ Ū (phase III) , ( 27
)
which is plotted in Fig. 1. Notice that the force remains constant in the phase of progressive debonding of the fiber.

The overall behavior of the composite bar in the Griffith model can be described by the language of the theory of phase transformations. In fact, the undamaged state with energy Ē1 ( Ū ) = Ě(1) = a Ū 2 (our phase I) and the completely debonded state with energy Ē2 ( Ū ) = Ě(0) = 1 + a Ū 2 /(1 + a) (our phase III) can be viewed as two different "phase modifications" of the "averaged material" of the composite bar. Now, since in the Griffith model these two "phases" can coexist without interaction, we can write for the mixed state (our phase II)

E(x, Ū1 , Ū2 ) = Ě( ) = x Ē1 ( Ū1 ) + (1 -x) Ē2 ( Ū2 ), (28)
where x = /L is the fraction of the length of the bar occupied by the first phase. The constraint of the total displacement of the fiber can be written in the form

x Ū1 + (1 -x) Ū2 = Ū. (29)
Now by minimizing the energy (28) under the constraint (29) we obtain the standard conditions of "phase equilibrium" furnishing the convexification of the total energy. The resulting overall response can be obtained through the common tangent (Maxwell) construction illustrated in Fig. 2. For a similar description of a broader class of partially damaged materials see [START_REF] Francfort | Stable damage evolution in a brittle continuous medium[END_REF]. Since in the Griffith problem the external length scale dominates, one can get an impression that in this theory the size effect is absent. This is not exactly so. Indeed, we have seen that the interface entirely debonds when the dimensionless force F reaches the value 1/a. In terms of the physical parameters this corresponds to a critical normal stress σ d given by (see [START_REF] Del Piero | Macro-and micro-cracking in one-dimensional elasticity[END_REF])

σ d = E f 4πR f γA A f A m = 2E f 1 + E f R 2 f E m R 2 m i R f . ( 30
)
Here the parameter i = γ/E f with a dimension of length represents the ratio between the energy of the surface (γ) and the energy of the bulk (E). The critical debonding stress σ d is therefore proportional to the dimensionless ratio i /R f of the two (internal) length scales: i , associated exclusively with toughness, and R f , giving the inner radius of the fiber. We conclude that the thicker the fiber, the lower the debonding stress.

Some important features of the debonding process will be different if instead of a hard device we consider a composite bar in a soft device. To illustrate the new possibilities we sketch in Sect. 6 the theory of a pull-out test for a Griffith fiber which is free at z = 0 and stretched by a force F at z = L. In particular, we find that there exists a limiting load above which no local minimizer exists. This means that the composite bar cannot support sufficiently large forces and, if overloaded, the fiber completely debonds. It is also interesting that if a < 1, the debonding is first initiated at z = 0 and if a > 1, it starts at z = 1.

Barenblatt model

Now assume that the surface energy is an arbitrary continuous and monotonically increasing function of δ satisfying the conditions κ(0) = 0, κ(δ) > 0 when δ > 0, and lim δ→∞ κ(δ) = γ > 0. To determine the local minimizers in the general Barenblatt model we cannot use Euler-Lagrange equations because they require the surface energy function to be differentiable, eliminating for instance an important case of the Dugdale theory (see below). An alternative approach is to use a Noether identity representing a first integral of the Euler-Lagrange equations.

Proposition 3. Let f : R → R be a continuous function, U ∈ R, and a ≥ 0. Then a field u ∈ W 1,2 (0, 1) such that u(0) = 0 is a local minimum (in the sense of the norm of W 1,2 (0, 1)) of the functional I defined by

I(u) = 1 0 u (x) 2 + f (u(x))dx + a(u(1) -U ) 2
only if it satisfies the following first integral condition and boundary conditions:

∃C ∈ R : -u (x) 2 + f (u(x)) = C, ∀x ∈ (0, 1), aU = u (1) + au(1).
Proof. The proof is given in the Appendix.

By applying Proposition 3 to a general Barenblatt model, we obtain that if ūd is a local minimizer of Ē it necessarily satisfies the following set of equations: ūd is increasing, (31a)

ū d (z) 2 -κ ūd (z) = C on (0, 1), (31b)
ūd (0) = 0, (31c) ū d (1) + aū d (1) = a Ū. (31d)
Due to the monotonicity of ūd , we know that C(ū) = [0, ]. It is convenient to consider the cases = 0 and > 0 separately.

• > 0. In this case, since ūd vanishes on [0, ], the constant C is necessarily zero. Then we obtain

ū d (z) = κ ūd (z) 1 2 on [ , 1]. (32)
This differential equation can be solved explicitly once is known. Define Ū1 = ūd [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF]. By using (32), (31c), and (31d), we obtain that and Ū1 must satisfy

1 -= g 0 ( Ū1 ) ≡ Ū1 0 κ u -1 2 du, (33a)
Ū = f 0 ( Ū1 ) ≡ 1 a κ Ū1 1 2 + Ū1 . ( 33b
)
Since the continuous function f 0 is strictly increasing from 0 to ∞ as Ū1 grows from 0 to ∞, (33b) has a unique solution Ū1 > 0 for any given Ū > 0. Concerning (33a) we must consider two possibilities depending on the behavior of the function κ(δ) -12 near δ = 0:

1. If the integral in (33a) is divergent (the case of a model with zero strength; see below), then (33a) has no solution and the minimizer must be in the other class = 0 (and C > 0). 2. If the integral is convergent, the function g 0 increases continuously from 0 to ∞ with Ū1 . Hence (33a) admits a unique solution ∈ (0, 1] provided that Ū1 ∈ [0, Ū c 1 ), with Ū c 1 = g -1 0 (1), which in turn requires that the loading Ū lies in [0, Ū c ) with Ū c = f 0 ( Ū c 1 ). Once Ū1 and are determined, the field ūd can be found explicitly: ūd (z) = 0 for z ∈ [0, ], and ūd (z) = g -1 0 (z -) for z ∈ [ , 1]. To summarize, for 0 < Ū < Ū c , there exists a unique function ūd satisfying the necessary conditions of optimality (31a)-(31d) and having length(C(ū)) > 0. When Ū ≥ Ū c , this type of a minimizer does not exist.

• = 0. Now consider local minimizers with an entire debonding of the interface. In this case, since ūd (0) = 0, we have C = ū d (0) 2 ≥ 0. We consider separately the cases C = 0 and C > 0. We still set Ū1 = ūd (1).

(i) C = 0. We can use (33a) and (33b) with = 0 and obtain that Ū must satisfy 1 = g 0 ( Ū1 ) and Ū = f 0 ( Ū1 ). That requires the integral of κ(δ) -1 2 near δ = 0 to converge. Then Ū = Ū c ≡ f 0 • g -1 0 (1) and the field ūd can be written explicitly as ūd (z) = g -1 0 (z) on [0, 1].

(ii) C > 0. In this case we get ū d (z) = κ ūd (z) + C 1 2

on [0, 1]. This differential equation can be solved as soon as C is known. By using (31b)-(31d), we obtain that C and Ū1 must satisfy

Ū = f (C, Ū1 ) ≡ Ū1 + 1 a C + κ Ū1 1 2 , ( 34a
)
1 = g(C, Ū1 ) ≡ Ū1 0 κ u + C -1 2 du. ( 34b
)
Consider the function g(C, Ū1 ) for C > 0 and Ū1 > 0. At fixed C, g is strictly increasing with Ū1 from 0 to ∞, while, at fixed Ū1 , it is strictly decreasing from the value g 0 ( Ū1 ) -which can be infinite -to zero. Now, at fixed C > 0, (34b) admits a unique solution Ū1 = G(C). Moreover, G(C) increases from the value g 0 -1 (1) -with the convention that g 0 -1 (1) = 0 when g 0 is not finite -to ∞ as C grows from 0 to ∞. Inserting the relation between Ū1 and C in (34a), we obtain the following equation for C:

Ū = f (C, G(C)). ( 35
)
Since, at fixed Ū1 , f increases from the value f 0 ( Ū1 ) to ∞ as C grows from 0 to ∞ and since, at fixed C, f increases from 0 to ∞ as Ū1 grows from 0 to ∞, we conclude that the function f (C, G(C)) increases from the value f (0, G(0)) = f 0 (g 0 -1 (1)) = Ū c to ∞ as C goes from 0 to ∞. Thus (35) admits one (and only one) solution C > 0 for a given Ū provided that Ū > Ū c .

Once Ū1 and C are determined, the field ūd can be obtained explicitly because ūd (z) is the unique solution of This solution is obviously a strictly increasing function of z. We have therefore found that for Ū ≥ Ū c there exists a unique function ūd satisfying the necessary conditions of optimality (31a)-(31d) and having length(C(ū)) = 0. Such a minimizer does not exist when Ū < Ū c . By comparing different cases, we can conclude that there exists at most one local minimizer. Since we have already shown that there is always a global minimizer, we may conclude that the only local minimizer is also the global one.

Explicit examples of the minimizers in the Barenblatt problem

To illustrate the general existence and uniqueness results for the Barenblatt model, here we present two particular examples. To cover different physical situations we consider two exactly solvable models with either "initially rigid" or "initially elastic" cohesive response (see Fig. 3). Specifically, we define 1. The Dugdale model ("initially rigid" response) :

κ(δ) =    γ δ δ c , 0 ≤ δ ≤ δ c γ , δ ≥ δ c . ( 37
)
2. The piece-wise linear model with zero ultimate strength ("initially elastic" response) :

κ(δ) =    γ δ δ c 2 , 0 ≤ δ ≤ δ c γ , δ ≥ δ c . ( 38
)

Dugdale model

Consider first the case where the surface energy is given by (37). The calculation of the integrals and solution of (34a) and (34b) is then straightforward. It remains to distinguish the points where the relative sliding ūd of the fiber is less than the critical value δ c (partial debonding) from those where it is larger (total debonding). The detailed succession of various phases of debonding is different depending on the value of , i.e. on whether the bar is long or short.

-Case of a long bar : < 1/2. We can distinguish three phases of debonding: 1. Phase I : 0 < Ū ≤ + 1/a. The bar is totally bonded along the interval [0, ] and partially debonded along the remaining part ( , 1). The relation between and Ū is given by

= 1 + 1 a - 1 a 2 + 4 Ū. ( 39
)
The relative sliding of the fiber takes the form

ūd (z) =    0 , 0 ≤ z ≤ (z -) 2 4 , ≤ z ≤ 1 . ( 40
)
At the end of this phase, i.e. when Ū = + 1/a, the relative sliding of the fiber reaches the critical value 1 at the loaded end (ū d (1) = 1), which marks the beginning of the second phase. During the first phase, the force F needed to sustain the displacement Ū is

F = 2 Ū 1 + √ 1 + 4a 2 Ū . ( 41
)
2. Phase II : + 1/a < Ū < 1 -+ 1/a. During this phase, the bar is split into three parts with perfect bonding along the interval [0, ], partial debonding along the interval ( , d), and total debonding along the remaining part [d, 1]. The two limit points are given by

= 1 -+ 1 a -Ū, d = + 2 . ( 42
)
Notice that the (dimensionless) length of the partially debonded part equals 2 . If we return to dimensional quantities, this length is independent of the length of the bar and is equal to

4γA R f E f A m -1 2 δ c . The relative displacement is now ūd (z) =        0 , 0 ≤ z ≤ (z -) 2 4 , ≤ z ≤ d z --, d ≤ z ≤ 1 . ( 43
)
At the end of this phase = 0, which means that the perfectly bonded domain ceases to exist, marking the beginning of the third phase. During the second phase, the force F remains constant:

F = 1 a . ( 44
)
3. Phase III : Ū ≥ 1 -+ 1/a. Since the end z = 0 of the fiber is fixed, there always remains a zone in which the debonding is partial. The length of this zone decreases to 0 as Ū goes to infinity. Outside this partially debonded zone, the debonding is total. From an analytical point of view, it is more convenient to use F as the parameter. The location of the limit point between the two zones is given by

d = 2 a F -a 2 F 2 -1 , ( 45
)
whereas the displacement field is now

ūd (z) =      (z -) 2 4 + z√ a 2 F 2 -1 , 0 ≤ z ≤ d a F z + d - ,d≤ z ≤ 1 . ( 46
)
Finally, the force F needed to sustain the displacement Ū is implicitly given by

Ū = (1 + a) F -+ 2 1 -a 2 F 2 + a F a 2 F 2 -1 . ( 47
)
One can show that F strictly increases from 1/a to ∞ as Ū grows from 1 -+ 1/a to ∞. The explicit inversion of this relation requires solving an algebraic equation of the third order. The three phases of the debonding process and the global response of the bar are illustrated in Fig. 4.

-Case of a short bar : ≥ 1/2. We can again distinguish three phases of debonding, the first and the last being the same as in the former case and the second corresponding to a case in which every point of the bar is partially debonded. The difference is due to the fact that the limiting point which separates the perfectly bonded zone from the partially debonded one reaches 0 before the totally debonded zone appears at z = 1.

               Phase I : 0 < Ū ≤ 1 4 + 1 2a
Phase II :

1 4 + 1 2a < Ū ≤ + a + 1 4a
Phase III :

+ a + 1 4a < Ū . ( 48
)
The relations giving , ūd , and F as functions of Ū during the first phase and those giving d, ūd , and Ū as functions of F during the third phase are the same as in the case of long bars. During the second phase, the whole bar is partially debonded, and the relative displacement of the fiber is now given by

ūd (z) = z2 4 + a 1 + a Ū - 1 4 - 1 2a z, (49)
whereas the force F is no longer constant but is strictly increasing:

F = 1 1 + a Ū + 1 4 . ( 50
)
The three phases of the debonding process and the global response of a short bar are illustrated in Fig. 4.

A piece-wise linear model with zero ultimate strength

Finally, consider the case in which the surface energy is given by (38) and is characterized by zero ultimate stress.

As we have already pointed out, in this model the integral g 0 in (33a) is divergent. As a result, the perfectly bonded phase does not exist: as soon as Ū > 0, all points of the fiber slide with respect to the matrix. As in the case of the Dugdale surface energy, we need to distinguish here the points where the relative displacement ūd is below the critical value δ c from those where it is above. We obtain two phases of the debonding process (Phases II and III):

ūd (z) = a Ū cosh 1 + a sinh 1 sinh z , (51)
F = Ū 1 + a tanh 1 . ( 52
)
Both increase with Ū . At the end of this phase, the relative displacement of the fiber at z = 1 reaches the critical value . At this point the totally debonded zone appears and the third phase begins. 2. Phase III : Ū > + 1/(a tanh(1/ )). Since the end z = 0 of the fiber is fixed, there always remains a zone in which the debonding is partial. The length of this zone simply decreases to 0 as Ū grows to infinity. Outside this zone, the debonding is total and the fiber does not interact with the matrix. From the analytical point of view, it is again convenient to use F as the parameter. The location of the limit point between the two zones is then given by

d = sinh -1 1 √ a 2 F 2 -1 , ( 53
)
whereas the displacement field takes the form

ūd (z) =        √ a 2 F 2 -1 sinh z , 0 ≤ z ≤ d + a F z -sinh -1 1 √ a 2 F 2 -1 , d ≤ z ≤ 1 . ( 54
)
During this last phase, the force F needed to sustain the loading Ū is implicitly given by

Ū = + 1 + a -a sinh -1 1 √ a 2 F 2 -1 F , (55)
One can show that F is a strictly increasing function of Ū . The explicit inversion of this relation requires the solution of a transcendental equation.

The two phases of the debonding process and the global response of the bar are illustrated in Fig. 5.

Soft device

For a bar in a soft device, the general energy functional takes the form

E(u) = 1 2 L 0 E f S f u f (z) 2 + E m S m (u f (z) -u d (z)) 2 dz + L 0 2πR f κ(|u d (z)|)dz -F u f (L), (56)
where the displacement of the matrix must satisfy the boundary condition u m (0) = 0, while the displacements of the fiber are free. For each value of F we are searching for the local minimizers of E satisfying the boundary conditions. By choosing the Griffith surface energy and following a procedure similar to the one described in Sect. 2, one can show that the dimensionless displacement discontinuity ū must minimize the energy Ē,

Ē(v) = 1 0 v (z) 2 dz -2 F (v(0) + av(1)) -a F 2 + 1 -length(C(v)). ( 57
)
Notice that, contrary to the case of a hard device, the configuration with F = 0 is not trivial. Since Ē(0) = 0 ≤ Ē(v), with equality only at v=0, the global minimum always corresponds to ū = 0. To see that this is not the only local minimizer, consider the fields ū = const = 0 corresponding to rigid translations of the fiber inside the unstretched matrix. All these fields have the same energy Ē(const) = 1. Since in any sufficiently small neighborhood of each one of these fields a perturbed field v does not vanish, the corresponding surface energy is equal to 1, whereas the bulk energy is non-negative. Consequently these "rigid translations" are all local minimizers and, as we prove below, there are no others. Moreover, although one can approach the trivial state by a sequence of such local minimizers, the limiting configuration v = 0 is separated from any of these local minimizers by an energy barrier of finite height (equal to 1). Another peculiar feature of the soft device is that the energy is no longer bounded from below when F = 0. Indeed, by taking v(z) = q F for all z, Ē(v) = -(2q + a) F 2 and tends to -∞ when q → ∞ provided that F = 0. Consequently, the global minimum does not exist as long as a nonzero force is applied. A minimizer then is necessarily a local one, if one exists at all. To find the local minimizers, we proceed as in Sects. 2 and 3. Proposition 4. Assume for determinacy that a < 1, meaning that the fiber is stiffer than the matrix. Then

ū (z) = F ( -z) , 0 ≤ z ≤ 0 , ≤ z ≤ 1 . ( 58
)
This family contains the state of perfect bonding ( = 0) as well as configurations with partial debonding of the fiber along the interval (0, );

4. If | F | > 1, there are no local minimizers.
Proof. The proof is given in the Appendix.

Discussion

In this concluding section we provide some additional arguments concerning the influence of the two main parameters characterizing the cohesive law -the toughness, measured by the internal length i , and the effective strength, measured by the internal length δ c -on the resulting equilibrium response of the composite bar.

The effective strength of the interface enters the description of the debonding process through the parameter . When the strength decreases, the parameter also decreases. Since this is equivalent to a transition from shorter to longer bars, we can conclude that the smaller the effective strength, the closer the response of the Barenblatt bar to that of a Griffith bar. The detailed character of convergence can be illustrated in the case of the Dugdale surface energy, for which the explicit relations for the minimizer ū and for the overall force F are available (see (39)-( 47)). Since we are interested in the limit of small toughness, it is sufficient to consider the case < 1/2. In particular, when 0 < Ū ≤ 1/a, we obtain from (39)-(41) that ū converges uniformly to 0 and that F converges to Ū , which corresponds to the first phase in the response of the Griffith bar (see [START_REF] Truskinovsky | Ericksen's bar revisited : energy wiggles[END_REF], (26), and (27)). When 1/a < Ū ≤ 1 + 1/a, we obtain from (42)-(44) that and d converge to 0 = 1 + 1/a -Ū , that ū converges uniformly to 0 on [0, 0 ] and to z -0 on [ 0 , 1], and that F remains constant and equal to 1/a. This corresponds to the second phase in the response of the Griffith bar. Finally, when Ū > 1 + 1/a, we obtain from (45)-(47) that d converges to 0, that ū converges uniformly to az/(1 + a), and that F converges to Ū/(1 + a). This corresponds to the third phase in the response of the Griffith bar. All these phases of the pull-out test are illustrated in Fig. 6, where the overall response of the Dugdale bars is plotted for different values of the effective strength of the interface against the response of the Griffith bar with the same toughness.

The value of the slope κ (0), characterizing the ultimate strength of the surface, also has an important influence on the character of the debonding process. Indeed, 1. In the Griffith model, in which κ (0) = ∞, there exists an interval of loading where the interface is entirely bonded; 2. In the model with zero ultimate strength, in which κ (0) = 0, a (partial) debonding appears along the entire bar as soon as the load is applied; 3. In the intermediate cases, when 0 < κ (0) < ∞, the debonding zone appears as soon as the load is applied, but it spreads only along a finite segment of the bar. The length 1of this segment increases with the loading. In particular, we can compute the derivative

d d Ū (0) = 2a Lκ (0) , (59)
confirming that the incipient debonding length is infinite when the slope κ (0) is zero ("initially elastic" model) and vanishes when the slope is infinite (Griffith model).

The above differences in the detailed character of the debonding process do not contradict our general conclusion concerning the convergence of the Barenblatt and Griffith models in the limit of zero effective strength. Thus, although in the piece-wise linear model with zero ultimate strength (but nonzero effective strength) a partial debonding appears along the entire bar as soon as the load is applied, one can show that the smaller the effective strength, the smaller the relative displacement of the fiber. In particular, in the limit of zero effective strength, the sliding disappears completely.

This discussion confirms that the important parameter determining the detailed structure of the response is the effective strength, not the ultimate strength. Notice, also, that if the pull-out test is performed for sufficiently long composite bars, then the only material parameter that can be extracted from the experiment is the (tangential) toughness of the interface which survives the transition from the Barenblatt to the Griffith description. The situation here is completely different from the one in a simple tension test, in which the structure of the set of local minimizers is much richer than in the present model and, therefore, in the absence of fluctuations (or imperfections) only the ultimate strength can be measured, the toughness playing a relatively minor role (see [START_REF] Charlotte | Revisiting brittle fracture as an energy minimization problem: comparison of Griffith and Barenblatt surface energy models[END_REF]). In that and related cases, the predictions of the Griffith model are not satisfactory and preference should be given to the Barenblatt model. while the bulk energies outside the interval [ 1 , 2 ] are equal, the energy stored in the interval [ 1 , 2 ] is less for the field vh than for the field ū except when ū = 0 in this whole interval, which means that C(ū) = [ 1 , 2 ]. Now examine the possibility that ū ∈ A. Consider vh = ū + hφ, φ ∈ W 1,2 (0, 1). For h sufficiently small, vh ∈ A and Ē(ū) ≤ Ē(v h ). This is possible only if the derivative of Ē at ū in the direction φ vanishes, meaning that 1 0 ū (z)φ (z)dz = F φ(0)+aφ [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF] . This in turn requires that F = 0 and ū = 0. We have therefore proved that a local minimizer describing the state of total debonding is possible only if F = 0 and that is necessarily a rigid translation of the fiber.

Finally, consider the possibility that C(ū) = [ 1 , 2 ]. For the corresponding displacement field v the energy

Ē(v) = 1 0 v (z) 2 dz + 1 2 v (z) 2 dz -2 F (v(0) + av(1)) -a F 2 + 1 -2 + 1
is a convex functional of v and therefore the local minimizer (ii) The local minimizers. Now we need to select among candidates the local minimizers. We have already seen that, when F = 0, all rigid translations of the fiber are local minimizers.

ū(z) =      F ( 1 -z) , 0 ≤ z ≤ 1 0 , 1 ≤ z ≤ 2 a F (z -2 ) , 2 ≤ z ≤ 1 (A6)
Consider now the candidates which belong to A. Recall that when | F | < 1 ū0 = 0 is the only candidate, and that when | F | = 1 all the candidates ū have the same energy as ū0 = 0: Ē(ū ) = Ē(0) = -a. Let x 0 ∈ [0, 1] and let B(x 0 ) = {v ∈ W 1,2 (0, 1) : v(x 0 ) = 0}. Notice that ū0 ∈ B(x 0 ) and that Ē admits a global minimizer in B(x 0 ): the functional is bounded from below and lower semicontinous in this closed subspace of W 1,2 (0, 1). If we call this (or one of these) minimizer ū(x 0 ) we can follow the steps of the analysis from (i) and obtain that C(ū(x 0 )) is necessarily of the form [ 1 , 2 ]. But we know also that ū0 is the global minimizer of Ē on the set B(x 0 ) when | F | ≤ 1 and therefore Ē(0) = Ē(ū(x 0 )). Since this property is independent of the choice of the point x 0 , it remains true on the union of the sets B(x 0 ), that is on the whole B. In other words, we can guarantee that when | F | < 1 Ē(0) ≤ Ē(v), ∀v ∈ B. Similarly, when | F | = 1 we obtain that Ē(ū ) ≤ Ē(v), ∀v ∈ B.

To complete the proof we need to compare the energies of the candidates from B and A. To prove that ū0 = 0 is really a local minimizer when | F | < 1, it suffices to show that for any sequence vn ∈ A converging to 0 we have lim n→∞ Ē(v n ) > Ē(0). Since for such sequences lim n→∞ Ē(v n ) = Ē(0) + 1 , ū0 is a local minimizer. To prove that ū is a local minimizer when | F | = 1, it suffices to show that for any sequence vn ∈ A converging to ū we have lim n→∞ Ē(v n ) > Ē(ū ). Since for such sequences lim n→∞ Ē(v n ) = Ē(ū ) + 1 -, ū is a local minimizer when 0 ≤ < 1. Now, for = 1, consider the sequence vn = ū1 + F /n ∈ A converging to ū1 . Since Ē(v n ) = Ē(ū 1 ) -2(1 + a) F 2 /n < Ē(ū 1 ), the candidate ū1 is never a local minimizer.

Fig. 1 .

 1 Fig. 1. Debonding process and global response in the Griffith model

Fig. 2 .

 2 Fig. 2. The Maxwell construction providing the overall response of the composite bar in the Griffith model.

Fig. 3 .

 3 Fig. 3. Different cohesive energy functions and the corresponding cohesive stress-opening relations

ab

  Case of a long bar : < 1/2. Case of a short bar : > 1/2.

Fig. 4 .

 4 Fig. 4. Debonding process and global response in the Dugdale model

Fig. 5 .

 5 Fig. 5. Debonding process and global response in the piece-wise linear model with zero ultimate strength

1 .

 1 If F = 0, the perfectly bonded state (ū = 0) and totally debonded states (ū = const) are the global and the local minimizers, respectively; 2. If | F | < 1, the perfectly bonded state (ū = 0) is the unique local minimizer; 3. If | F | = 1, there is a one-parametric family of local minimizers ū 0≤ <1 ,

Fig. 6 .

 6 Fig. 6. Convergence of the global response of the Dugdale model to the Griffith model as the effective strength of the interface decreases

  is necessarily the global one. Being projected on the family (A6) the energy becomes an affine function of 1 and 2 :Ē( 1 , 2 ) = (1 -F 2 ) 1 + (1 -a 2 F 2 )(1 -2 ) -a F 2 .Since ū depends continuously on 1 and 2 , it is a local minimizer only if the pair ( 1 , 2 ) delivers a global minimum to Ē( 1 , 2 ). Now, by recalling the assumption a < 1, we obtain that, if| F | < 1, then 1 = 0, 2 = 1; if | F | = 1, then 1 ∈ [0, 1], 2 = 1; and finally, if | F | > 1, then 1 = 2 = 1.In other words, we obtained that the local minimizers ū with C(ū) = ∅ must necessarily belong to the following set:ū0 if | F | < 1, the family {ū } 0≤ ≤1 if | F | = 1, and ū1 if | F | > 1.

du = z. (36)

Appendix

Proof of Proposition 1. First recall that (being in W 1,2 ) u d is continuous and that u d (0) = 0. Assume now that u d is not monotonic. Then there exist z 1 , z 2 , and u 0 such that

Let h be such that 0 < h < 1 and define v h d by

.

This displacement field is admissible: it satisfies the boundary condition v h d (0) = 0, it is as smooth as u d in each subinterval delimited by z 1 and z 2 , and it is continuous at z 1 and z 2 , where it takes the same value u 0 as the field u d . Let us note also that u

d tends to u d when h tends to 0. Hence, for sufficiently small h, it must have a larger energy than u d . But a straightforward computation shows that

Since κ is increasing, both terms on the right hand side are not positive and so they must vanish for the inequality to hold. That requires u d = 0 on (z 1 , z 2 ) and therefore u d (z) = u 0 on (z 1 , z 2 ). This is a contradiction, so u d is monotonic.

It remains to be proven that u d increases or decreases depending on the sign of U . Since u d is monotonic and vanishes at 0, it has the same sign on the whole interval [0, L]. Thus it suffices to prove that U and u d (L) have the same sign.

Let h > 0 and

If Uu d (L) < 0, the four terms on the right hand side are not positive and we obtain a contradiction. If Uu d (L) = 0 the inequality is satisfied if and only if u d = 0. Consequently, when U > 0, then u d (L) ≥ 0, which means that u d is increasing. Similarly when U < 0, then u d (L) ≤ 0, meaning that u d is decreasing, and when U = 0, then u d (L) = 0, and u d = 0 is the only possible local minimizer. In this latter case, it is easily checked that 0 is really the (unique) global minimizer.

Proof of the Theorem. The proof is divided into four steps. Essentially we are establishing a point-wise convergence of the Barenblatt functional to the Griffith functional and dealing with an elementary example of Γ -convergence (see, for instance, [START_REF] Braides | Gamma-convergence for beginners[END_REF]).

(i) We first verify the basic inequality. Let v be an element of W 1,2 (0, 1) such that v(0) = 0 and let v be a sequence of such elements weakly converging to v when → 0. Then we need to show that

To prove (A1), we must consider separately the three terms in the energy expression. For the bulk elastic energy, since 0 ≤

2 dz, we easily obtain that

For the boundary term, since v (1) converges to v(1) we get (v(1) -Ū ) 2 = lim →0 (v (1) -Ū ) 2 . To prove the required inequality for the surface energy let C(v) be the set of points where v vanishes and D(v) the set of points where it does not. Due to the positivity of κ, we have

This confirms that (A1) holds.

(ii) Now we need to ensure the existence of a recovery sequence. It suffices to show that, for every v in W 1,2 (0, 1), that

Notice first that only the surface energy depends on . Since the surface energy function vanishes at 0, we need to prove that lim →0 D(v) κ(|v(z)|/ )dz = length(D(v)). But this follows from the assumption that κ(|v(z)|/ ) converges to 0 if v(z) = 0 and to 1 if v(z) = 0. Since κ is bounded, we can apply the Lebesgue dominated convergence theorem and the result follows.

(iii) We can now prove that a global minimizer ū of Ē weakly converges to the global minimizer ū0 of Ē0 . Since Ē (ū ) ≤ Ē (0) = a Ū 2 , the sequence ū is bounded in W 1,2 (0, 1) and hence we can extract a subsequence (also denoted ū ) which weakly converges to a certain ū. Then, due to (A1), (A2), and the fact that ū is a minimizer of Ē , we obtain for any v

which ensures that ū is the unique global minimizer of Ē0 , and therefore that ū = ū0 . (iv) By choosing v = ū0 in (A3) we obtain that Ē0 (ū 0 ) ≤ lim →0 Ē (ū ) and hence that 1 0 ū (z) 2 dz converges to 1 0 ū 0 (z) 2 dz. That ensures that ū converges strongly to ū0 . Proof of Proposition 3. The analysis below is based on a standard use of inner variations (see [START_REF] Truskinovsky | Ericksen's bar revisited : energy wiggles[END_REF] for a somewhat different derivation in a more general case).

The minimizers must belong to the closed subspace of W 1,2 (0, 1) satisfying u(0

h , and u h converges to u when h goes to 0. Consequently, if u is a local minimizer of I, the inequality I(u h ) ≥ I(u) must hold for h sufficiently small. A straightforward calculation gives

By noticing that I(u • φ -1 h ) is now differentiable with respect to h and that it has a minimum at h = 0, we can conclude that its derivative at 0 must vanish. This leads to

Now it is straightforward to show that there exists a constant C such that u (x) 2 -f (u(x)) = C for all x ∈ (0, 1).

Let us derive the boundary condition. Let u be a local minimizer, h > 0, θ ∈ R, and define u h by

Notice that for h sufficiently small, u h is admissible and converges to u as h goes to 0; then we can obtain the following inequality:

Dividing by h and passing to the limit when h goes to 0 we obtain that 0 ≤ 2θ u (1) + au(1) -aU + θ 2 . This last inequality will be satisfied for all θ if and only if u (1) + au(1) = aU , which constitutes the desired boundary condition.

Proof of Proposition 4. The proof is divided into two steps. We first determine the short list of candidates by using necessary conditions of equilibrium and then select among them local minimizers. We adopt the following notation:

(i) The candidates. Suppose that ū is a local minimizer of (57). We first show that C(ū) is either empty or a closed interval of the form [ 1 , 2 ], with 0 ≤ 1 ≤ 2 ≤ 1. Indeed, if C(ū) is not empty, define 1 = min{z ∈ C(ū)} and 2 = max{z ∈ C(ū)}, C(ū) being closed because ū is continuous. Let h > 0 and choose vh to be the following admissible displacement field:

.

Since vh tends to ū when h goes to 0, we must have Ē(ū) ≤ Ē(v h ) for h sufficiently small. Notice, however, that the surface energy is equal for both vh and ū because they vanish at the same points. The energies of the loading device for the two fields ū and vh are also the same, because the fields are necessarily equal at 0 and 1. Now,