Zeina Azmeh
email: azmeh@lirmm.fr

Maha Driss
email: mdriss@irisa.fr

Marianne Huchard
email: huchard@lirmm.fr

Naouel Moha
email: moha@irisa.fr

Chouki Tibermacine
email: tibermacin@lirmm.fr

QoS-Driven Selection of Composable Web Services

Web Services are universally accessible software units that are advertised, discovered, and invoked over the internet. As web service technology is becoming widely adopted by organizations that need to integrate their information systems within and across organizational boundaries, it is important to identify automatically relevant composable services. In this paper, we propose an approach based on a variant of Formal Concept Analysis for the identification of composable services. It allows the classification of web services according to their QoS and composability levels and proposes a query mechanism that allows users to specify their required QoS and composability levels. In the case of multiple choices that satisfy users' requirements, we enhance our approach with a mechanism based on vectors for identifying the optimized choice. We validate our approach on a set of real-world web services obtained from Service-Finder.

Introduction

Service-Oriented Computing (SOC) is an emerging paradigm for developing lowcost, flexible, and scalable distributed applications based on web services [START_REF] Huhns | Service-oriented computing: Key concepts and principles[END_REF]. Web services are autonomous, reusable and independent software units that can be accessed through the internet. SOC is becoming broadly adopted and in particular by companies, which are more and more willing to open their information systems to their clients and partners over the Internet. The reason comes from the fact that SOC offers the ability to build efficiently and effectively addedvalue service-based applications by composing ready-made services. Web service composition addresses the situation when the functionality required by users cannot be satisfied by any available web service, but by assembling suitably existing web services [START_REF] Rao | A Survey of Automated Web Service Composition Methods[END_REF]. Discovering and selecting relevant services that closely fit users' functional and non-functional requirements is an important issue highly studied in the literature [START_REF] Zeng | Qosaware middleware for web services composition[END_REF][START_REF] Aggarwal | Constraint driven web service composition in meteor-s[END_REF][START_REF] Yu | Efficient algorithms for web services selection with end-to-end qos constraints[END_REF]. Functional requirements define functionalities provided by web services and non-functional requirements define Quality of Service (QoS) criteria such as availability, response time, security, and throughput, etc. [START_REF] Menascé | Qos issues in web services[END_REF]. However, discovering and selecting relevant composable services is another important issue that still need to be investigated, since few approches focused on service composability problem [START_REF] Medjahed | A multilevel composability model for semantic web services[END_REF][START_REF] Ernst | Detection of web service substitutability and composability[END_REF][START_REF] Lécue | Semantic and syntactic data flow in web service composition[END_REF]. Relevant composable services represent services minimize the amount of adaptation among them while best fitting requirements.

In this paper, we propose an approach for the identification of composable services that best fit QoS and adaptation requirements. This approach is based on the formal concept analysis (FCA) and its variant, the relational concept analysis (RCA). FCA has been successfully applied in the past as a formal framework for service substitution [START_REF] Peng | Management and retrieval of web services based on formal concept analysis[END_REF][START_REF] Aversano | Using concept lattices to support service selection[END_REF][START_REF] Azmeh | Wspab: A tool for automatic classification and selection of web services using formal concept analysis[END_REF][START_REF] Azmeh | Using concept lattices to support web service compositions with backup services[END_REF][START_REF] Chollet | Heterogeneous service selection based on formal concept analysis[END_REF]. Thus, we use FCA and RCA to classify services according to their QoS and composability level and suggest substitutable services and composable services. The approach is also based on a query mechanism that allows users to specify their required QoS and composability levels. In particular, users can specify their interests for identifying composable services that maximize QoS requirements and minimize their adaptation level. In the case of multiple choices that satisfy users' requirements, we enhance our approach with a mechanism based on vectors for identifying the optimized choice.

The paper is organized as follows. Section 2 describes our approach along with a motivating example. Section 3 describes the experiments performed on a real case study for validating our approach. The paper ends with the related work in Section 4 and the conclusion in Section 5.

Approach

We explain our approach along with an abstract orchestration of three sequential tasks. We try to identify the concrete services offering the needed functionalities, in order to instantiate the defined orchestration. Our objective is to select composable services with the highest affordable QoS values, and the least required adaptations. We have 3 sets of services: S 1i , S 2j and S 3k corresponding to the three tasks. For each service of the three sets, we assume certain values of QoS (for availability (A) and response time (RT)). We categorize the QoS values into 5 levels: very low, low, medium, high and very high, for each of the A and RT values. A service that offers a specific QoS level, does also offer QoS levels that are lower. For each pair of consecutive services (according to our defined orchestration), we assume three composability values: exact, partial and disjoint.

RCA-Based Classification. We base our classification on Formal Concept Analysis (FCA, [START_REF] Ganter | Formal Concept Analysis: Mathematical Foundations[END_REF]), which is a data analysis method that aims at extracting concepts from entities described by attributes. The description of entities is encoded in a formal context. For example, in Fig. 1, the third table S 3k is a formal context where entities are services for the Task 3, and are described by the values of availability and response time (attributes). Concepts are maximal sets of entities (extent) sharing a maximal set of attributes (intent), and they are organized in a classification provided with a lattice structure. The lattice associated with the context S 3k is shown in Fig. 5. Each concept in this lattice inherits the attributes of its ascendants (super-concepts), and has lesser attributes than its descendants (sub-concepts). Reversely, each concept inherits the services of its descendants (sub-concepts). For example, the concept c 0 represents the services s 31 and s 33 , which covers the QoS values of the services in the upper concepts, and offers better ones (like: low A and low RT). The concept c 5 indicates that s 33 has very low RT and high A. The concept c 5 is a sub-concept of the concept c 0 , which shows that s 33 has better QoS values than s 31 . Since we also want to capture the composition relations, existing between pairs of services in an orchestration, we use a variant of FCA called Relational Concept Analysis (RCA), which includes relations in the attributes used for the concept mining [START_REF] Huchard | Relational concept discovery in structured datasets[END_REF]. RCA data are given in the form of a Relational Context Family (RCF). An RCF is a set of binary relationships (tables): some (non-relational) tables describe entities by attributes, while some (relational) tables describe relations between entities. The non-relational part of the RCF contains one table for each task: entities are services that can make the task, and attributes are QoS values. Thus, our RCF is composed of three non-relational contexts: S 1i * QoS, S 2j * QoS and S 3k * QoS, and two sets of three relational contexts: composition of S 1i * S 2j and composition of S 2j * S 3k , according to the three considered levels.

As we can see in Fig. 1 We run the RCA algorithm on the RCF, illustrated in Fig. 1. The RCA takes as input all the tables and iterates on two steps: (1) building a concept lattice on each non-relational context of services concatenated with the corresponding relational contexts that express composition, that is on S 1i * QoS concatenated with Exact 12 , P artial 12 and Disjoint 12 , on S 2j * QoS concatenated with Exact 23 , P artial 23 and Disjoint 23 and on S 3k * QoS; (2) transform the six relational contexts to integrate the concepts found at that current step (to use this knowledge in the next iteration). During this transformation, we use the existential scaling operator which is one of the scaling operators provided by RCA [START_REF] Huchard | Relational concept discovery in structured datasets[END_REF]. Services that form the columns of a composition context are replaced by concepts that group services. If (S xy , S zt ) initially holds in a composition context (with composition level cl) and S zt is in the extent of a concept C at the current step, thus in the current version of the composition context we add (S xy , C). This can be interpreted as S xy can be composed with at least one service from the set of services grouped in C with the composition level cl. For example, since we have (s 13 , s 25 ) ∈ Exact 12 and since s 25 is in the extent of Concept c 1 (in lattice S 2j ), the scaled Exact 12 table will contain (s 13 , c 1 ). By regarding the lattices in Fig. 2, we can find for example, that the service s 13 (in the leftmost lattice) has in its intent Exact 12 : c1. This means that it can be composed with services contained in the concept c1 (in the middle lattice). Thus, s 13 can be composed at least with s 25 . At a next step this will allow to group a set of services (implementing a first task) because they can be composed with services of another group of services (implementing another task). The process stops when no new concept emerge during the FCA analysis. At the end of the process, a concept lattice is generated for each set of services. In each lattice, services are classi-fied into concepts, showing their QoS levels as well as their composability levels with other services, according to the orchestration. Concept c 7 in the S 1i lattice groups services s 13 and s 11 , and indicates that they can be partially composed with services of concept c 3 of the S 2j lattice, which groups s 23 and s 25 . Navigation by Query Integration. We define a query mechanism that enables us to navigate into a lattice, in order to choose a combination of services that satisfies the aimed orchestration, with our aimed levels of QoS and composability. A query determines the required QoS for each requested task, as well as, the needed composability level (required adaptations) between the services. We define a QoS query as a new line (a new entity), integrated into the corresponding service*QoS binary context, specifying the least expected QoS level.

In our example, we perform three QoS queries: query 1 , query 2 and query 3 on S 1i * QoS, S 2j * QoS and S 3k * QoS respectively. They are as follows: query 1 looks for a service from S 1i , with a low A and a low RT; query 2 tries to retrieve a service from S 2j , with a high A and a low RT; query 3 demands a service from S 3k , with a high A and a medium RT. In order to select services that satisfy a specific composability level, we merge the relational contexts that express levels equal or greater than the requested level. In other words, if a partial composability level can be accepted, this means that if we find an exact level, it is even better. In this example, we demand an "exact" composability level between S 1i and S 2j , and a "partial" level between S 2j and S 3k . This gives us a variation of our RCF that includes queries, where the relational context compose1 is the context exact12 and compose2 is the fusion of exact23 and partial23, as shown in Fig. 3. This latter RCF generates three lattices, from which we can extract the services combinations that meet our requested QoS and composability. These lattices are illustrated in Fig. 4, from which we can identify the services that satisfy the QoS queries, which appear in the sub-concepts of the concept where a query n appears, and are as follows: services s 13 and s 14 for query 1 ; services s 23 , s 24 and s 25 for query 2 ; services s 32 and s 33 for query 3 . The Optimal Services Composition. From the extracted services, we can have several combinations that satisfy our desired orchestration. If we are searching for an optimized selection, this means that we have to find the set of services that meet the optimal compromise of QoS levels and composability. Thus, services that can be coupled according to our acceptable levels of required adaptations and QoS. In order to meet this issue, we propose to represent the services by vectors. A vector per set of services, on which, services are ordered according to their QoS. We also define a vector for composability levels. Each composition can be regarded as a triangle, having a head corresponding to the composability level, and the two other heads corresponding to the pair of services to be composed, as shown in Fig. 6. The optimized choice would be the triangle that has the minimal area, in case of a two services composition. Otherwise, it will be the minimal sum of the services triangles according to an orchestration. Thus in our example, by regarding the triangles in Fig. 6, we notice that the triangle (E, s 25 , s 33 ) represents an optimized composition between S 2j and S 3k (corresponding to Task2 and Task3, respectively). Accordingly, if we consider the triangle (E, s 13 , s 25 ) that shares an edge with the previous triangle, it represents a composition between S 1i and S 2j (corresponding to Task1 and Task2, respectively). These two triangles together may be an optimized combination for the required orchestration, after comparing them with all the existing triangles4 .

Experiments

This section describes data and results of the experiments that we have conducted.

Experiments Data

To validate our approach, we consider an orchestration composed of three sequential tasks. The first task IPToCity returns the city name for a given IP address. The second task CityToZipCode returns the zip code for a given city name. The third task ZipCodeToWeather returns the weather information for a given zip code.

We discover a set of WSs that reply to each of these tasks by querying Service-Finder [START_REF] Valle | Realizing service-finder, web service discovery at web scale[END_REF] [19]. Service-Finder is a Web 2.0 platform for WSs discovery that indexes almost 25.000 WSs. For example, we use keywords 'IP + City' to search WSs that achieve the first task. Service-Finder returns a result set of 36 WSs5 .

To reduce this set, we process three-step filtration. In the first step, we eliminate repeated WSs (i.e., WSs that offer exactly the same operations). This first filtration generates a new reduced set that contains 29 WSs. These WSs are passed to a second filtration step, in which we check whether the endpoint URI exists or not. This produces a new reduced set that contains 21 WSs. These WSs are passed to a third filtration step which leaves only services offering operations that satisfy the IPToCity task. This is done by parsing every service WSDL in order to check if the required operation exists or not. Finally, we obtain a set of 5 WSs. Table 1 gives for the three tasks, the keywords used to query Service-Finder, the number of WSs returned by it, and the number of the obtained WSs after each filtration step. The total number of WSs that reply to the three tasks is 21 WSs6 . Table 2 is a synopsis of the 21 WSs used in our experiments, their operations, inputs, and outputs. For the sake of simplicity, we use the abbreviated notation 'S11' to refer to the service '1' that replies to the task '1' and 'Op11a' to refer to the operation 'a' from the service '1' that replies to the task '1', and so on. We consider 2 QoS properties to characterize WSs: availability (A) and response time (RT). Real time monitoring information of service availability and response time are provided by Service-Finder. To spread out availability and response time values, we use the boxplot statistical technique [START_REF] Chambers | Graphical methods for data analysis[END_REF].

Tasks

Table 3 gives availability and response time values for the 21 WSs. To characterize the composability between two WSs, for example S11 and S21, we compute the semantic and structural similarity between the output and input parameters of their operations. We obtained values (Sim) ∈ [0,1]. Then, we used a standard reasoning inferences [START_REF] Paolucci | Semantic matching of web services capabilities[END_REF], common in semantic web services: Exact(Sim = 1), Plugin (2/3 <= Sim <1), Subsume (1/3 <= Sim < 2/3), Intersection (0 < Sim < 1/3), and Disjoint (Sim = 0). Computation of semantic and structural similarity between parameter names is done using the Java WordNet Similarity Library (JWSL) [START_REF] Pirró | Design, implementation and evaluation of a new semantic similarity metric combining features and intrinsic information content[END_REF]. Due to the limited paper space, we show in Table 4 the obtained similarity values between Op11a of the service S11 and operations of the services that reply to CityToZipCode task.

Experiments Results

We use the previously obtained data in building an RCF, as described in Section 2. Thus, our non-relational contexts are composed of the obtained services together with their QoS values, illustrated in Table 3 constructed using the composability values in Table 4. Having two composable services, means that the first service contains at least one operation that can be composed with one or more operations from the second one. The maximal level of composability between the operations belonging to the two services, is used as the level of composability of the two services. For example, from the Table 4, we can see that the service s 11 has two values of composability with service s 29 , which are "Exact" and "Plugin". Accordingly, we consider the composability level between s 11 and s 29 to be "Exact". After building the RCF, the next step is to integrate queries, in order to specify the required QoS of the service corresponding to each task in the orchestration. We define the following queries: query1 demands services from S1i having a "Low A"and a "High RT"; query2 demands services from S2j having a "High A" and a "Medium RT"; query3 demands services from S3k having a "High A" and a "Low RT". We also specify the minimal accepted composability level, we set it to be "Plugin". Thus, the relational contexts become one single context, which is the fusion of the "Exact" and "Plugin" contexts for the two sets of relational contexts (S 1i ,S 2j ) and (S 2j ,S 3k ). Executing RCA on the built RCF gives us the relational concept lattices (RCL) in Fig. 7, in which, we have composable services with a "Plugin" composability level (or "Exact"), and having the expected QoS levels. By looking to query3 in L3, we notice that it asks for services having a "Medium RT" and "Medium A". While the original query3, requests a "High A" and a "Low RT", as indicated above. Since there were no matches for the original query, it switches to the next best level of QoS. The services returned by the queries are as follows: query1 returns s 12 in c3, and s 14 in c5 (lattice L1); query2 returns s 23 , s 26 in c13, s 28 in c21, and s 212 in c8 (lattice L2); query3 returns s 31 in c0 and s 34 in c4 (lattice L3).

In order to clarify the results obtained by the RCL, we show the composability values between the returned services in Table 5. We notice that, the services s 12 and s 14 have a "Plugin" composability with the services s 23 ,s 26 , and s 28 . Thus, s 12 and s 14 must have the concepts c13 and c21 in their attributes sets. By regarding the lattice L1 in Fig. 7, we notice that c3 and c5 inherits Compose(L2) : c13 and Compose(L2) : c21 from c2. This means that s 12 and s 14 are composable with s 23 ,s 26 , and s 28 , which is true. In the same way, we notice that s 23 ,s 26 , and s 28 are composable with s 31 , and s 34 in lattice L3. On the other hand, we can also notice s 212 in L2, which can be composed with s 31 , and s 34 , but can not be composed with neither s 12 nor s 14 . This explains the absence of Compose(L2) : c8 as an attribute for c3 and c5 in L1.

Related Work

The related work can be organized in three categories:

Web Service Composability and Substituability. Many works in the literature addressed the composability and substituability of web services. Ernst et al. [START_REF] Ernst | Detection of web service substitutability and composability[END_REF] present a method based on syntactic descriptions of Web services. This Table 5. Composability values between the services of query1 and query2, then services of query2 and query3.

approach analyzes multiple invocations of these services. The input and output parameter values are compared and matchings are deduced. Contrarily to our approach, this work is based on the experimental usage of Web services, and is perfectly complementary to ours. In this research area, much work has been done on semantic web services [START_REF] Sycara | Automated discovery, interaction and composition of semantic web services[END_REF][START_REF] Medjahed | Composing web services on the semantic web[END_REF][START_REF] Medjahed | A multilevel composability model for semantic web services[END_REF][START_REF] Lécue | Semantic and syntactic data flow in web service composition[END_REF]. In [START_REF] Medjahed | Composing web services on the semantic web[END_REF] and [START_REF] Medjahed | A multilevel composability model for semantic web services[END_REF], the authors present a model for checking the composability of semantic web services at different levels: syntactic, semantic and quality of service. They first proposed a model for describing the static (non-computational) and dynamic (business) semantics of web services. Then, they formalized a set of rules that allow the computation of a composability degree.This work deals with two kinds of composability: horizontal (normal composability) and vertical (substituability). In our approach, web service semantics are extracted from service's signatures and documentation. We use tools (based on WordNet) to compute the similarity between names or documentation keywords of different web services. This allowed us to deduce the composaibility of services. We do not use a high-level model of web service semantics as in [START_REF] Medjahed | A multilevel composability model for semantic web services[END_REF]. In addition, in our work we experimented the approach on real-world (not semantic) web services obtained from Service-Finder. The authors of [START_REF] Medjahed | A multilevel composability model for semantic web services[END_REF] made however a simulation-based experimental study. In [START_REF] Medjahed | Composing web services on the semantic web[END_REF], the same authors deal with three quality attributes: fees, security and privacy.

In our approach, we focus on other web service qualities, which are availabil-Fig. 7. The RCL corresponding to the services listed in Table 2.

ity and efficiency (response time). In [START_REF] Lécue | Semantic and syntactic data flow in web service composition[END_REF], the authors present a method which combines semantic and static analysis of messages in semantic Web services. Data types (XML schemas) and semantic description (domain ontologies and SA-WSDL [START_REF] Sivashanmugam | Adding semantics to web service standards[END_REF] specifications) of parameters (parts of web service messages) are used to deduce mappings. These mappings are then transformed into adapters (XSL documents). In our work, we concentrated on the selection of services which offer maximal QoS and which require, when composed, minimal adaptations. The two approaches are complementary. defines the problem as a multi-dimension multi-choice knapsack problem and the graph model defines the problem as a multi-constrained optimal path problem. A heuristic algorithm is introduced for each model. Despite the significant improvement of these algorithms compared to exact solutions [START_REF] Zeng | Qosaware middleware for web services composition[END_REF] [1], both algo-rithms are not suitable for large service-based systems. This type of systems (i.e., composed of a large number of services) may be composed easily and efficiently using our approach.

Web Service Selection using Concept Lattices. Many works in the literature have addressed the classification of web services using concept lattices [START_REF] Azmeh | Wspab: A tool for automatic classification and selection of web services using formal concept analysis[END_REF][START_REF] Azmeh | Using concept lattices to support web service compositions with backup services[END_REF][START_REF] Aversano | Using concept lattices to support service selection[END_REF][START_REF] Peng | Management and retrieval of web services based on formal concept analysis[END_REF][START_REF] Chollet | Heterogeneous service selection based on formal concept analysis[END_REF].

In [START_REF] Azmeh | Wspab: A tool for automatic classification and selection of web services using formal concept analysis[END_REF], Azmeh et al. presented a tool which helps in classifying directories of web services from Seekda into concepts lattices. The classification is built upon contexts which associate services to operation signatures. In [START_REF] Azmeh | Using concept lattices to support web service compositions with backup services[END_REF], the used contexts formalize relations between services and keywords extracted from their documentation, and similarity between operations. In these works, the classification is not based on QoS (filtering in [START_REF] Azmeh | Wspab: A tool for automatic classification and selection of web services using formal concept analysis[END_REF] is done on only one quality attribute, availability). In addition, these approaches are based on FCA, and not on RCA, and don't base the classification on an orchestration. In [START_REF] Chollet | Heterogeneous service selection based on formal concept analysis[END_REF], the authors present an approach where formal contexts associate web services and functionalities provided by these services as well as their non-functional properties (security). This interesting work addresses a different kind of quality attributes (security) that is critical in pervasive computing. Besides, it does not deal with service composition, and uses FCA (not RCA).

Conclusion

We presented an approach for assisting the development of a web-service orchestration based on a workflow including interacting abstract tasks, and on readymade services from public registries. Relational Concept Analysis (RCA) helped us to take into account both QoS requirements and levels of composability in classifying candidate web services for implementing the tasks. The result is a set of web service lattices (one for each task), where the designer can see potentially collaborating services, as well as their QoS values. Inside each lattice, services are ordered by QoS and composability values. We have also proposed a querying mechanism that, given an expected level of QoS and composability, guides the designer in the lattices and highlights parts of the lattices where services should be chosen. We led an experiment on a sequential orchestration composed of three tasks. This illustrated and showed the feasibility of the approach. As a future work, we plan to lead other experiments and to study other encodings of the problem with RCA, to give complementary views. In addition, we intend to study more deeply other composability computation techniques (based on other semantic descriptions of web services). This will allow us to obtain more precise compositions.

Fig. 1 .

 1 Fig. 1. The relational contexts family (RCF) corresponding to the three sets of services.

  , the three relational contexts for composition of S xi * S yj are called Exact xy , P artial xy and Disjoint xy .

Fig. 2 .

 2 Fig. 2. The services lattices for Task 1 and Task 2.

Fig. 3 .

 3 Fig. 3. The RCF after integrating the QoS queries and specifying the composability level.

Fig. 4 .

 4 Fig. 4. The service lattices with the integrated queries.

Fig. 5 .Fig. 6 .

 56 Fig. 5. The lattice associated with the context S 3k of Task 3. Fig. 6. Optimizing the service selection.

Table 1 .

 1 Summary of the obtained sets of services.

		Keywords	# WSs	# WSs	# WSs	# WSs
		used to query returned by obtained by the obtained by the obtained by the
		Service-Finder Service-Finder	1 st step	2 nd step	3 rd step
				filtration	filtration	filtration
	IPToCity	'IP +	36	29	21	5
		city'
	CityToZipCode	'City +	166	140	112	12
		Zip +
		Code'
	ZipCodeToWeather	'Zip +	9	9	7	4
		Code +
		Weather'
					Total	21

Table 2 .

 2 . They are S 1i * QoS, S 2j * QoS, and S 3k * QoS. The services inter-relation (relational) contexts are WSs used in our experiments.

	Tasks	WSs names	Operations names	# Inputs # Outputs
		&	&
		identifiers	identifiers
	IPToCity	IP	S11 GetCityNameByIP	Op11a 1	1
		DOTSGeoPinPoint	S12 GetLocationByIP	Op12a 2	9
			GetCountryByIP	Op12b 2	6
			GetLocationByIP V2	Op12c 2	12
		WsIp	S13 GetLocation	Op13a 1	3
		GeoCoder	S14 IPAddressLookup	Op14a 1	10
			PhysicalAddressLookup	Op14b 1	10
		XigniteHelp	S15 WhatsMyIP	Op15a 0	33
	CityToZipCode	GeoPlaces	S21 GetPlaceDetails	Op21a 2	12
		MediCareSupplier	S22 GetSupplierByCity	Op22a 1	14
		DOTSGeoCoder	S23 GetGeoLocationByCityState Op23a 3	5
			GetGeoLocationWorldWide	Op23b 4	5
		ClientScriptServices S24 GetZipCodeCompletionList Op24a 3	5
		FedAch	S25 GetAchByName	Op25a 3	5
		DOTSAddressValidate S26 ValidateCityStateZip	Op26a 4	4
			S26 ValidateAddressSingleLine Op26b 2	19
		ProMilesWebService S27 GeoCodeLocation	Op27a 6	6
			GetZipsInRadius	Op27b 2	4
			ValidateLocation	Op27c 4	4
		CityZip Service	S28 GetZipCode	Op28a 4	3
		USZip	S29 GetinfoByCity	Op29a 1	1
			GetInfoByState	Op29b 1	1
		DotSGeoCash	S210 GetATMLocationsByCityState Op210a 2	8
		LPWebService	S211 GetLocationCities	Op211a 2	1
			GetLocationStates	Op211b 1	1
			GetLocationByCity	Op211c 3	1
		GeoServiceWebService S212 GetZipCodesFromCityState Op212a 3	1
	ZipCodeToWeather DOTSFastWeather	S31 GetWeatherByZip	Op31a 2	18
			GetWarningsByZip	Op31b 2	11
			GetFiveDayForecastByZip	Op31c 2	8
			GetWeatherHistoricalByZip Op31c 4	18
		USWeather	S32 GetWeatherReport	Op32a 1	1
		WeatherForecast	S33 GetWeatherByZipCode	Op33a 1	13
		WeatherService	S34 GetWeatherData	Op34a 1	15

Table 3 .

 3 QoS values of the WSs used in our experiments.

	Task	Services	Availability Response time
	IPToCity	IP	High	Very high
		DOTSGeoPinPoint	Medium	Low
		WsIp	High	Very high
		GeoCoder	High	High
		XigniteHelp	Low	Medium
	CityToZipCode	GeoPlaces	High	Medium
		MediCareSupplier	Very low	High
		DOTSGeoCoder	High	Medium
		ClientScriptServices	High	High
		FedAch	Very low	Low
		DOTSAddressValidate	High	Medium
		ProMilesWebService	Medium	Medium
		CityZip Service	High	Medium
		USZip	Very low	Low
		DotSGeoCash	Medium	Medium
		LPWebService	High	Medium
		GeoServiceWebService	High	Low
	ZipCodeToWeather	DOTSFastWeather	High	Medium
		USWeather	Very low	Low
		WeatherForecast	Very Low	Low
		WeatherService	Medium	Medium

Table 4 .

 4 Similarity values and inferences between Op11a and operations achieving CityToZipCode task.

	Operation Operations achieving Similarity values
		CityToZipCode task	& inferences
	Op11a	Op21a	0.682	Plugin
		Op22a	1	Exact
		Op23a	0.802	Plugin
		Op23b	0.789	Plugin
		Op24a	0.291 Intersection
		Op25a	0.292 Intersection
		Op26a	0.699	Plugin
		Op26b	0.534 Subsume
		Op27a	0.460 Subsume
		Op27b	0.643 Subsume
		Op27c	0.424 Subsume
		Op28a	0.766	Plugin
		Op29a	1	Exact
		Op29b	0.710	Plugin
		Op210a	0.903	Plugin
		Op211a	0.599 Subsume
		Op211b	0.377 Subsume
		Op211c	0.772	Plugin
		Op212a	0.614 Subsume

  Web Service Selection wrt QoS. In[START_REF] Zeng | Qosaware middleware for web services composition[END_REF], Zeng et al. present a middleware platform that enables the quality-driven composition of web services. In this platform, the QoS is evaluated by means of an extensible multidimensional model, and the selection of web services is performed in such a way as to optimize the composite services QoS given a set of constraints and a set of candidate web services. Two service selection approaches for constructing composite services have been proposed: local optimization and global planning. Their study shows that global planning is better than local optimization. Aggarwal et al.[START_REF] Aggarwal | Constraint driven web service composition in meteor-s[END_REF] present a constraint driven web services composition tool that enables the selection and the composition of web services considering QoS constraints. Like Zeng et al.[START_REF] Zeng | Qosaware middleware for web services composition[END_REF], a linear integer programming approach is proposed for solving the optimization problem. In[START_REF] Yu | Efficient algorithms for web services selection with end-to-end qos constraints[END_REF],Yu et al. propose heuristic algorithms to find a near-to-optimal solution more efficiently than exact solutions. Yu et al. models the QoS-based service selection and composition problem in two ways: the combinatorial model

We haven't shown all the triangles in Fig.6, for the sake of simplicity.

The discovery result set is obtained on June 18 th , 2010.

The WSDL URLs are available on: http://www.lirmm.fr/∼azmeh/wise10/expe/wsdls.html


Martin Larcher
email: martin.larcher@jrc.it

Development of Discrete Cracks in Concrete Loaded by Shock Waves

Keywords: Element-free Galerkin method, Strain rate effect, Hugoniot curve, Detonation load, Discrete cohesive cracks

The non-linear behaviour of statically loaded concrete results mainly from the development of cracks. In addition, highly dynamically loaded concrete is influenced by the failure of micro pores. Several authors use a damage model to consider both effects. The objective of this work is to test the capability of a discrete crack model to simulate shock wave-loaded concrete. Discrete cracks are modelled using the element-free Galerkin method, which allows cracks to be implemented independently of the mesh. The discrete crack model is used in combination with a fracture process zone to simulate the cohesive behaviour of the cracks in concrete. The behaviour of the micro pores is described by a compression modulus, which depends on hydrostatic strain (Hugoniot curve). The strain rate effect is considered by modifying the failure surface for calculating crack development. Static and dynamic benchmarks with several concrete beams show good agreement with experimental results. Experiments with concrete loaded by contact detonation show the development of cracks below the crater. The simulation model used shows good agreement with the experimental crack pattern.

T P I R C S U N A M D E T P E C C A T P I R C S U N A M D E T P E C C A 1 Introduction
Security-sensitive buildings have to be designed against extreme dynamical loads. In addition to loads due to earthquakes and the impact of vehicles, restraining these structures against highly dynamic loads, mostly from terrorist attacks becomes more and more important. Detonation loads can be separated in contact detonations resulting from the detonation of an explosive which is in contact with the material and air blast loads resulting from a detonation at some distance. The main objective of this paper is the description of the behaviour of the concrete loaded by contact detonations.

Highly dynamic loading results in a high strain rate (up to 10 6 sec -1 ) and in high pressures. The aim of this work is to develop a simulation model for the failure of highly dynamic loaded concrete. A combination of a discrete crack model implemented with the element-free Galerkin method and a cohesive crack model is used.

The development of cracks results in the non-linear behaviour of concrete. A highly dynamic loading of concrete produces additional non-linear effects such as the failure of micro-pores and the increasing of the strength by increased strain rates.

Crack models for concrete

There are several methods that can be used to model the cracking of concrete. All can be divided into two main procedures: Discrete crack models and continua mechanical models. A comparison between the two procedures is given by Mosler [START_REF] Mosler | Embedded crack vs. smeared crack models: a comparison of elementwise discontinuous crack path approaches with emphasis on mesh bias[END_REF].

Examples of continua models are smeared cracks, plasticity and damage models, although only the latter two models are used often to describe highly dynamic loaded concrete. The cracks are represented in these models by decreasing stiffness. A classical combination of a plastic material law and the strain rate effect is given by Malvern [START_REF] Malvern | The propagation of longitudinal waves of plastic deformation in a bar of material exhibiting a strain-rate effect[END_REF] and is, for example, used by Tedesco [START_REF] Tedesco | A strain-rate-dependent concrete material model for ADINA[END_REF]. More sophisticated is the plasticity model of Malvar [START_REF] Malvar | A plasticity concrete material model for DYNA3D[END_REF], which is implemented in LS-DYNA and uses different failure surfaces for the failure modes. Several material parameters are used to fit the model to concrete. The hydrocode model of Holmquist, Johnson and Cook [START_REF] Holmquist | A computational constitutive model for concrete subjected to large strains, high strain rates, and high pressures[END_REF] looks only at highly dynamic effects and splits the strains into a hydrostatic and a deviatoric part. The deviatoric stresses are calculated by damage evolution, whereas the Hugoniot effect is considered by the calculation of the hydrostatic stress. The resulting stresses are multiplied with a strain rate factor. Similar models are proposed by Ockert [START_REF] Ockert | Ein Stoffgesetz für die Schockwellenausbreitung in Beton[END_REF], Riedel [START_REF] Riedel | Beton unter dynamischen Lasten -Meso-und makromechanische Modelle und ihre Parameter[END_REF] and in combination with SPH by Rabczuk [START_REF] Rabczuk | Simulation of high velocity concrete fragmentation using SPH/MLSPH[END_REF], [START_REF] Rabczuk | A three-dimensional large deformation meshfree method for arbitrary evolving cracks[END_REF].

T P I R C S U N A M D E T P E C C A

ARTICLE IN PRESS

All these methods need additional approaches for failed concrete. Several investigations use element erosion to avoid disturbed elements. Alternatively, the failed parts can be switched in a discrete crack model. A comparison of the methods used after failure is presented by Song [START_REF] Song | A comparative study on finite element methods for dynamic fracture[END_REF].

Other possibilities of considering the non-linear behaviour are discrete crack models, which implement the discontinuity with a real crack. Therefore, the position and the size of a crack have to be calculated. Discrete cracks can be implemented in the finite element method by using, for example, remeshing algorithms. Comacho [START_REF] Comacho | Computational modelling of impact damage in brittle materials[END_REF] uses a fixed mesh of triangles with cohesive cracks between the elements. An effective stress intensity factor for mixed-mode fracture acts as a criterion for the development of the crack.

In contrast with these methods there are further models that allow cracks between the nodes. The eXtended finite element method (X-FEM) proposed by Moes [START_REF] Moes | A finite element method for crack growth without remeshing[END_REF] uses additional degrees of freedom for the description of the crack. To obtain better results at the crack tips, Belytschko [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF] uses an asymptotical crack tip function with X-FEM. Pedersen [START_REF] Pedersen | Continous-discontinous modelling of dynamic failure of concrete using a viscoelastic viscoplastic damage model[END_REF] uses X-FEM with a viscoelasticviscoplastic damage model. The viscoplasticity describes the strain rate effect. In this study the element-free Galerkin method is used in combination with a cohesive crack model to describe the discontinuity of the cracks.

There are several ways of implementing cracks in mesh-free methods. The extended element-free Galerkin method (XEFG) proposed by Rabczuk [START_REF] Rabczuk | A meshfree method based on the local partition of unity for cohesive cracks[END_REF] is a combination of the local partition of unity method X-FEM and EFG. Degrees of freedom are added for all nodes concerned of a crack. Another possibility is the cracking particle method (Rabczuk [START_REF] Rabczuk | Cracking particles: a simplified meshfree method for arbitrary evolving cracks[END_REF], [START_REF] Rabczuk | A three-dimensional large deformation meshfree method for arbitrary evolving cracks[END_REF]), which adds crack planes through particles. The crack nodes are enriched to discretise the discontinuity. The method of phantom nodes (see Song [START_REF] Song | A method for dynamic crack and shear band propagation with phantom nodes[END_REF]) implements additional (phantom) nodes in each cracked element to describe the crack. A crack is then represented by two superimposed elements. The method can be used to describe shear bands. In this study, the element-free Galerkin method will be used to implement discrete cracks.

Shock wave loaded concrete

There are different possibilities of implementing the behaviour of the micropores in a constitutive relation for concrete. Some authors implement this effect in the form of an equation of state. Riedel [START_REF] Riedel | Beton unter dynamischen Lasten -Meso-und makromechanische Modelle und ihre Parameter[END_REF] proposes a constitutive law based on the model of Holmquist [START_REF] Holmquist | A computational constitutive model for concrete subjected to large strains, high strain rates, and high pressures[END_REF] which expands the hydrostatic part with a p-δ equation of state (see Herrmann [START_REF] Herrmann | Constitutive equation for the dynamic compaction of ductile porous materials[END_REF]) and combines the mesomechanical hydrostatic behaviour of the aggregates and the cement.

Another possibility is to use a nonlinear elastic material law, as it is presented [START_REF] Schmidt-Hurtienne | Ein dreiaxiales Schädigungsmodell für Beton unter Einschluss des Dehnrateneffekts bei Hochgeschwindigkeitsbelastung[END_REF]. The bulk modulus is changed depending on hydrostatic strain.

T P I R C S U N A M D E T P E C C A from Schmidt-Hurtienne

Element-free Galerkin Method

The proposed simulation model uses discrete cracks to model the non-linear behaviour of concrete. The cracks are implemented using the element-free Galerkin method (EFG). In contrast to the finite element method, the interactions between the nodes are defined only by the distance between the nodes.

MLS Interpolation

Belytschko [START_REF] Belytschko | Element-free galerkin methods[END_REF] proposes the element-free Galerkin method (EFG), which approximates a field by using moving least-squares interpolation (MLS Interpolation). The following equation is used to approximate the displacement field:

u h (x) = n i=1 φ k i (x)u i = p T • a (1)
The shape function φ is built from monomial functions p. A linear twodimensional example for p is:

p(x) = (1 x y) T (2)
The vector a is calculated by the minimization of the interpolation error J i

J i (a, x) = n i=1 w i (x) • u i -p T i (x) • a 2 (3)
The weight function w i depends on the distance x -x i . This weight function can be written as:

w i (x) = w i (s) with s = x -x i h i (4)
A common spline function is used as the weight function:

w(s) = 1 -6s 2 + 8s 3 -3s 4 (5)
T P I R C S U N A M D E T P E C C A

ARTICLE IN PRESS

The shape functions φ k i can be calculated by using the derivation of equation (3). We define a matrix P with the coordinates of the nodes:

P =           p 1 (x 1 ) p 2 (x 1 ) . . . p m (x 1 ) p 1 (x 2 ) p 2 (x 2 ) . . . p m (x 2 ) . . . . . . . . . . . . p 1 (x n ) p 2 (x n ) . . . p m (x n )           (6)
and a matrix W with the weight functions of the nodes w i :

W(x) =           w 1 (x) 0 . . . 0 0 w 2 (x) . . . 0 . . . . . . . . . . . . 0 0 . . . w n (x)           (7)
By using the auxiliary matrices A and B:

A = P T W(x)P (8) B = P T W(x) (9)
The derivation of J(a, x) can be written as:

∂J i (a, x) ∂a = A • a -B • u i = 0 (10)
Then, the interpolation can calculated with:

u h (x) = p T (x) • A -1 • B • u i (11)
The size of the radius of influence h i (see equation ( 4))is a variable parameter, which defines the number of nodes in support. If a too big value is chosen for this parameter, a large number of nodes lie in a supported area and the difference between the shape functions of two neighbouring nodes is therefore small. This results in an ill-conditioned interpolation; the possibility of numerical instability is high.

For the interpolation, the auxiliary matrix A has to be inverted. If the number of nodes in a supported area is smaller than the dimension of the problem, matrix A becomes singular. A is not singular if the nodes span two independent directions (for the 2D case). Experience shows that the size of the

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
radius of influence h i should be chosen with approximately 4 to 10 nodes in the supported area.

Discretisation of the EFG method is the same as used for the finite element method with the changed shape functions in equation [START_REF] Mosler | Embedded crack vs. smeared crack models: a comparison of elementwise discontinuous crack path approaches with emphasis on mesh bias[END_REF]. Discretisation in the model presented is Lagrangian. There are several descriptions of the dimensionality to discretise a model in two dimensions. The EFG model described is used here with plain strain, plain stress and axial symmetric discretisations.

The procedure for EFG is the same as for finite elements (e.g. described by Bathe [START_REF] Bathe | Finite Element Procedures[END_REF]). A three-dimensional implementation shows that the calculation time increases rapidly due to the much bigger search effort.

The integration can be done with background rather than nodal integration. Background integration uses the well-known Gauss integration of cells that are developed by a quadtree structure (see Belytschko [START_REF] Belytschko | Element-free galerkin methods[END_REF]) or triangulation or are defined manually by the user.

Nodal integration is proposed by Beissel [START_REF] Beissel | Nodal Integration of the element-free Galerkin method[END_REF] and uses as the integration factor the portion of the volume (or area) associated with a node. The combination of nodal integration and explicit time integration sometimes results in problems with under-integration. To prevent this instability Beissel [START_REF] Beissel | Nodal Integration of the element-free Galerkin method[END_REF] proposes a stabilisation term. Another possibility to stabilise the nodal integration is to add stress points that are used additionally for the integration (see Belytschko [START_REF] Belytschko | Stability analysis of particle methods with corrected derivatives[END_REF]).

Rabczuk [START_REF] Rabczuk | Stable particle methods based on Lagrangian kernels[END_REF] compares both integration methods using a static solution as well as dynamic solutions. The best results in comparison with analytical results were obtained with background integration but this method is significantly more expensive than nodal integration.

Background integration is mostly used for EFG. This integration method will also be used here to prevent stability problems. In the case of an explicit time integration the computational times for background integration are in the implementation used not significantly higher than for nodal integration.

MLS Interpolation and Cracks

In this work, the cracks are implemented by cutting off the weight functions (and hence the shape functions) at the location of the crack (see figure 1, Belytschko [START_REF] Belytschko | Element-free galerkin methods[END_REF]). This method is widely used in combination with the EFG method and has the advantage of being relatively robust. The geometry of the cracks is stored here with the level set method (see Ventura [START_REF] Ventura | A vector level set method and new discontinuity approximations for crack growth by EFG[END_REF]).

The domain of the node is divided into two subdomains: subdomain B beyond

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
the crack and subdomain A on the side of the node. In subdomain A the spline function w is the same as before; in subdomain B w is set to 0. Ventura [START_REF] Ventura | A vector level set method and new discontinuity approximations for crack growth by EFG[END_REF] proposes different methods for dividing the domain at a crack tip (e.g. visibility criterion). Herein, the diffraction method is used.

The implementation of a new crack results in smaller support for the nodes near to the cracks. To avoid a singular matrix A (see section 2.1) the radius of influence is updated in the case of a new or a longer crack. It could also be seen that the level set method used works very well with a small number of cracks. When using large numbers of cracks the effort of calculating the points affected by a crack is expensive.

Explicit Time Integration

The calculation of highly dynamic loaded structures requires such small time increments that an explicit time integration is reasonable. An explicit time integration is stable if the time step size is smaller then a critical value.

Menouillard [START_REF] Menouillard | Efficient explicit time stepping for the extended finite element method (X-FEM)[END_REF] shows that the stable time step of X-FEM elements is of the same order as that of the corresponding element without X-FEM extension. He suggests a time step of 2/3 of the stable time step of a corresponding element without X-FEM extension.

Menouillard [START_REF] Menouillard | Efficient explicit time stepping for the extended finite element method (X-FEM)[END_REF] also shows that the critical time step of X-FEM elements depends on the position of the crack. The critical time step increases when the crack is near a node. Figure 3 shows the influence of the location of the crack on critical time step. The location of the crack is described by using the distance to the next node s divided by the distance between the nodes l.

Lu [START_REF] Lu | Element-free Galerkin method for wave propagation and dynamic fracture[END_REF] shows that for the EFG method the time step is up to 1.7 times longer than that obtained by the finite element method.

The critical time step for the EFG method with cracks is calculated here by using the eigenfrequencies of the model. The critical time step Δt can then be calculated with the stability condition following the Newmark scheme:

Δt ≤ Δt c = 2 ω max (12)
The resulting critical time step for the EFG method is shorter than for a calculation without cracks. The time step for a crack near a node correlates with the resulting critical time step for an X-FEM element shown by Menouillard (see figure 3).

Another problem when using discrete cracks with explicit time integration is the occurrence of node oscillations on the crack surfaces by the time the

T P I R C S U N A M D E T P E C C A
cracks develops. Global artificial viscosity can damp these oscillations. The parameters for artificial viscosity are similar to the parameters used for the finite element method.

Constitutive Law for Concrete

A combination of a discrete and a cohesive crack model is used to calculate the concrete. The cohesive crack model takes account of the non-linear behaviour of crack development. The reason for this non-linear behaviour is the lack of homogeneity at meso level (the dimension of the visible components of concrete). The cohesive part replaces the damage model, which is used by most authors to model the smoothing behaviour after cracking. Instead of a damage surface a failure surface is used for crack development. With this procedure, an elastic material model can be used for the concrete. In the case of highly dynamically loaded concrete, the elastic material law is extended with descriptions of the behaviour of the micro pores.

Table 1 shows the parts of the used constitutive law.

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS 3.1 Constitutive Model -Bulk Modulus
The behaviour of concrete under hydrostatic load is affected by the behaviour of the micro-pores. Following a linear elastic hydrostatic loading the micropores fail (plastic loading). The compression stiffness decreases. Further compaction of the micro-pores results in increased stiffness. This effect is the reason if the wave is steeping up to a shock front.

But in most cases of shock waves the concrete is still loaded by a shock front.

The amplitude of a shock wave decreases very fast due to enormous local energy dissipation through the destruction of the micro-pores. This destruction results in a granular material.

Here, a nonlinear elastic material law should be used. The bulk modulus is changed depending on hydrostatic strain. Schmidt-Hurtienne [START_REF] Schmidt-Hurtienne | Ein dreiaxiales Schädigungsmodell für Beton unter Einschluss des Dehnrateneffekts bei Hochgeschwindigkeitsbelastung[END_REF] presents a volume function that considers the variation of the stiffness as a function of hydrostatic strain v . This volume function is used in this instance to calculate the compression modulus:

K tot = Y ( v ) • K (13)
with:

K = E 3(1 -2ν) (14)
The Poisson's ratio ν is chosen to 0.2. The Young's modulus E is chosen depending on the model.

The shape of the Y-function (figure 5) is calculated by

Y =            1 -a v • 1 -e -| v |-e v,th ev • 1 + | v |-e v,th bv•ev 2 for v < -e v,th 1 for v ≥ -e v,th (15)
The following parameters are used:

e v,th = 0.008 (16) e v = 0.0223 a v = 0.7 b v = 3.5
At a certain point, all pores of the concrete are damaged and additional compaction of the micro pores is no longer possible. Ruppert [START_REF] Ruppert | Zur numerischen Simulation von hochdynamisch beanspruchten Betonstrukturen[END_REF] defines the density of the unloaded concrete without micro pores as ρ ref = 2.181 g/cm 3 .

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
This corresponds to a value of Y = 4.8, which is used as a limitation of the Y-function.

Equation [START_REF] Rabczuk | A meshfree method based on the local partition of unity for cohesive cracks[END_REF] represents only the loading part. If unloading occurs after destruction of the micro-pores the maximum stiffness due to the increased stiffness of the existing granular material is used. The unloading volume function in plastically loaded concrete is higher than the initial stiffness. The function for the unloading path can be written as:

Y unloading =        1 for max Y loading = 1 AND min Y loading = 1 Y unloading,plast for max Y loading = 1 AND min Y loading < 1 max Y loading
for max Y loading > 1 (17) Investigations (see chapter 4.5) show that the numerical model represents the material by using Y unloading,plast = 1.5 (see also figure 21).

The failure of concrete under high hydrostatic loads also has to be considered. If the micro-pores of the concrete are destroyed, the concrete cannot react to a tension load. Furthermore, the shear modulus is reduced. If these effects are not considered, the dissipation of energy is too small and the amplitude of the wave does not decrease rapidly enough.

The reduced shear modulus is implemented with damage evolution as proposed by Ruppert [START_REF] Ruppert | Simulation of concrete structures under highly dynamic loadings[END_REF]

D z = v v,max γ (18)
with the hydrostatic strain v and the hydrostatic strain capacity v,max = 0.09. The parameter γ defines the shape of the damage. A parameter study shows that a value of γ = 0.1 is suitable to represent the behaviour of the concrete.

After the crushing of the micro-pores the concrete cannot sustain a tension load. The concrete reacts like a granular material. Therefore, if crushed concrete is loaded by tension, the integration points and the integration cells are deleted. The mass of these cells is removed from the mass matrix due to the fact that the fragment is no longer connected with the remaining concrete. The energy of these cells is therefore also removed.

Failure Surface

A failure surface can be used for different objectives. Here, the failure surface proposed by Hsieh [START_REF] Hsieh | A plasticity-fracture model for concrete[END_REF] should effect the development of the cracks.

Hsieh [START_REF] Hsieh | A plasticity-fracture model for concrete[END_REF] has converted the failure surface of Ottoson [START_REF] Ottosen | A failure criterion for concrete[END_REF] into a form which

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
allows easier adaptation to experimental results. The failure surface can be written as:

aJ 2 + b J 2 + cσ 1 + dI 1 -1 = 0 (19)
The parameters a, b, c, d can be calibrated with the biaxial experiments of Kupfer [START_REF] Kupfer | Behavior of concrete under biaxial stresses[END_REF] and the triaxial experiments of Mills [START_REF] Mills | Compressive Strength of Plain concrete under multiaxial loading conditions[END_REF]. The following values are used in this approach for equation [START_REF] Herrmann | Constitutive equation for the dynamic compaction of ductile porous materials[END_REF]:

a = 2.0108; b = 0.9714; c = 9.1412; d = 0.2312 (20)
Because these parameters are developed by scaling with the compressive strength f c , the values are independent of the strength of the concrete. This failure surface will therefore be used in this investigation.

Strain Rate Effect

The strength of concrete depends on the strain rate. The strain rate effect is shown by experiments of different authors. An arrangement of several experiments is shown by Bischoff [START_REF] Bischoff | Compressive behaviour of concrete at high strain rates[END_REF] (figure 4).

The strain rate effect has different reasons:

• The strain rate effect with small strain rates is due to the inertia of the concrete around the crack. Bischoff [START_REF] Bischoff | Impact behavior of plain concrete loaded in uniaxial compression[END_REF] shows that the inertia effects linked to crack propagation also predominate for high strain rates. • Humidity in the micro-pores results in a deceleration of micro-crack development both under tension loading and under compression loading. Reinhardt [START_REF] Reinhardt | Joint investigation of concrete at high strain rates[END_REF] shows the influence of humidity on the tensile strength by increased strain rates by using SHB tests with wet and dry specimens. Under compression loading the humidity has a reinforcing effect. Under tension loading capillary forces in the micro-pores resist displacement. Therefore, the velocity of crack development is reduced. • Another reason for the increase in strength is the homogenisation of the micro-structure. Static loading results in stress peaks at the edges of the aggregates. These stress peaks are followed by local failure of the structure.

If the structure is loaded in the range of the wave velocity, however, the stress distribution is smoothed and thus the structure can sustain a higher load.

Experimental data are used in the CEB-bulletin 187 [START_REF] Ceb | Concrete Structures under Impact and Impulsive Loading -Synthesis Report[END_REF] to develop a bilinear function for the increasing of strength. Extrapolation of the bilinear function in a range without experimental data (over approx. 200sec -1 ) results in physical unrealistic values. Gebbeken [START_REF] Gebbeken | A new material model for concrete in high-dynamic hydrocode simulations[END_REF] introduces a hyperbolic function to regard the limitation of the increased strength in this range depending on the damage level (see figure 4).

T P I R C S U N A M D E T P E C C A

ARTICLE IN PRESS

Hypothetically the maximum experimental value is set as limit. The dynamic strength is calculated by:

f incr = ˙ ˙ 0 δ for ˙ ≤ 30 sec -1 f incr = η ˙ 1 3 for ˙ > 30 sec -1 with ˙ 0 = 3 • 10 -6 sec -1 , δ = 1 10+ 1 2 f c , log η = 7δ -0.492 (21)
f dyn f stat = min(f incr , 8.0) (22)
The increasing strength is implemented by an expansion of the failure surface (for tension as far as for compression). Therefore, the strain rate has an influence of the initiation of the cracks and the development of the cracks. The increased strength is also regarded by the calculation of the additional forces in the fracture process zone.

Crack Development

Crack development proceeds in the following steps:

• Crack initialization:The basis for this decision is a failure surface. A Rankine criterion is sufficient if the concrete is loaded only by tension stresses, but it cannot represent the behaviour of concrete under triaxial loading. Therefore, in the proposed simulation model the failure surface of Hsieh, Ting, Chen [START_REF] Hsieh | A plasticity-fracture model for concrete[END_REF] is used (see section 3.2). • Decision whether a crack can growth: The same criterion as for the initialization is used. • Determination of the length of crack development: the simulation model uses a constant crack development length. The investigation in section 4.1 shows that a length of half the distance between the nodes is a good choice. • Determination of the direction of the crack: The direction orthogonal to the principle stress is used. This assumption is only true for MOD1 cracks. The examples (see chapter 4) show that the crack direction of MOD2 cracks is with this assumption represented precisely enough.

Knowledge of the stresses at the crack tip is necessary for the crack development criteria. The higher order continuity of the EFG method used obtains smoothed stress distributions at the crack tip. Therefore, the solution at the crack tip is similar to a non-local approach.

Dynamic loads sometimes cause a crack to close, making a contact formulation necessary. In this case the weight function has to be set to continuous again.

T P I R C S U N A M D E T P E C C A

ARTICLE IN PRESS

The closed crack can then transmit all loads between the crack sides.

In order to find a solution for MLS interpolation there are certain limits for the cracks. Cracks near the edges, cracks parallel to the edges, parallel cracks with a low distance and single nodes cut by a crack result in a divergence of the solution. Crack geometries of this kind should be avoided in such algorithms for crack development.

Maximum Crack Velocity

The velocity of crack development in concrete is lower than the Rayleigh wave velocity. Experimental results of different authors are shown in figure 6. The results show a limit in the crack velocity in a range of 500 m/s. This is the same limit like postulated by Curbach [START_REF] Curbach | Crack Velocity in Concrete[END_REF]. This limit is also validated by new experiments by Steiner [START_REF] Steiner | Beton unter Kontaktdetonation -Experimentelle Untersuchungen[END_REF]. These experiments with concrete loaded by contact detonations show a crack velocity in a range of 300 m/s.

The numerical investigations presented use a crack velocity capacity of v crack,max = 500 m/s. By using a constant crack development length l crack (section 3.4), crack development is possible after a time period of:

t crack = l crack v crack,max (23)

Fracture Process Zone

The process of crack development in concrete can be shown with a tension experiment (figure 2a). Shrinkage and creep generate micro-cracks in the cement paste. These micro-cracks grow when an external tension load is applied. Reaching 60% of the tensile strength, these micro-cracks are placed more and more in the direction normal to the loading (A). On reaching the tensile strength (B), most micro-cracks have joined. The tensile forces are only transmitted by a certain number of aggregates. After exceeding the tensile strength (C), a macro-crack has developed. The rough surfaces of the crack sides can still transmit forces between the crack surfaces (softening).

To describe this behaviour, Hillerborg [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF] proposes a "'fictitious crack model"', which hypothesises that the material behaviour is elastic until the strength is reached. After exceeding the failure limit, the strain increases in a local zone.

In contrast with brittle failure, this zone can still transmit stresses between the crack sides.

The stress-strain relationship can be split into two parts: the part until the

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
failure limit is reached and the part after exceeding the tensile strength. Because the position of the cracks is not known, the fracture process zone is not initialised before the crack is developed. Instead, the proposed simulation model uses an elastic material model (with a variable bulk modulus, see section 3.1) until the failure limit is reached (rigid cohesive model, see figure 2b). After cracking, additional forces are applied which represent the transmitted stresses of the fracture process zone (FPZ).

These so-called crack closure stresses represent the forces transmitted between the crack sides and are implemented as external loads at the nodes. A bilinear shape of the crack closure stresses is used:

σ ccs = (β t -1) • f t α t • w max w + f t for w < α t • w max ( 24
)
σ ccs = - β t • f t w max (1 -α t ) w + β t • f t 1 -α t for w > α t • w max
The parameter α t is proposed by Roelfstra [START_REF] Roelfstra | Numerical method to link strain softening with failure of concrete[END_REF] in the following form:

β t = 0.25 for f c < 30 MPa (25) β t = 0.25 -0.0015 • (f c [M P a] -30) for f c > 30 MPa
The coefficient α t is set to 0.14 as also proposed from Roelfstra.

Then, the crack opening, which corresponds to a stress-free surface (w max ), can be calculated with:

w max = 2 • G f f t (α t + β t ) ( 26
)
This form of the crack closure stresses provides a good representation of the experimental data.

Experiments with a Split Hopkinson Bar by Schuler [START_REF] Schuler | Spall experiments for the measurement of the tensile strength and fracture energy of concrete at high strain rates[END_REF] show that the values of the fracture energy depend on the strain rate. The fracture energy is duplicated at a strain rate of 100 sec -1 . This increase corresponds approximately to the increase in the strength at this strain rate. Experiments by Weerheijm [START_REF] Weerheijm | Concrete under Impact Tensile Loading and Lateral Compression[END_REF] do not show a dependency of fracture energy on strain rate.

Since detailed values are not available for the strain rate dependency of the fracture energy, the following assumption is used: the values of the strength in equation ( 24) are multiplied by the strain rate factor (see equation ( 21)); the values for the fracture energy are set as constant, whereas the maximum crack opening is calculated with the increased strength using equation [START_REF] Ventura | A vector level set method and new discontinuity approximations for crack growth by EFG[END_REF].

T P I R C S U N A M D E T P E C C A

ARTICLE IN PRESS 4 Validation and Examples

The simulation model is validated by the following examples.

Static Three-Point Bending Test

The static routine of the simulation model can be verified with a static threepoint bending test. An unreinforced concrete beam with an artificial crack on the bottom side of the specimen was tested by Körmeling [START_REF] Körmeling | Determination of the fracture energy of normal concrete and epoxy modified concrete[END_REF] to identify fracture energy. The geometry of this experiment is shown in figure 7.

The following material parameters are used: Young's modulus 20000 N/mm 2 , Poisson's ratio 0.20, tensile strength 2.4 N/mm 2 , fracture energy 113 N/m. The thickness of the specimen is 100 mm.

Numerical discretisation (plain strain) is achieved with three different meshes. Mesh1 uses a uniform mesh with a size of 5 mm. Mesh2 and mesh3 use a basis mesh of 10 mm. These meshes are refined in the region of the crack to 2.5 mm (mesh2) and 2 mm (mesh3). The meshes therefore use different numbers of integration cells through the height. While mesh1 is built with 25 elements through the height (26 nodes), mesh2 and mesh3 are built with, respectively, 40 and 50 integration cells through the height (41 resp. 51 nodes). By using the EFG method the refinement can be performed easily by cutting the integration cells. Triangular intermediate cells are not necessary.

The numerical investigations with bilinear softening in the fracture process zone show that there is a dependency on the mesh size (see figure 8). The coarser mesh shows a much higher load but a smaller softening part. The finer meshes represent the experiments much better.

Crack development is performed stepwise with a fixed crack development length (see chapter 3.4). The crack development length influences not only the development of the crack; it also has an influence on the fracture process zone. A shorter crack development length results in a better description of the cohesive behaviour due to the frequent updating of the crack. External forces resulting from the fracture process zone are also bigger. Therefore, the load capacity for a given displacement should be higher.

Due to this fact, the crack development length should be as small as possible to describe the crack geometry and the fracture process zone as accurately as possible. The crack development length should also consider other magnitudes used for the calculation or related to the concrete material. One magnitude used in the calculations is the distance between nodes and integration points.

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
The shape functions are an approximation between the nodes. Because of this, a much shorter crack development length than the distance between the integration points is not suggested.

Concrete is a composite material made up of a cement paste and aggregates. Usually, the maximum size of the aggregates is in a dimension of 2 cm. Due to the fact that the presented approach uses a homogenized model for the concrete, a crack development length much smaller than the maximum size of the aggregates is not sensible. The investigations also show (figure 9) that the crack development length has a non-negligible influence on the numerical result. A length of half of the integration cell is suggested and used here. This choice leads to a good representation of the experimental results.

The load-displacement curve is not continuous. This is due to the fact that the discrete cracks develop with a constant crack development length, resulting in steps in the load-displacement curve. This corresponds with the failure of the aggregates in the experiment.

Shear Test

The development of curved cracks can be tested using the shear experiment of Arrea [START_REF] Arrea | Mixed-Mode Crack Propagation in Mortar and Concrete[END_REF]. The geometry and the material parameters are shown in figure 10.

The specimen has a thickness of 156 mm. The mesh is build with a distance of the nodes of 17 mm with an refinement near the proposed crack to 9 mm. The crack development size is chosen to 4.5 mm. The numerical crack pattern resulting by a calculation with a plain stress model compares well with the experimental pattern (figure 11). The crack patterns resulting from the use of a damage model (e.g. Akkermann [START_REF] Akkermann | Rotationsverhalten von Stahlbeton-Rahmenecken[END_REF]) or other discrete crack models (e.g. Prasad [START_REF] Prasad | Fracture mechanics model for analysis of plain and reinforced high-performance concrete beams[END_REF]) are steeper than the results with the proposed simulation model.

Dynamic Three-Point Bending Test

Experiments run by Du [START_REF] Du | Fracture process zone for concrete for dynamic loading[END_REF] show the behaviour of dynamically loaded cracks. Unreinforced concrete beams (geometrical and material data see figure 12, thickness 50.8 mm) are loaded with a falling weight (9.68 kg, altitude of fall 1 m). The beam has an implemented crack on the underside of the specimen. Strains are measured at three locations.

The results of these experiments can be used to validate the simulation model for dynamic crack development. The calculations are done with a plain strain discretisation. The mesh has a size of 6 mm; the crack development size is 3 mm. The numerical results of dynamic crack development show oscillations

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
around the crack. The reason is discontinuous crack development. These oscillations can be reduced by using artificial viscosity.

The numerical results reflect the experiments with the measuring sensor on the bottom (SG1). The strains of the numerical calculation at the measuring sensor SG2 are somewhat delayed, but representation of the results is still good. The crack edges at the measuring sensor SG3 in the calculation are unloaded after the experiment and therefore the strain in the simulation is too large.

Shearing of a Plate under Detonation Load

Plates loaded by heavy air blast waves fail through shearing rather than bending. Experiments conducted by Albritton [START_REF] Albritton | Response of deep two-way reinforced and unreinforced concrete slabs to static and dynamic loading[END_REF] with reinforced concrete show the crack pattern illustrated in figure 14.

To simulate these experiments a specimen is loaded by a triangular load-time function with a maximum load of 20 MPa and a rising time of 0.05 msec. A background mesh with a distance between the nodes of 10 mm is used (plain strain). The strain rate effect is considered.

The following material parameters are used: Young's modulus 34500 N/mm 2 , Poisson's ratio 0.20, fracture energy 141.75 N/m, static tensile strength 4.14 N/mm 2 . The crack development length is chosen to 3 mm considering the distance between the nodes and the size of the aggregates. The maximum crack velocity is chosen to 500 m/sec.

Figure 15 shows the numerical crack pattern. The crack runs from the support into the specimen. Additional cracks develop on the right-hand side of the support due to the disregard of reinforcing in the numerical simulation.

Contact Detonation

The aim of this study is to calculate concrete loaded by contact detonation. The Institute for Reinforced Concrete Structures and Building Materials -Universität Karlsruhe (TH) has carried out experiments to identify the parameters of highly dynamic loaded concrete. The test specimens have a size of 1 m x 1 m x 0.5 m. The region of interest in the concrete body has no reinforcement. The specimen was loaded by a plane wave generator with a mass of 640 g TNT. This loading introduces a plane wave across a large charged area.

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
Loading with this kind of charge does not destroy the test specimen. Different failure ranges can be identified in figure 16. The detonation results in a crater beneath the explosive. The concrete underneath the crater is highly compacted. Below this region concrete is damaged by cracks. The concrete between these cracks is not visibly damaged.

The costs of the element-free Galerkin method are much higher than those of mesh-free particle methods (e.g. Rabczuk [START_REF] Rabczuk | A three-dimensional large deformation meshfree method for arbitrary evolving cracks[END_REF]) due to the search effort, especially with cracks. Therefore, three-dimensional numerical modelling of the test specimen with the EFG method was not possible with reasonable effort. Therefore, the EFG method is used with an axial symmetrical approach (figure 18).

The convergence of the mesh is investigated with a smaller model (geometrical data see figure 17). A comparison of this model with experimental values is not possible. The results show (see also figure 17) that a distance from the node of 3 mm results in similar pressures to finer meshes. Therefore, the distance between the nodes in the region beneath the explosive is chosen for the full model as 3 mm and in the other regions as 6 mm. The model uses 10786 degrees of freedom.

The model is loaded with a pressure-time function applied to the nodes beneath the explosive according to the results of Plotzitza [START_REF] Plotzitza | Ein Verfahren zur numerischen Simulation von Betonstrukturen beim Abbruch durch Sprengen[END_REF], see figure 19. The boundary conditions are presented in figure 18. The following material parameters are used: Young's modulus 30000 N/mm 2 , Poisson's ratio 0.2, static tensile strength 3.6 N/mm 2 , maximum crack velocity 500 m/sec, crack development length 3 mm. The cracks are built by using the failure surface of Hsieh (see section 3.2). Cells and nodes, which are loaded over the pressure capacity and then loaded with hydrostatic tension, are deleted (see also section 3.1). This loading results in fragments which will not be considered here.

The variation in crack direction can be determined in a section through the specimen (figure 20). This section shows a maximum crack direction variation of 50 o , which is used as the limit for the calculations.

The experiments of Steiner [START_REF] Steiner | Beton unter Kontaktdetonation -Experimentelle Untersuchungen[END_REF], Herrmann [START_REF] Herrmann | Experiments on concrete under shock loading[END_REF] and Ockert [START_REF] Ockert | Ein Stoffgesetz für die Schockwellenausbreitung in Beton[END_REF] can be used to compare maximum pressures as a function of the distance from the top of the specimen. The experiments show fast decreasing pressure.

Figure 21 shows the maximum pressures using different unloading functions (parameter Y unloading,plast , equation 17). The numerical results with a value of Y unloading,plast = 1.5 represent the experimental results.

Rabczuk [START_REF] Rabczuk | A three-dimensional large deformation meshfree method for arbitrary evolving cracks[END_REF] also investigated these experiments using the cracking particle method (see Rabczuk [START_REF] Rabczuk | Cracking particles: a simplified meshfree method for arbitrary evolving cracks[END_REF]). The resulting maximum pressures as well as the resulting crack pattern are similar. The difference with this method lies in

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
the fact that the EFG method results in a discrete crack pattern whereas the cracking particle method results in regions with cracked particles.

The numerical and experimental failure patterns show a region with deleted nodes beneath the explosive, representing the crater (figure 16). The amplitude of the pressure wave decreases very fast. The crack pattern resulting from the detonation load is affected by transverse tensile cracks. After reflection of the wave on the underside of the specimen, the cracks rotate and a spalling region develops. The experimental crack pattern compares well with the numerical crack pattern (figure 16, main cracks are highlighted).

Conclusion

This work presents a new approach which uses a discrete crack model (instead of a damage model) to calculate shock wave-loaded concrete. Discrete cracks are modelled using the element-free Galerkin method. The discrete cracks are combined with a strain rate-sensitive fracture process zone. The behaviour of the micro pores is considered by using a new description of the Y function which concerns also the maximum possible compression of the pores and the unloading part.

The results presented show that crack geometry for static and dynamic loading of concrete structures can be represented with the proposed model. The maximum pressures and the crack geometry resulting from loading concrete with contact detonation can be calculated using the proposed simulation model.

The simulation model presented here uses only two-dimensional approaches.

The extension of the element-free Galerkin method to a three-dimensional approach is possible. The computational costs increase when using threedimensional calculations due to the fact that the search effort is much more expensive. The most complex enlargement is the development of the crack in three dimensions. The crack line becomes a surface; the crack tip becomes a crack tip line. Gravouil [START_REF] Gravouil | Non-planar 3d crack growth by the extended finite element and level sets -part 2: Level set update[END_REF] shows that the level set method can also be used for the development of three-dimensional cracks. Also the fracture process zone has to be enlarged to the third dimension. To simplify the implementation of the three-dimensional fracture process zone, a specialised transparency or diffraction criterion for the crack tip could be used, which considers also the cohesive behaviour of the crack.

The advantage of the simulation method presented is the fact that the cracks are available as cracks and not as a zone of damaged elements. However, the calculation time seems to be much longer than a comparable model which using damage material.

Fig. 1 .Fig. 3 .Fig. 4 .Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .Fig. 10 .Fig. 11 .Fig. 12 .Fig. 13 .Fig. 14 .Fig. 15 .Fig. 16 .Fig. 17 .Fig. 18 .Fig. 19 .Fig. 20 .Fig. 21 .

 1346789101112131415161718192021 Fig. 1. Weighting functions a) Without crack b) Subdomains in the case of a crack c) Weighting function with crack

Acknowledgement

The support of the Institute for Reinforced Concrete Structures and Building Materials, Universität Karlsruhe (TH), Germany, is gratefully acknowledged.