Brian Patrick Kernan

Andrew Lynch

Dr Con Sheahan

Dr Dang Tung

A Novel Metric for Determining the Constraining Effect of Resources in Manufacturing Simulation

Keywords: SIMULATION, BENCHMARKING manufacturing, enterprise, resource, machine, skill, starvation

This paper provides a novel method for determining the constraining effect of resources in a manufacturing system using discrete event simulation. Traditionally manufacturing systems are constrained by one or more bottlenecks. Eliminating or mitigating the bottleneck will speed up the system throughput. However, bottlenecking resources generally only refer to machines, and primarily focus on flow-shops not job-shops. One important resource we believe that is often overlooked is workers and their associated skills, and we propose that a particular skill could be flagged as a bottleneck resource. We define new metrics known as Resource Constraint Metrics (RCM) for measuring the constraining effect of a resource on the entire manufacturing system. These metrics are flexible and differentiate between the constraining effects of machines and their requested skills. The metrics can also deal with complex workflows with alternative routing, alternative resources, calendars (a necessary consideration when dealing with workers), worker performance, and multiple modes of operation of machines (e.g. run, setup, and maintenance). The use of RCMs in simulation aids in realworld decision-making, by determining which resource should be focussed on and improved to reduce the overall system constrainedness. This will have the effect of increasing throughput or at least providing the capacity for increased throughput.

Introduction

The work presented in this paper is founded primarily on the area of research of bottleneck detection in manufacturing simulation. There are numerous definitions as to what constitutes a bottleneck [START_REF] Lawrence | Shifting Production Bottlenecks: Causes, Cures, and Conundrums[END_REF]. In Roser's paper [START_REF] Roser | Shifting Bottleneck Detection[END_REF] a bottleneck is defined as a stage in a production system that has the largest effect on slowing down or stopping the entire system. Finding the bottleneck is no trivial task, and Cox and Spencer [START_REF] Cox | The Constraints Management Handbook[END_REF] for example simply recommend that '…the best approach is often to go to the production floor and ask knowledgeable employees.' There are however a number of systematic methods available to locate the bottleneck for production systems however. One approach is to measure the utilisation of the different machines of the production system (Law 2000). The machine with the highest utilisation is considered to be the bottleneck. However, the utilisations of different machines are often very similar, so relatively long simulation run-times may be required to determine the system bottleneck with a high degree of confidence.

Another frequently used method analyses the queue lengths of the machines in the production systems. In this method, either the queue length or the waiting time is determined, and the machine with the longest queue length or waiting time is considered to be the bottleneck. This method of detection is similar to looking at machine utilisation, since a machine with a large queue should also have commensurately high utilisation.

The shifting bottleneck method [START_REF] Roser | Shifting Bottleneck Detection[END_REF] further distinguishes between a momentary bottleneck, describing the bottleneck at any given point in time, and an average bottleneck, describing the bottleneck behaviour over a selected period of time. The shifting bottleneck approach is based on the theory of constraints [START_REF] Blackstone | Theory of Constraints -A Status Report[END_REF] (Goldratt 1992). Although most manufacturing systems usually have one main bottleneck, in all but the simplest applications, bottlenecks are not static, but rather shift between different machines. These shifts may for example be due to the sequence of random events or due to a gradual change in the manufacturing system. A non-bottleneck machine may become a bottleneck, for example due to a machine failure, and similarly a bottleneck machine may become a non-bottleneck machine. Over longer periods of time, a system therefore may not only have one primary bottleneck, but also secondary and tertiary bottlenecks, i.e. machines which are also occasional bottlenecks, yet to a lesser extent that the primary bottleneck. This highlights the fact that it may be possible to improve throughput by not just focusing on the primary bottleneck. Roser [START_REF] Roser | Shifting Bottleneck Detection[END_REF] deals with shifting bottleneck detection and uses an example with a job-shop arrangement with alternative routing, and alternative resources.

In Muthiah et al's paper [START_REF] Muthiah | Overall throughput effectiveness (OTE) metric for factory-level performance monitoring and bottleneck detection[END_REF]) the authors define a metric known as Overall Throughput Effectiveness (OTE) which builds upon the well established Overall Equipment Effectiveness (OEE), which is widely accepted as a primary performance metrics in manufacturing [START_REF] Hansen | Overall Equipment Effectiveness -A Powerful Production/Maintenance Tool for Increased Profits[END_REF]. While OTE has the benefits of monitoring total factory performance, being a useful aid in machine bottleneck detection in a varied selection of workflow sub-systems (e.g. machines in parallel, multiple routing) workers and skills were not considered during the development of the metric. [START_REF] Li | Data Driven Bottleneck Detection of Manufacturing Systems[END_REF]) use analytical methods for bottleneck detection, based on historical data on machine blockages and starvations and also take into account machine failure times. The authors also note that using analytical methods for bottleneck detection are generally only appropriate for relatively simple workflows comprised of machines in sequence. However in Li's other paper (Li 2009a) an analytical data-driven method is used for bottleneck detection in a complex manufacturing system. Again however, workers are not dealt with in these two of Li et al's papers. For larger more complex systems, that include worker resources simulation would need to be used.

All of the aforementioned approaches for bottleneck detection generally assume heavily automated flow-shop production times where the human element can be ignored. However, in a large proportion of enterprises, workers are an important resource and cannot be disregarded when looking at bottleneck detection. The reason they tend to be ignored in simulation, is that they add an extra layer of complexity and they also necessitate the modelling of calendars, shifts, and breaks. Once workers are included in bottleneck analysis, it muddies the water and the conventional analyses may not suffice. Consider a simple production line with four machines and a single worker. Each machine has the same processing time and requires the worker to operate it. In this case the bottleneck resource is not any of the machines, but is in fact the worker. If we were to add a second worker, we will double the throughput.

Another important issue which is not dealt with are machine modes of operation, such [START_REF] Lin | Manufacturing Process Analysis with Support of Workflow Modelling and Simulation[END_REF]. Our workflows can deal with alternative routings such as OR-Splits and AND-Splits and the corresponding Joins, and also deals with alternative ranked resources. In addition the workflows can also contain subcontracting activities. On the resource level, machines can have multiple modes of operation such as setup and preventive maintenance, sequence dependent setup times, and workers can possess multiple skills, with an efficiency level being defined for each worker skill combination.

For the rest of this paper, in section 2 we define the metrics which we use to measure the contributing effect each resource has on the overall system constrainedness, and in section 3 we present some computational examples for calculating RCM values. Section 4 covers the determination of RCM values through a simulation of a real-life production plant.

Resource Constraint Metrics

Modelling Workers, Skills, and Efficiencies

Firstly we define a machine as an internal locus where a manufacturing activity is performed. A machine can process orders as well as preventive maintenance requests (PMRs), and emergency maintenance requests (EMRs). A skill meanwhile, is a worker ability that is requested by a machine in its particular mode of operation. For example a lathe may require a skill 'Can operate Lathe' in its run mode of operation, and the skill 'Can repair Lathe' in its emergency maintenance mode of operation. A worker can have a set of skills, and each skill for a particular worker can have an efficiency between 0% and 100%. A skill can be regarded as an abstract resource, and a worker can only provide one of their skills at a time. The efficiency increases the processing time for a machine as follows: If the nominal processing time for a machine is 10 minutes, and if a worker with an efficiency of 75% is assigned to the machine, then the actual processing time will be 13.3 minutes. We must also define the concept of skill starvation. We say a machine is neglectfully starved for skill when capable workers are busy elsewhere on the shop-floor, and therefore unavailable to satisfy a request from the machine.

Orders and PMRs may also be in queue for a machine that is suffering from absenteeism induced skill starvation, which means that there are no workers on the shop-floor for the current shift capable of operating the machine in either its standard modes of operation (run, setup, load, unload, teardown) or its relevant preventive maintenance mode. Our definition of starvation is not to be confused with the concept of a machine being starved for parts. In traditional simulation analysis, a machine that is idle and starved for parts is often frowned upon, and it is usually pointed out that its utilisation should be increased. However, since our simulator, models the entire enterprise, all parts/orders can be traced back to external customer demand. A machine that is idle is not necessarily a negative thing; if there is no associated customer demand there may be no point providing the machine with extra parts to produce (the exception being if the enterprise is building to stock). However, a machine that is starved for skill is impeding the movement of parts through the system, parts which will generally have actual associated customer demand.

Formulae for Resource Constraint Metrics

For each machine and skill the simulation model will calculate a value known as the Resource Constraint Metric (RCM) to quantify how much of a limiting factor that resource is to the manufacturing operations, and thus help identify it as a resource bottleneck. The higher the RCM value, the more of a limiting factor that resource is.

Shown below are the actual formulae for calculating the RCMs for machines and skills. for a more detailed explanation.

Machine

2 EMRs will never be in queue for a machine. When a machine fails, an EMR is generated and instantly replaces the order that was on the machine at the time of failure.

Explanation of Resource Constraint Metric Formulae

Excess processing time is incurred due to worker inefficiency. This will only be relevant for machines that have workers assigned to them, and those workers have efficiency levels of less than 100% for a particular skill. For example, consider the nominal run-time for a machine is 10 minutes. If a worker with 80% efficiency for the requested skill by the machine is operating the machine, then the actual run-time will be 12.5 minutes. Other orders in queue will be waiting and thus contribute to the RCM, but for each order we add only 10 minutes for the RCM for the machine and add the other 2.5 minutes to the RCM for the skill. The 2.5 minutes is the excess processing.

We also separate out the different RCMs based on the shift when they were recorded. Separating out RCMs based on shift is useful in determining not only where the problem skill is, but when it is. For example, the main contributing factor to a high RCM for a particular skill could be due to a higher absenteeism rate during the night shift than the day shift. Sample RCMs are given in (Figure 1). We can see that we have a higher RCM value for skill 2 during the night shift than the day shift.

Therefore, one possible recommendation to reduce the SCM is to increase the skill-set of a worker who works the night shift, by providing them with skill 2.

It should be pointed out that the equation on the previous page for calculating skill RCMs is applicable to multiple skill requests. For example a machine is neglectfully starved for skill for 10 minutes and is requesting 2 workers with skill A and the only two workers who possess this skill are on their lunch-break. If the machine has been assigned an order, and there are no other orders or PMRs in queue for the machine, then a value of 20 minutes will be appended to the RCM for that machine. If one worker had been idle and the other on their lunch-break, then the recorded value would be only 10 minutes. Another thing to note is that even though we model machine alternatives, an order can only ever be in queue for a single machine at any instant in time. When an order moves onto its next activity, and there are multiple machine alternatives, the order is moved to the queue which will result in it being processed fastest3 . A simple analogy is a person in a supermarket; they look at their options in terms of check-outs, and they will choose the check-out which they believe will result in them leaving the supermarket quickest. It is possible for an order to switch to a queue for a different machine however. For example if a machine fails suddenly, the orders that were in queue for it may be moved over to a queue for an alternative machine, if the time to repair the first machine is expected to be longer than setting up the alternative machine.

The calculations do not attempt to weight orders based on priority, quantity, or cost. We do not weight by priority, since higher priority orders should spend proportionally less time in queue for machines, and on machines that are starved for skill than lower priority orders anyway. Also, we do not weight by quantity or cost, since the calculations also include maintenance requests; if we were to try and weight an order, we were have to apply some arbitrary weighting variable to maintenance requests. However, future work may dictate that it is necessary to weight each order by the total euro-cost associated with each order. In this way a high cost order will contribute proportionally more to the RCM for the constrained resource, than lower cost orders. The weighting variable for maintenance requests would be a euro-cost that captures the delay of a PMR in queue for a machine, or a starved machine that is processing a PMR or EMR. Assigning a weighting variable to maintenance requests is in the realm of fuzzy logic however; intuitively we can say that the RCM from a machine that fails unexpectedly and is waiting for a worker rated to repair it, should be worse that the RCM for a machine that is starved for skill while processing a low priority work order.

When we run a simulation each machine and skill will have a generated RCM value, which will generally have units of hours. We are only interested in relative RCM values, and the resource with the highest RCM value will generally be the one we should focus on if we want to reduce the system impedance. At this stage we also define the concept of the System/Shop-Floor Constraint Metric (SCM) which is simply the sum of all the individual RCMs. We may find that if we reduce the RCM of one resource (e.g. reducing the processing time of a machine) then the RCM of another resource may increase. In this situation we could fine that the SCM remains the same, indicating that there is no overall net improvement of the system. This will be the case if in a simple production line we reduce the processing time of a nonbottleneck machine.

Advantages of Resource Constraint Metrics

The advantage of using RCMs over simply trying to locate the bottleneck machine, is that in any manufacturing systems which are other than the most simple cases, we need a way of ranking the contribution of each resource to the system constrainedness. The case of performing analysis on a manufacturing system and claiming that 'machine X is the bottleneck' will generally only be tenable for a single sequential production line, with no workers modelled. Roser's paper [START_REF] Roser | Shifting Bottleneck Detection[END_REF] has expanded the research by creating a methodology to identifying secondary and tertiary bottlenecks, but the RCM is a more flexible way of quantifying (relatively)

the contribution of each machine to the overall system constrainedness. Once we factor in workers and skills, which are often neglected in simulation, but a critical resource in most enterprises, we need a way to differentiate between the contributing effect of machines and workers. The RCMs do this in a intuitive manner, and can be used to aid in real-world decision-making. By looking at the status quo or as-is situation on the manufacturing shop-floor, the RCM values will inform the user where best to focus their efforts in order to reduce the overall system constrainedness e.g.

add a new machine, reduce the processing time of an existing machine, add a new worker, or upskill an existing worker. This is particularly relevant for small enterprises with relatively small workforces, and where workers must be included as a resource in the simulation of the system; upskilling an existing worker could potentially have an equivalent rate-of-return as purchasing a new machine but without the capital investment. Subsequent simulations with the changes in place should support the RCM data i.e. we should see a reduction in SCM, and we may see a reduction in the number of late satisfied customer demands (a positive consequence of orders flowing through the system quicker). We may also see an increase in the finished goods inventory, since a decrease in SCM may result in reduced work order fulfilment lead time. It is important to note however, that a high RCM for a skill associated with a particular machine can be tackled by reducing the processing time for the machine. This is because by reducing the processing time, there is less opportunity and time for the machine to be starved for skill.

Another advantage of using these metrics is that the coding for the calculations is incorporated directly into a database driven general simulation model that has been developed in eM-Plant. Our approach instantiates a model of an enterprise from a RDB which results in significantly reduced build effort. Since the coding for calculating RCM is build into the simulation model in eM-Plant the coding for the metrics are not specific to a single case study. Multiple instances of enterprises can be represented in relational database format, and the RCM values calculated, once the model is automatically instantiated in eM-Plant. The common approach to building simulation models is to build them ad hoc; from the ground up, using components and elements in a simulation application user interface. For example the paper by Alfieri [START_REF] Alfieri | Workload simulation and optimisation in multi-criteria hybrid flowshop scheduling: a case study[END_REF]) describes a relatively complex simulation of a cardboard company which contains elements such as multi-machine stations, sequence-dependent setup times, resource calendars, and external operations. Similarly the paper by Nacsa [START_REF] Nacsa | SIMAN Simulation of Assembly Operations in a Furniture Company[END_REF]) demonstrates a model that includes components such as machines, shifts, and workers. These models are however each specific to the shop-floor of a single manufacturing firm. In the case of the paper by Nacsa, the model was built using Siman, a simulation application.

Our simulator is capable of instantiating models from a database at least as complex as described in [START_REF] Alfieri | Workload simulation and optimisation in multi-criteria hybrid flowshop scheduling: a case study[END_REF]) and (Nasca 2008) and will produce RCMs, and other decision oriented KPIs. The RCMs could therefore potentially be used for controlling resources on an operational level, in a similar manner to the paper by Li et al. [START_REF] Li | Real time production improvement through bottleneck control[END_REF] where the authors have achieved real time production improvement through bottleneck control.

Computational Examples

Example One

This section will describe three computational examples. The first example deals with a simple sequential workflow with four machines each with a constant processing time of one minute. We start the simulation with 200 orders in queue for the first machine, and we end the simulation when all the orders have left the system.

At the end of the simulation we find the following RCM values as shown in (Table 1).

As expected, the only non-zero RCM value for the machines is for the first machine. Since the first machine feeds machine B at a rate of one order per minute, and machines B, C, D have a processing time of one minute, no orders will be in queue downstream of machine A. Also, in this simple example we are not dealing with workers, hence none of the machines will undergo skill starvation. The reason machine A has a non-zero RCM value is because at the start of the simulation there are 200 orders in queue for the machine. From the formula for RCM for machines, the RCM value for a machine is proportional to the queue size or WIP (work-inprocess) for the machine; the more parts in queue for a machine, the higher the RCM for the machine, and also for the skill associated with the machine if the machine undergoes skill starvation. The RCM for example one can be calculated by noting that at time 1, there will be 199 orders that have been in queue for machine 1 for 1 minute. Hence we add 199 minutes to the RCM for machine 1. At time 2, there will be 198 orders that have been in queue for machine 1 for 1 minute, so we add on 198 minutes to the RCM value. Therefore, the total RCM value at the end of the simulation should be given by: Which for a value of n = 199 yields an RCM value of 19,900. The reason for the discrepancy between the calculated value of 19,900 and the simulated value of whereby events which are technically supposed to occur at the same time they were trigged, must be scheduled for 1 second in the future. Events must be dealt with in this manner in eM-Plant in order to avoid events being inadvertently bypassed by the simulator. However, the error in RCM values for this simple model is only approximately 1.6%, and since the simulator is aimed at simulating enterprises over relatively long periods of time, this incurred error is negligible, and should not negatively affect the output results in a decision making context. Next, if we reduce the processing time of machine A to 30 seconds (New Scenario 1), we obtain the following results as shown in (Table 2).

Because machine A can process orders twice as fast as before, a queue builds up in front of machine B. Hence machine B will now have a non-zero RCM. As expected however, the SCM stays the same at 20,225 minutes. This indicates that halving the processing time of machine B has had no effect on the overall system constraint metric. Only if we reduce the processing time of all four machines will we notice a reduction in the SCM. Shown in (Table 3) are the results when the processing time of all the machines is 30 seconds (New Scenario 2). The value of 10,278 minutes is approximately half of 20,225 minutes; we have reduced the SCM by 50% which corresponds to a doubling of throughput.

Example Two

In the second example, we use the model from first example with some alternations.

The processing times of all machines are now stochastic, uniformly distributed between 55 and 60 seconds, and orders are now released to the system at a rate of one order per minute. This will have the effect of reducing the queue in front of the first machine. Since the processing times are stochastic we run the simulation for 10 replications until all the orders have left the system. Shown in (Figure 2) are the RCM values obtained for this example.

For each machine, from the X-axis to the thick black line in the middle of the grey box, the corresponding Y value is the mean RCM. From this value, to the top and bottom of the grey box, represents one standard deviation4 above and below the mean. We see that machine A has a zero RCM value. Since orders are now released at a rate of one per minute, and the longest processing time capable for the machine is one minute, there will never be orders in queue for machine A. The RCM values increase the further downstream we go. This is due to the cumulative effect of the variability of the processing times. Next if we change the processing time of machine C to a uniform distribution of between 80 and 85 seconds we obtain the following results for the RCMs, for 10 replications as shown in (Figure 3).

We see that machine C dominates. Again this is to be expected since it is now the bottleneck because it has the longest processing time of the four machines, and thus will have the largest queue in front of it. Machine B has a non-zero RCM but it is dwarfed by the RCM for machine C. The RCM value for machine D is now zero, since machine C will release orders at a rate of less than one per minute, and since the processing time of machine D will never exceed one minute, there will never by any orders in queue for machine D. If we want to make improvements to the system, we should concentrate on machine C. Also it is worth noting that the RCM value for machine D in (Figure 3) is 130 hours, whereas when all the machines had the same processing time, the highest RCM value for machine D is around 0.015 hours (54 seconds) (Figure 2). There is a significant difference in scale.

Example Three

In the third example we introduce workers and skills, job dependent processing times, and multiple workflows. We have three sequential workflows in total

• Workflow 1 = {Activity A → Activity B → Activity C → Activity D} • Workflow 2 = {Activity B → Activity A → Activity D → Activity C} • Workflow 3 = {Activity C → Activity B → Activity A → Activity D}
Activities are performed at machines as follows:

• Activity A performed at {Machine A} We assume machines only have a single mode of operation "Run" and machines request the following skills for their "Run" mode of operation. Also each machine sorts jobs in the queue by their priority, and then by shortest processing time. 4). The information on jobs is shown in (Table 5). Each job is given the same arbitrary scheduled release date, which corresponds to time zero in the simulation model, and the same internal due date. The priority of an order dictates its queuing priority and worker dispatching priority i.e. a high priority order will be given higher preference to be processed ahead of a lower priority order.

An overview of the relationships between workflows and skills is given for workflow 1 in (Figure 4). The linkage type between activities is 'sequence', and simply defines a precedence constraint dependency between activities i.e. we can not start activity B until activity A is completed. Looking at this diagram, we see that Con will need to share his time between machines C and D. Therefore there is a potential for these two machines to be starved for skill.

We run the simulation for a simulated time of 45 minutes, after which all the jobs have been processed. Shown in (Figure 5) is a Gantt chart for the machines. We use a different colour to denote each job, and when a job is on a machine that is starved for skill, we use a thin binary colour strip with the same colour as the job plus the colour red. This Gantt chart does not try to represent jobs which are queued for machines. For example, job 5 (yellow) between time 20min and 25min is in queue for machine D, but during this time period machine D is busy processing job 4 (grey); job 5 is not shown in the Gantt chart for this period. Also, machine C starts processing job 4 ahead of job 6 because job 4 is higher priority than job 6. Machine A processes job 1, job 5, and job 3 in that order because they have increasing processing times. If we look at machine C at time 14 we see that job 5 is loaded onto the machine, but the only worker capable of operating machine C is Con and he is currently busy with job 2 on machine D. Therefore machine C is neglectfully starved for skill C for two minutes. When Con finishes job 2 he moves up to machine C and begins processing job 5. Also, as expected at any instant in time, machines C and D will never be simultaneously processing. The worker Andrew only possess skill A, and since machine A is the only machine that requests this skill, Andrew is confined to machine A. This implies that it is impossible, under the current conditions, for machine A to be starved for skill. It is the same situation with the worker Brian and machine B. Shown in (Figure 6) is a Gantt chart for the jobs. We use a different colour to denote each machine, and when a machine is starved for skill, we use a thin binary colour strip with the same colour as the machine plus the colour red.

We will now show how we calculate the RCM values for machine D and skill D which are values of 17 minutes and 10 minutes respectively. In (Table 6) the second column Q D contains, for each job, the total time the job spent in queue for machine D. The third and fourth columns break this time down, for each job, into the time when machine D was processing, and when machine D was starved for skill (skill D in this case). For example, looking at (Figure 6) for job 6 we see that it spent a total of 11 minutes in queue for machine D. However, for 3 of these 11 minutes machine D was starved for skill D when processing job 5. So this 11 minutes is separated out to 8 minute for Q(p) D and 3 minutes for Q(p) DD . In the fifth and final column, we record the time for each job, when the job was on machine D and the machine was starved for skill D. When we sum up the third column in (Table 6) we get a value of 17 minutes which is the RCM value for machine D. When we sum up the fourth and fifth columns we get 3 minutes and 7 minutes respectively, and when we add these together we get a value of 10 minutes which is the RCM for skill D. We do not have to consider the other machines in the calculation of the RCM for skill D since they will never be starved for this skill.

dShown in (Figure 7) are the RCMs for all the machines and skills. Looking at this graph, some recommendations could be giving skill C or skill D to the workers Andrew or Brian, or reducing the processing time of machine B. These changes could reduce the SCM and also reduce the makespan of the jobs i.e. the time when the last job leaves the system. However, it should be pointed out that our model here is relatively simple, and over a very short period of time. Decisions made based on RCM values should be for real world systems that are in steady state and over relatively long periods of time. In this way we can be confident in the relativity of RCM values and the dominance of the constraining power of one resource over another on the shop-floor.

Real Life Example

In this section we use the simulation model to determine the RCM values for an entire factory. The factory under examination is a SchuF Fetterolf production plant in Ireland that designs and manufactures industrial valves that control, isolate, divert, or sample gases, liquids, slurries, or powders. The main characteristics of the plant being simulated are as follows:

• Simulation time of one year This information was stored in the companies ERP (Enterprise Requirements Planning) system and was translated and stored in the Relational Database for the simulation model. The primary source of variability being modelled is worker absenteeism and tardiness. Based on the following papers [START_REF] Barmby | Worker Absenteeism: An Analysis Using Microdata[END_REF], (Bytestart 2008), (Irish Times 2008) we assume the following parameters:

• Probability of being absent is 3.3%

• Probability of being late is 4.4% (However, 65% probability that when a worker is late this will occur on a Monday)

• Average Lateness is 10 minutes, with an upper limit of 18 minutes. Hence, assume actual lateness time is log-normally distributed with mean 10 minutes and standard deviation 2 minutes, where probability of being more than 18 minutes late is less than 5%.

The main caveat with modelling absenteeism and tardiness in this manner is that we assume all workers have an equal probability of missing a day of work, or being late for work. In reality of course, some workers will be more prone to missing work than others. A more accurate way to model this variability is to look at historical data for each worker and assign a probability for each worker. However, this data was not currently available at the time this paper was written. It should be noted that with a probability of being absent of 3.3% and with 24 workers the chances of at least one worker being absent on a given day is greater than 55%.

Another source of variability in the model is that equally ranked alternative resources are randomly selected. For example, if a machine requires a worker with the skill 'Can Setup' and two equally ranked workers possess this skill, then one of the workers will be selected at random. The same rule applies to alternative machines for a particular activity. We run the simulation for ten replications, once with worker This means that when we eliminate worker absenteeism and tardiness this translates to a 79% reduction in skill RCM.

When we perform a similar statistical analysis on the total machine RCM we find that there is no evidence of a significance difference in means of total machine RCM between the two experiments, and this can be seen in Figure 9. This result was to be expected since eliminating worker absenteeism and tardiness should not significantly reduce the machine RCMs. If instead the experiments were only looking at the effect of machine break-downs, which can be regarded as the machine equivalent of worker absenteeism, then we should find a significant difference between machine RCMs and a possible significant difference between skill RCMs. Machine RCMs and skill RCMs in this context would be dependent on the mean time to repair the machine (MTTR), and the availability of skills to repair the machine.

When we look at the System Constraint Metric (SCM) which is simply the sum of skill RCMs and machine RCMs we see that machines are the main contributor to SCM (See Figure 10.). The final result is that based on our simulation model of the SchuF plant, when we eliminate worker absenteeism and tardiness the SCM reduces Analysis of the different skill RCMs could also inform us which skills are critical, and simulation could also inform us the effect of absenteeism from different workers on the system constrainedness.

Summary

The metrics presented in the paper are an intuitive and highly flexible way of determining the constraining power of a resource in a manufacturing system, via simulation. They are superior to traditional bottle neck detection methods by taking into account, many aspects of a manufacturing that are usually ignored in simulation models, such as workers, calendars, multiple modes of machine operation, and complex workflows. Resource constrains metrics are a way of quantifying the constraining power of worker-skills as well as machines, and can aid in real-world decision making, by determining which resource should be focussed on and improved, in order to reduce the overall system constrainedness. This will have the effect of increasing throughput or at least providing the capacity for increased throughput. Additionally it can result in improved customer demand satisfaction, since reduced system impedance can mean a commensurate reduction in the number of incidents where a customer's sales order is satisfied late. Further work using these metrics and one that the authors are currently working on is using these metrics to identify the most constrained skill, and simulating worker trainer, and how this improves the economic performance of the organisation e.

  handling, and maintenance. Each mode may require a different set of skills to operate e.g. a technician may be required to perform maintenance on the machine, or a fitter and an electrician to setup the machine to process a new product type.In our research, we have developed a database driven simulation model that is applicable to a generic manufacturing enterprise. All the necessary information about the enterprise such as resources, orders, workflows, and work calendars, is stored in a relational database (RDB), and a simulation engine, eM-Plant, instantiates a model of the enterprise from the database. The workflows we have developed are based on work by Russell et al. (Russell 2004) and Lin et al.

  or Average Efficiency of Team of Workers

  is due to a condition known as event time shifting in the simulation model,

  in queue for M D Q(p) D = Time job in queue for M D when M D is processing Q(s) DD = Time job in queue for M D when M D is starved for S D H DD = Time when M D was neglectfully starved for S D

  average of 4.1%. Further work would involve analysing historical data and obtaining absenteeism and tardiness percentages based on each worker.

Figure 1 .Figure 3 .Figure 4 .Figure 5 .Figure 6 .Figure 7 .

 134567 Figure 1. Sample Resource Constraint Metrics (By Shift) Figure 2. Resource Constraint Metrics for Computational Example 2 (Baseline) Figure 3. RCMs for Computational Example 2 (New Scenario 1) Figure 4. Breakdown of Workflow 1 for Computational Example 3 Figure 5. Machine Gantt Chart for Computational Example 3 Figure 6. Job Gantt Chart for Computational Example 3 Figure 7. Resource Constraint Metrics for Computational Example 3 Figure 8. Skill RCMs for Production Plant Figure 9. Machine RCMs for Production Plant Figure 10. SCMs for Production Plant

  of data to see if there is sufficient evidence that the data deviates from the normal model. We get p-values greater than 5% meaning we can assume the data is normally distributed. Next we check the data sets to see if we can assume equal variances. In this case we get a p-value less than 5% indicating that the variances are unequal. This result was to be expected since the lions-share of variance is from modelling absenteeism and this is only present in the first experiment. The only source of variance present in the second experiment is from the random selection of equally ranked alternative resources. Finally we perform a Welch test for testing for significance difference between means with unequal variances and we get a p-value less than 5% indicating that the means are not equal. This is evident from Figure8.

	absenteeism and worker tardiness being modelled, and for another ten replications
	with absenteeism and tardiness switched off.
	We test for a statistically significant difference (with a significance level of 5%)
	between the mean total skill RCM for the two different experiments using R, a
	statistical analysis software package. Firstly we perform a Kolmogorov-Smirnov test
	on both sets
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Table 5 .

 5 g. upskilling worker Andrew with Skill B with yield the best economic performance compared with all other worker-skill training combinations. Russell, N., Van Der Aalst, W. (2004). "Workflow Resource Patterns: Identification, Representation and Tool Support", Eindhoven University of Technology. Jobs for Computational Example 3

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Table 6 .

 6 Breakdown of Queue/Starvation Time for Computational Example 3

					Replication Number
		1	2	3	4	5	6	7	8	9	10
	Skill RCM (Hours)	64	61	37	25	52	39	12	15	15	30
	Machine RCM (Hours) 592 592 591 598 591 592 592 598 588 583
	SCM (Hours)	656 653 628 623 643 631 604 613 603 613

Table 7 .

 7 Simulation Results for Production Plant with Absenteeism and Tardiness

					Replication Number
		1	2	3	4	5	6	7	8	9	10
	Skill RCM (Hours)	5	6	8	10	6	7	7	11	7	7
	Machine RCM (Hours) 594 595 596 588 600 594 593 587 595 595
	SCM (Hours)	599 601 604 599 606 602 601 598 603 602

Table 8 .

 8 Simulation Results for Production Plant with no Absenteeism and Tardiness

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.ukInternational Journal of Production Research

This is a prediction within the simulation at run-time however, and may actually turn out to be wrong. For example if the selected machine fails unexpectedly, the alternative machine may have been a better choice.

We use sample standard deviation, with Bessel's correction for measuring the spread of the results, which we store in a separate RDB. We do not use confidence intervals, since confidence intervals assume that the outputted results follow the normal model distribution which may not be the case.

Machine

Resource Constraint Metric / Minutes Machine A


Prof. Ryszard Gawinecki

Ryszard Gawinecki

Doctor; Piotr Borowski

Anna Miłaczewska

Doctor; Agnieszka Skotnicka

Professor; Krzysztof Woliński

Agnieszka Brzyska

acetic acid. Tautomerization versus decarboxylation

Keywords: Carboxylic acids, decarboxylation, tautomerization, reaction path, transition states, DFT calculations Carboxylic acids  Decarboxylation  Tautomerization  Reaction path  Transition states  DFT calculations

The DFT calculations at the B3LYP level with 6-311G** basis set were carried out in order to reveal whether tautomerization or decarboxylation is responsible for the instability of 2,2-di(pyridin-2-yl)acetic (DPA) and 1,8-diazafluorene-9-carboxylic (DAF) acids. The carboxyl protons in both compounds are involved in the intramolecular hydrogen bonds (the pyridine nitrogen atoms are the hydrogen bond acceptors). Although formation of two intramolecular OH•••N hydrogen bonds in the enols of both carboxylic acids enables effective electron delocalization within the quasi rings (•••HO C=C C=N), only ene-1,1-diol of DAF has somewhat lower energy than DAF itself (ΔE is ca 7 kcal/mol). DPA and its enediol have comparable energies. Migration of the methine proton towards the carbonyl oxygen atom (to form enediols) requires overstepping the energy barriers of 55-57 kcal/mol for both DPA and DAF. The enaminone tautomers of the acids, formed by migration of this proton towards the pyridine nitrogen atom, are thermodynamically somewhat more stable than the respective enediols. The energy barriers of these processes are equal to ca 44 and 62 kcal/mol for DPA and DAF, respectively. Thus, such tautomerization of the acids is not likely to proceed. On the other hand, the distinct energetic effects (ca 15 kcal/mol) favour decarboxylation. This process involves formation of (E)-2-(pyridin-2(1H)-ylidenemethyl)pyridine and its cyclic analogue followed by their tautomerization to (dipyridin-2-yl)methane and 1,8-diazafluorene, respectively. Although the later compound was found to be somewhat thermodynamically more stable, kinetic control of tautomerization of the former is more distinct.

Therefore, I think the paper would benefit if the authors elaborate more on the specific (medicinal) properties of these compounds (in the Introduction section), and draw more general conclusions about the chemical behavior of this class of molecules by contrasting the findings of the paper with those obtained for other similar systems (in the Conclusion section). One minor thing: misprint on p.3 (bottom): Fig. 11 should be Fig. 5.

2,2-Di(pyridin-2-yl)acetic acid and its cyclic analogue are supposed to be unstable compounds. We found DFT procedure to be very helpful to show what is the reason of this instability. Susceptibility of these compounds to undergo tautomeric rearrangements and decarboxylation are the only properties discussed in our paper (no other features of unstable compounds can be studied: they decompose during the measurement). On the other hand, explanation what are the products the labile compounds are transformed to, would be very helpful to other chemists trying to use them. Despite their instability, 2,2-di(pyridin-2-yl)acetic acids may appear as the intermediates in numerous (bio)chemical processes. Conclusions of the present paper may also refer to other substituted acetic acids (and their homologues) that carry electronegative heteroatoms in the gamma positions, such as diacylacetic acids and their derivatives (to be also known as unstable compounds). "Fig. 11" at the former p. 3 (bottom) is now changed to "Fig. 5". New fragments added in the text are in yellow.

Reviewer #2: This is a well written manuscript reporting theoretical study undertaken with the hope to clarify the reasons why the title compound is unstable. My major concern is about the very high barriers for TS1 and TS4, both of which are connected with transition states featuring 4-membered rings. The presence of such strained rings in TS is probably a major reason for the high activation barriers. If the hydrogen atom transfers could take place between two molecules of species 1 or 4, then the strained rings would be avoided in the transition states.

For example, b3lyp/6-311G**//b3lyp/6-31G calculations for 4a dimer shows that the activation energy (no ZPE included) for transformation into 5a dimer is only 17.9 kcal/mol, which is much lower than the reported value of 42.4 kcal/mol. Thus, I strongly recommend the Authors consider bimolecular reactions for steps 1 -> 2 and 4 -> 5 in the revised manuscript. below are the coordinates (in Angstroms) for structures optimized by the Reviewer: 4a dimer E b3lyp/6-311G** = -1069.61126123806 N 1. 5489423635 We would like to thank the reviewer for his/her valuable comments. We admit, we were not considering the bimolecular mechanisms from the very beginning. This was definitely the point worth further investigations. The reaction like (4a)2 → 2 5a involves the proton transfer between the heteroatom (nitrogen, in this case) and the carbon atom. Indeed, we found the stable (4a)2 dimer at the B3LYP level, just as the reviewer. The counterpoise (CP) corrected value of the binding energy turned out to be somewhat above 4 kcal/mol. We also found the transition state for the N1(A)…H9(A)…C7(B) proton transfer (A and B denote the different molecules of a dimer). The energy barrier, being close to 20 kcal/mol, was lower than that found for the unimolecular mechanism. Similar calculations were also carried out in the case of b-type systems (DAF), and the conclusions were basically the same as in the case of DPA; all data are now reported in the revised manuscript. However, in the meantime, we recalled the calculations we carried out before on systems interacting via the dispersive-type forces, the forces that may play a crucial role in the case of two 4a molecules at some of their relative orientations. Note, that in the case of N(A)H(A)…C(B) system the electrostatic interactions are quite weak. In such cases DFT (as also known from the literature) is not sufficient for the accurate determination of the molecular properties. We therefore decided to carry out the MP2 calculations, which are well suited for such problems. They are really time-consuming -the PQS quantum chemistry package does not allow the parallel MP2 calculations (this is actually the reason for the delay in answering the reports sent to us). MP2 results did not confirm the DFT predictions. The "MP2 dimer" we found turned out to be stabilized by two NH(A)…N(B) hydrogen bonds. Such bonds enable the chemical exchange of protons between the heteroatoms rather than the proton transfer from N to C. We are inclined to trust the MP2 calculations in this particular case; the bimolecular mechanism seems to be the artifact of the DFT calculations. The relevant comments regarding the bimolecular mechanism are now in the revised manuscript. On the other hand, the remaining DFT results (i.e. these, regarding the unimolecular processes) are expected to be accurate. All processes investigated take place within the single, covalently bonded systems, for which the DFT calculations provide very accurate results.

Finally, we would like to emphasize that we are grateful to the reviewer for this comment, as well as for his/her effort to carry out the initial DFT calculations. The additional calculations we carried out were very educative to us. We are very much grateful for the referee for its/their valuable suggestions. All of us (the authors) feel that the manuscript sounds much better now.

Introduction

Ene-1,1-diols, often called enols of carboxylic acids, are intermediates in hydration of ketenes [1,2]: R 2 C=C=O + H 2 O  R 2 C=C(OH) 2  R 2 CHCO 2 H. These compounds are less stable then the corresponding carboxylic acids being energetically favored by delocalization of the lone electron pair of the hydroxy oxygen atom (Fig. 1) [1]. It is known, however, that ene-1,1diols that contain bulky aromatic groups in the molecule, are more stable [1][2][3][4]. The hindered protonation of the  carbon atom in Ar 2 C  =C(OH) 2 caused by ortho alkyl groups (Ar = 2,4,6trimethylphenyl, 2,4,6-triisopropylphenyl or 2,3,4,5,6-pentamethylphenyl) is responsible for such a behavior. Enediol of cyclopentadiene-1-carboxylic acid (Fig. 2) was found to be relatively stable [5]. Two intramolecular hydrogen bonds of the RAHB (Resonance-Assisted Hydrogen Bond) type [6][7][8][9] are expected to be present in its derivative annulated with two pyridine rings (Fig. 3). Unfortunately, neither 1,8-diazafluorene-9-carboxylic (9Hcyclopenta[1,2-b;3,4-b']dipyridine-9-carboxylic) acid nor its enol are known, but close analogue of the former compound, 2,2-di(pyridin-2-yl)acetic acid, was claimed to be obtained earlier from 2-diazo-1,2-di(pyridin-2-yl)ethanone [10]. Since numerous heterocyclic acetic acids are potentially interesting as anti-arthritic agents [11,12], the compounds mentioned above seem worthy to be studied from point of view of their stability. <Figure 1> <Figure 2> <Figure 3> Loosing of the CO 2 molecule by 2-and 4-pyridylacetic acids proceeds smoothly even below 100ºC (this process is much more difficult for 3-pyridylacetic acid) [11][12][START_REF] Katritzky | Handbook of Heterocyclic Chemistry[END_REF][START_REF] Doering W Von | [END_REF][15][16][17]. 2-(Pyridin-2-and 4-yl)phenylacetic acids are also exceptionally unstable: benzylpyridines were found to be the only products of the acid hydrolysis of 2-(pyridin-2-and 4-yl)phenylacetonitriles [18]. On the other hand, the enediol and enaminone tautomeric forms of the closely related 2,2-di(pyridin-2-yl)acetic acid (Fig. 4) can be stabilized by two intramolecular hydrogen bonds. Its tendency to decarboxylate is not reported in the only paper devoted to synthesis and properties of this compound [10].

<Figure 4>

Decarboxylation of carboxylic acid involves formation of carbon dioxide and an organic residue [16]. Organic product is usually stabilized by delocalization of its unshared electron pair. Formation of the intermediate zwitterion HOCOC 5 H 4 N  Ө OCOC 5 H 4 N  H is a common step in the mechanism of enzyme-catalyzed decarboxylations of 2-and 4pyridylacetic acids [12,[START_REF] Doering W Von | [END_REF]17]. Electron delocalization is particularly important in enzymatic reactions, where an "electron sink" is generally provided by a coenzyme [16].

Despite their instability, the labile compounds may appear as intermediates in numerous (bio)chemical processes. Explanation what are the products these derivatives are transformed to, seems very important. Numerous attempts done in our laboratory to obtain 2,2-di(pyridin-2-yl)acetic acid (1a, Fig. 5) from 2-diazo-1,2-di(pyridin-2-yl)ethanone according to the known procedure [10] as well as by hydrolysis of 2,2-di(pyridin-2-yl)acetonitrile (the standard procedure) were unsuccessful. Irrespective of the reason of inaccessibility of this compound to us, its properties can still be discovered by quantum chemical calculations. In the present paper 2,2-di(pyridin-2-yl)acetic (1a) and 1,8-diazafluorene-9-carboxylic (1b) acids are considered to be susceptible to tautomerization and decarboxylation. The related pyridylacetic acids were earlier used as the models when studying decarboxylation mechanism of aminoacids [19]. It is noteworthy that numerous heterocyclic derivatives of acetic acid are potentially interesting as anti-arthritic agents [20].

Computational details

Density functional theory (DFT, see e.g. [START_REF] Parr | Density-Functional Theory of Atoms and Molecules[END_REF]) has been recognized for years to be capable of predicting the very accurate conformations of the covalently bonded systems. Thus, its use in the present paper is justifiable. All calculations were carried out using B3LYP hybrid functional [22,23] with the 6-311G** basis set [24].

Formulas of 2,2-di(pyridin-2-yl)acetic (1a) and 1,8-diazafluorene-9-carboxylic (1b) acids, their enol (2a,b) and enaminone (3a,b) tautomeric forms as well as the (expected) intermediates (4a,b) and final products of the decarboxylation process (5a,b) are depicted in Fig. 5, which also shows numbering of heavy atoms. It should be noted that the carboxyl hydrogen atom H9 (in 1, 2 and 3) has the same number even if it is bound to N1 (in 4) or C7 (in 5). Formulas of the transition states that are (or may be) formed during decarboxylation of 1a and 1b can be seen in Fig. 6. The geometry optimization of all substrates, products and transition states was carried out first. Then, the harmonic frequencies were calculated to establish the types of stationary points, and to find the zero-point energy (ZPE) corrections for the energy barriers and the energetic effects of the processes investigated in the present work.

All frequencies for systems 1-5(a,b) were real. On the other hand, one imaginary frequency was found for transition states TS(1-4)a,b. Analysis of the eigenvectors corresponding to the imaginary frequencies provided us with the information about the possible products "on each side" of a transition state at a given stage of the reaction.

<Figure 5> <Figure 6>

Geometry optimization for the transition state is much more complex than for the stable system. This prompted us to carry out some initial, point-wise calculations, corresponding to the so-called distinguished reaction coordinate (DRC). This procedure is not recommended for the accurate description of the reaction path, for example on account of possible discontinuities of energy (and geometrical parameters) as a function of DRC (i.e. some selected geometrical parameter, e.g. bond length). It may, however, provide a reasonable initial structure for which typical algorithms are capable of converging quickly to the required transition state. Therefore, the C7C8 bond was selected as the DRC in reactions 1 → 4 + CO 2 .

Its length, R, was gradually increased from the equilibrium value found for 1 (the applied increment was 0.1Å), the remaining parameters were optimized, and the geometry corresponding to the highest energy was determined. Since H9 turned out to be still bound to the carboxylic oxygen atom, it was placed in the middle of the O9 ... N1 distance prior to running the transition state search. On the other hand, the two internal coordinates r 1 = H9N1 and r 2 = H9C7 (cf. Fig. 5) were selected for the description of the proton migration in the 4 → 5 tautomerization. Energies of structures optimized for the remaining geometrical parameters were then computed at a grid of selected linearly independent points. Then, the cross-section of PES (see Fig. 7 for an example), approximated by the 3 rd degree polynomial, was subjected to straightforward search of the saddle point. All initial structures obtained in this way turned out to converge quickly to the corresponding transition states by following (maximizing along) the lowest (negative) Hessian eigenmode. In order to follow the reaction paths, concept of the so-called intrinsic reaction coordinate (IRC) [25] in mass-weighted Cartesian coordinates was used. All calculations were carried out with the parallel version of the PQS quantum chemistry package [26,[START_REF]PQS version 3.3, Parallel Quantum Solutions[END_REF]. <Figure 7>

Results and discussion

Molecular geometries

The optimized structures of carboxylic acids as well as these of the enol and enaminone tautomeric forms and products of their decarboxylation (dipyridyl-2-ylmethane and its 1,2dihydro tautomeric form) are depicted in Fig. 8. Their most important geometrical parameters, as well as geometrical parameters of the transition states TS(1-4)a,b are presented in Tables 1 and2. <Figure 8> <Table 1> < Table 2> The most significant conclusions regarding the molecular geometries are as follow:

1. Arrangement of the N1 and N1' atoms in the most stable conformer of 2,2-di(pyridin-2yl)acetic acid 1a is trans-like relative to the C2C7C2' link (Fig. 8). Turning of one pyridine ring in 1a by ca. 90° around the C7C2' bond increases the energy by only 0.3 kcal mol -1 . Linking of the pyridine rings via the C3C3' bond in 1b prevents their free rotation and makes the molecule rigid. Note that the N1C2C7C2' and N1'C2'C7C2 torsion angles in 1b are close to 180°, i.e. the pyridine rings in the molecule are nearly coplanar.

The hydroxyl groups both in 1a and 1b are involved in the intramolecular hydrogen bond with the N1 atom being the hydrogen bond acceptor. The OH bond in 1a is somewhat longer than 0.99 Å (a typical OH bond length in the dimers of ordinary carboxylic acids where very strong intermolecular hydrogen bonds are present). On the other hand, the OH bond in 1b is shorter (ca. 0.98 Å). Rigidity of the 1b molecule and, consequently, steeper energy increase along the interring torsional motions, prevent the H9 and N1 atoms to approach each other as close as in 1a. This results in elongation (weakening) of the H•••N1 hydrogen bond (1.918 Å in 1b vs. 1.738 Å in 1a) and shortening (strengthening) of the O9H bond. The carbon-nitrogen bond lengths in the pyridine rings of 1a and 1b are close to these in pyridine itself (1.337 Å; geometry optimized at the B3LYP/6-311G** level).

The same applies to the ring carbon-carbon bond lengths which are comparable to the average bond length in pyridine (1.393 Å). Slight elongation of the carbon-nitrogen bonds in 1a and 1b is observed for the ring involved in the hydrogen bonding. The valence CNC angles in these systems are also comparable to those found in the molecule of pyridine (117.2°).

Formation of the H•••N1 and H•••N1

' hydrogen bonds in 2a and 2b affects the hybridization of C7 and shortens the C7C8 bond lengths (1.419 Å in 2a and 1.376 Å in 2b) as compared to these in 1a and 1b. These bonds resemble more the carbon-carbon bonds in the benzene ring (~1.39 Å) rather than those in alkenes (~1.33 Å). It should be noted that C7C8 bond in 2a is longer, while this in 2b is shorter than 1.39 Å. Effective electron delocalization in the quasi rings that involve OH•••N fragments results also in shortening of the C7C2 and C7C2' bonds both in 2a and 2b (they are still significantly longer than the aromatic carbon-carbon bonds) as compared to these in 1a and 1b (the C7C8, C7C2 and C7C2' bonds in their molecules are typical carbon-carbon single bonds).

The H9O9 and the hydrogen H•••N1 (') bond lengths in 2a,b follow the same pattern as these in the molecules of 1a,b. Moreover, C3C3' bond in 2b is much longer than C7C2 ( ' ) and C2 ( ' ) C3 ( ' ) bonds (which have comparable lengths, cf. Table 1). Thus, aromaticity of the five membered ring in 2b is doubtful what negates stabilization of this form by resonance suggested in Fig. 3.

The NH•••O hydrogen bonds are present in the enaminone tautomeric forms 3a and 3b

(these bonds are longer than the OH•••N hydrogen bonds in 2a and 2b). Contribution of the zwitterionic structures of 3 (cf. Fig. 5) is manifested in shortening of the C7C8 bonds (this effect is not as significant as for 2) and elongation of the C8O10 bonds as compared to 1.

4. 4a and 4b seem likely to be the intermediate products of decarboxylation of 1a and 1b, respectively. Both these molecules are planar. Hydrogen atom involved in the H•••N1 hydrogen bonding in 1a and 1b is bound to N1. The N1H bond lengths (1.009 Å in both systems) are in between these found for pyrrole (1.006 Å) and aliphatic amines (1.015 Å; geometry optimized at the B3LYP/6-311G** level). The C7C2 and C7C2' bond lengths are intermediate between these of the single and double carbon-carbon bonds. In both 4a and 4b the C7C2 bonds are shorter, while C7C2' are longer than 1.39 Å (Table 1).

5. 5a and 5b are probably the final products of decarboxylation of 1a and 1b. Their geometrical parameters remain similar to these found for 1a and 1b, as well as those in pyridine itself. Note that the pyridine rings in 5b are coplanar.

Reaction paths and energetic effects

Calculations of the energetic effects of the most intuitive reactions 1 → 5 + CO 2 reveal that products are by ca. 13.6 (case a) and ca. 15.8 kcal mol -1 (case b) more stable than substrates.

All energetic effects |E products -E substrates | and energy barriers E TS -E substrates/products discussed (denoted as ΔE) were corrected for the zero-point energy, ZPE. They are depicted in Fig. 9.

The data used in order to obtain these values (i.e. total energies of the molecules and ZPE corrections of all stable structures and transition states) are reported in Table 3. Note that for the brevity reason the ΔG values were calculated for decarboxylation processes only.

Comparable entropies of the substrates and products for the tautomerization processes suggest that the  ΔE and ΔG values also approximate each other.

<Figure 9> <Table 3>

Ene-1,1-diols 2 have lower energies than the respective acids 1: the energy lowering is equal to approximately 0.7 and 7 kcal mol -1 for 2a and 2b, respectively (Fig. 9). The energy barriers for the proton migration from C7 towards the carbonyl oxygen atom O10 (this processes proceed via TS1 transition states) are nearly the same for a and b and amount to ca. 55 kcal mol -1 . The enaminone tautomeric forms 3a,b are thermodynamically somewhat more stable than 2a,b (by nearly 4 kcal mol -1 in the case of 3b). However, the energy barriers for the migration of proton H7 towards N1', associated with the presence of TS2a and TS2b transition states, differ significantly: the calculated ΔE values are equal to ca. 44 and ca. 62 kcal mol -1 for 1a → 3a and 1b → 3b, respectively. In a view of finding more energetically favorable processes from the thermodynamic and kinetic point of view (vide infra) we conclude that tautomerization of 1a and 1b is not likely to proceed. Among all investigated processes, namely 1 → 2, 1 → 3, and 1 → 5 + CO 2 , the most distinct energetic effect can be seen for 1 → 5 + CO 2 (for both a and b, cf. Fig. 9). These processes will be discussed in details.

Rupture of the C7C8 bond in the 1 → 4 + CO 2 reaction is probably followed by migration of the carboxylic hydrogen atom to C7 (this results in formation of 5). The concerted and stepwise mechanisms have to be considered here. In the first one, the eigenmode (atomic displacements) corresponding to the imaginary frequency should correspond to simultaneous detachment of CO 2 and migration of the carboxylic proton H9 towards C7. The later mechanism is expected to involve a few different transition states. The case a will be considered first. A quick glance at the structure of 1a points at the second, i.e. step-wise mechanism. The reason for this is a significant energy increase when H9 approaches C7, while C7 is still bound to C8. It is shown in Fig. 10 as a steep, violet (turning to blue at larger values of R) "slope" on a cross-section of the potential energy hypersurface of 1a. In addition, elongation of the C7C8 bond reveals that there is a clear trend for the H9 atom to approach N1, i.e. to form 4a. Therefore one should search for a transition state that appears in the course of reaction 1a → TS3a → 4a + CO 2 . It was found to be nearly 20 kcal mol -1 above 1a.

The calculated ΔG value corresponding to the reported ΔE amounts to 18.33 kcal mol -1 . This value is known to be low enough for the process that is fast under standard conditions. Analysis of geometrical parameters shows that C7C8 bond in TS3a is significantly longer than that in 1a, and that the H•••N1 bond remains close to that of 4a (cf. Tables 1 and2). The structure of TS3a is depicted in Fig. 11 (the eigenmode, corresponding to the imaginary frequency equal to i320cm -1 , is also shown there). Following the reaction path from TS3a along the intrinsic reaction coordinate (IRC), 1a or 4a + CO 2 were obtained. No energy lowering is observed at this stage of decarboxylation. However, value of the Gibbs free energy of the products calculated at 298 K is by more than 12 kcal mol -1 lower than that of 1a.

Although rotational and vibrational entropies were also changed (ΔS rot = 11.9 cal mol -1 K -1 and ΔS vib = -12.2 cal mol -1 K -1 ), significant increase in the translational entropy (ΔS trans = 36.6 cal mol -1 K -1 ) gives the major contribution to the reported value of ΔG. The low energy barrier being < 20 kcal mol -1 (corresponding to energy of quanta in the near IR region) as well as negative value of ΔG clearly show that 1a decomposes spontaneously (and rapidly) to form 4a and CO 2 . This proves that Eistert and Schade [10] could not obtain 1a. It is probably why we were not able to obtain this compound starting from 2-diazo-1,2-di(pyridin-2-yl)ethanone [10] as well as by hydrolysis of 2,2-di(pyridin-2-yl)acetonitrile (the standard procedure). <Figure 10> <Figure 11> 4a is probably an intermediate in the process of decarboxylation of 1a. In order to obtain 5a, which is the most stable system studied (at least at the B3LYP level, cf. Fig. 9 and Table 3), the C7H9 bond has to be formed and, at the same time, the N1H9 bond has to be broken. It was not clear to us what is the mechanism of this step (i.e. unimolecular or bimolecular). It was initially assumed that the rearrangement takes place within the single molecule. This process is accompanied by breaking of the planar symmetry of 4a (note that C7 atom has to change its hybridization from sp 2 in 4a to sp 3 in 5a).The search for another transition state (TS4a) was carried out in the way described in Section 2. Its structure along with the imaginary eigenmode (i1834cm -1 ) is shown in Fig. 11. Location of H9 more or less half way between N1 and C7 (cf. Table 2) is the characteristic feature of this transition state. The energy barrier ΔE = 42.37 kcal mol -1 that has to be overstepped, corresponds to ΔG value of ca. 38 kcal mol -1 . It is two times higher than that for the CO 2 detachment. This seems to be mostly due to the high strain within the 4-membered ring formed in the TS4a transition state.

Thus, this step of the decarboxylation process, though spontaneous from the thermodynamic point of view, is expected to be much slower than the former one.

An alternative mechanism involves two molecules of 4a. The B3LYP calculations show, that they are capable of forming a relatively stable dimer possessing the C 2 symmetry. It is depicted on Fig. 12a; the binding energy, counterpoise (CP) corrected for the basis set superposition error [28] and ZPE corrected, is equal to 4.26 kcal mol -1 . The transition state, corresponding to the imaginary frequency of i1281 cm -1 , was found to be only 14.58 kcal mol - 1 above the dimer (the ZPE corrected value; the corresponding G value is 11.41 kcal mol -1 ).

Analysis of the imaginary eigenmode shows that this is the transition state for the migration of the H9 atom from N1 atom of one 4a molecule to the C7 atom of the other 4a molecule (actually, two molecules of 5a were finally obtained by following the reaction path along the intrinsic reaction coordinate). The B3LYP energetic effect of the overall reaction, i.e. (4a) 2 → 2 5a, is 10.85 kcal mol -1 per molecule (G=17.14 kcal mol -1 ). However, it should be noted that the dimer predicted by B3LYP calculations is stabilized by the interactions of the NH fragment of one molecule (a rather weak dipole) with the non-polar C … C … C system of the other molecule (the quadrupole, cf. Fig. 12a). Thus, the dispersive-type interactions that show the U disp ~ r -6 dependence on distance, are not negligible as compared to the weak electrostatic dipole-quadrupole interactions U dip-quad ~ r -4 . The same type of interactions are expected in the transition state. Since DFT is known to poorly describe the dispersive interactions, we decided to investigate this problem using the MP2 method with the same basis set. The geometry optimization of the dimer of 4a was carried out at the MP2/6-311G** level, using the B3LYP structure as an initial guess. This time we did not obtain the analogous dimer; two NH•••N hydrogen bonds were formed instead. This complex, for which the CP corrected binding energy is 14.29 kcal mol -1 (the MP2 frequencies were not calculated for practical reasons, but the ZPE estimate obtained at the DFT level, would further reduce this value by ca. 0.7 kcal mol -1 ), is shown in Fig. 12b. Thus it seems, that the DFT calculations overestimate the dispersive interactions when two molecules of 4a approach each other.

<Figure 12>

In a view of the MP2 results we conclude that the unimolecular mechanism is responsible for the formation of 5a. Note that the NH•••N links are typically responsible for the chemical exchange of the labile protons. Actually, compounds 4a and 5a were obtained by following the reaction path from TS4a along the intrinsic reaction coordinate. It should be noted that the reported energy barrier of ca. 40 kcal mol -1 is comparable to that found for the 1a → TS2a → 3a reaction. However, low energy barrier of the CO 2 detachment, and significantly higher stability of 5a as compared to 3a, proves that synthesis of 1a [10] could not be successful: this compound spontaneously decarboxylates when it is formed.

Structures of the TS3b and TS4b transition states in decarboxylation of 1b as well as imaginary eigenmodes (corresponding to frequencies equal to i332cm -1 and i2050cm -1 , respectively) are depicted in Fig. 13. Weaker intramolecular hydrogen bond in 1b as compared to that in 1a (cf. Table 1), is to some extent responsible for the lower energy barrier (ΔG ≈ 13 kcal mol -1 , cf. Fig. 9) of the dissociation process 1b  4b + CO 2 , and for the formation of the N1H bond in 4b. This stage of the reaction is expected to be spontaneous as well: energy lowering (ΔE = 4.8 kcal mol -1 , and ΔG = 17.24 kcal mol -1 ) was observed. The most significant difference between this (b) and previous (a) decarboxylation processes is an increase of the energy barrier associated with the migration of H9 towards C7 in the 4b  5b tautomerization by more than 20 kcal mol -1 . This value refers to both ΔE and ΔG. Note, that we assume the unimolecular mechanism of tautomerization, in spite of finding a relatively stable dimer of 4b at the B3LYP level, analogous to the dimer of 4a, for which the CP corrected binding energy is equal to ca. 6 kcal mol -1 (the imaginary eigenmode for the transition state found for (4b) 2  2 5b reaction corresponds to i1149cm -1 , the E/G energy barriers are ca. 14.7/12.6 kcal mol -1 , and the E/G energetic effects are ca. 7.5/13.6 kcal mol -1 per molecule, respectively). However, this time the lengthy MP2 calculations were not carried out. The reason for this can be easily understood when structures involved in the described processes are analyzed. The change in the hybridization of C7 to sp 3 when going from 4b to 5b via TS4b is associated with breaking of the planar symmetry of 4b. Similar behavior is also observed in reaction 4a  5a. One should keep in mind, however, that the C3C3' bond makes the molecules 4b and 5b rigid (cf. Section 3.1), and this results in a steeper energy increase under torsional distortions. As a consequence, significantly risen energy barrier of the 4b  5b step (ca. 60 kcal mol -1 ) may preclude tautomerization to proceed.

Thus, 4b may appear to be the final product of decarboxylation of 1b. In addition, it can be seen from Fig. 9 and Table 3 that although 5b is thermodynamically more stable than 5a, kinetic control of the 4a  5a process is more distinct than that of 4b  5b.

<Figure 13>

Finally we would like to point out that the bimolecular mechanism was also considered in the case of the enolization process of 1a. The stable dimer, that could initiate the H proton transfer from C7 atom of one molecule of 1a to O10 atom of the other molecule of 1a, was found at the B3LYP level. It is shown in Fig. 14. This time no transition state associated with such a migration was found (the MP2 calculations were not carried out).

<Figure 14>

Conclusions

Tautomerization and decarboxylation of 2,2-di(pyridin-2-yl)acetic acid were considered to be responsible for its instability. DFT/B3LYP/6-311G** calculations show that the products of decarboxylation of this acid and its rigid cyclic analogue, 1,8-diazafluorene-9-carboxylic acid, are more energetically favorable than the corresponding ene-1,1-diols (products of migration of the methine proton towards the carbonyl oxygen atom) and enaminones (products of migration of the same proton towards the pyridine nitrogen atom). The later tautomers were found to be thermodynamically more stable than enediols. Aromatization of the five membered ring has no effect on stabilization of ene-1,1-diol of 1,8-diazafluorene-9-carboxylic acid. Detachment of CO 2 from 2,2-di(pyridin-2-yl)acetic and 1,8-diazafluorene-9-carboxylic acids to give (E)-2-(pyridin-2(1H)-ylidenemethyl)pyridine and its cyclic analogue, respectively, may be followed by tautomerization of these by-products into (dipyridin-2-yl)methane and 1,8-diazafluorene. The calculated energy barriers for the CO 2 detachment correspond to quanta of IR radiation, so they can be easily overstepped in the course of reactions expected to afford 2,2-di(pyridin-2-yl)acetic and 1,8-diazafluorene-9-carboxylic acids. Thus, spontaneous decarboxylation of the former is in line with the unsuccessful synthetic efforts.

Despite their instability, the labile compounds may appear as the intermediates in numerous (bio)chemical processes. Conclusions of the present paper may be valid also for other substituted acetic acids (and their homologues) that carry electronegative heteroatoms in the gamma positions, such as diacylacetic acids and their derivatives (to be also known as unstable compounds).

Figure captions

Fig. 1 Relation between carboxylic acids and ene-1,1-diols

Fig. 2

Cyclopentadiene-1-carboxylic acid and its enol

Fig. 3 Fig. 4 Fig. 5 Fig. 6 Fig. 7 Fig. 9 Fig. 10 Fig. 11 Fig. 12

 345679101112 Fig. 3Relation between 1,8-diazafluorene-9-carboxylic acid and its enol

Fig. 13 Fig. 14

 1314 Fig. 13 Transition states TS3b and TS4b. The arrows show the atomic displacement corresponding to the imaginary frequencies

  Figure 1 Click here to download line figure: Fig. 1 eps.eps

  Figure 2 Click here to download line figure: Fig. 2 eps.eps

  Figure 3 Click here to download line figure: Fig. 3 eps.eps

  Figure 4 Click here to download line figure: Fig. 4 eps.eps

  Figure 5 Click here to download line figure: Fig. 5 eps.eps

Figure 6

 6 Figure 6 Click here to download line figure: Fig. 6 eps.eps

9

 9 Figure 7 Click here to download high resolution image

  Figure 10 Click here to download high resolution image

Table 2

 2 Selected optimized (B3LYP/6-311G**) bond lengths [Å] as well as valence and torsion angles [deg] in the transition states

	Bond/angle	TS1a	TS1b	TS2a	TS2b	TS3a	TS3b	TS4a	TS4b
	O9H9	1.009	0.984	0.994	0.982	1.829	1.549	-	-
	N1H9	1.754	2.078	1.721	1.883	1.030	1.093	1.263	1.356
	O10H7	1.246	1.345	2.881	2.861	-	-	-	-
	N1'H7	2.829	2.596	1.271	1.358	-	-	-	-
	C8O10	1.265	1.260	1.216	1.208	1.192	1.204	-	-
	C8O9	1.290	1.307	1.338	1.350	1.221	1.262	-	-
	C7H7	1.526	1.455	1.583	2.115	1.085	1.093	1.091	1.086
	C7H9	2.459	2.464	1.362	2.355	-	-	1.608	1.609
	C7C2	1.485	1.456	1.463	1.417	1.432	1.464	1.467	1.435
	C7C2'	1.493	1.488	1.490	1.434	1.492	1.499	1.455	1.460
	C7C8	1.489	1.452	1.495	1.475	1.951	1.692	-	-
	N1C2	1.354	1.336	1.354	1.357	1.371	1.335	1.359	1.343
	N1C6	1.337	1.342	1.335	1.322	1.356	1.351	1.343	1.338
	N1'C2'	1.345	1.328	1.351	1.427	1.345	1.327	1.353	1.341
	N1'C6'	1.335	1.338	1.341	1.344	1.333	1.339	1.335	1.331
	C2C3	1.403	1.418	1.413	1.437	1.416	1.411	1.399	1.411
	C3C4	1.387	1.391	1.384	1.403	1.374	1.388	1.387	1.399
	C4C5	1.395	1.396	1.397	1.381	1.414	1.406	1.405	1.401
	C5C6	1.388	1.396	1.387	1.415	1.367	1.383	1.384	1.403
	C2'C3'	1.407	1.427	1.390	1.483	1.405	1.422	1.413	1.444
	C3'C4'	1.386	1.391	1.392	1.382	1.388	1.391	1.383	1.394
	C4'C5'	1.395	1.392	1.399	1.400	1.392	1.392	1.398	1.388
	C5'C6'	1.390	1.398	1.386	1.394	1.392	1.398	1.390	1.404
	C3C3'	-	1.455	-	1.420	-	1.460	-	1.452
	C2C7C2'	117.8	104.6	119.9	105.4	122.7	101.4	121.3	102.0
	C2N1C6	119.7	116.0	119.9	116.7	124.1	120.8	123.3	118.4
	C2'N1'C6'	118.9	116.2	123.1	115.6	118.8	116.2	118.5	116.5
	C2C7C8	113.7	121.0	123.1	127.0	104.2	107.7	-	-
	C2'C7C8	112.3	127.7	111.0	127.2	101.7	122.7	-	-
	N1C2C7C8	-30.2	10.3	23.2	-5.2	54.7	-41.3	-	-
	N1'C2'C7C8	-1.7	-22.8	-95.8	13.2	104.5	62.8	-	-
	N1C2C7C2'	104.3	163.7	173.6	-178.8	-59.6	-171.3	138.3	-157.3
	N1'C2'C7C2	-136.7	-173.7	110.5	-173.2	-139.9	-177.4	-12.6	-176.9

Table 3

 3 Total energies (in a.u.) and ZPE corrections (in kcal mol -1 ) of the substrates,

	products and transition states
		E	ZPE
	2a	-723.447573	125.63
	3a	-723.449268	126.14
	↑
	TS1a	-723.350792	122.21
	TS2a	-723.370685	122.49
	↓
	1a	-723.447339	126.17
	↓
	TS3a	-723.413217	124.57
	↓
	4a + CO 2	-723.443569	123.81
	↓
	TS4a + CO 2 -723.369568	119.72
	↓
	5a + CO 2	-723.465199	123.80
	2b	-722.252866	113.08
	3b	-722.259203	113.25
	↑
	TS1b	-722.146161	108.67
	TS2b	-722.135826	108.65
	↓
	1b	-722.240663	112.39
	↓
	TS3b	-722.213498	110.27
	↓
	4b + CO 2	-722.245186	110.43
	↓
	TS4b + CO 2 -722.136557	105.69
	↓
	5b + CO 2	-722.262235	110.10


P Bailly

F Delvare

J L Hanus

Experimental study of the behaviour of building structural elements under soft impact

Dynamic three-point bend tests are performed on small beams specimen made of a quasi-brittle material. The most classical hypotheses needed to analyse the tests with SHPB, which are based on a state of quasi-static balance of the body, cannot be used. In order to get the behaviour it is necessary to perform a transient analysis of the specimen response. In most cases, the useful duration of the test -the elapsed time between the beginning of the load and the total failure -is lower than the time needed by the waves to reach the supports. The test is thus called a "one point bending test". In the case of bending failure, an analysis of the tests is possible by using only the measures given by the input bar. The technique is based on the knowledge of the analytical solution of the transient elastic response of the specimen. The failure mechanisms involved are specific features of the dynamic response. Different failure modes can be observed according to the loading rate: bending or shear fracture, single or multiple fractures. The determination of the bending moments and the rotations of the beam section, where the failure occurs, lead to an estimation of the strength and the failure energy.

INTRODUCTION

Structural elements such as walls, beams and columns can be damaged by shocks. The study of the behaviour of structural components under soft shocks is possible by using the Hopkinson bars (SHPB). Dynamic three-point bend tests are carried out on specimen made of brittle materials (like concrete, bricks ...). The experimental results show that, in many cases, the failure occurs for small strains and before static equilibrium is achieved in the specimen.

EXPERIMENTAL DEVICE

In order to look at the behaviour of quasi-brittle materials under dynamic loads, there are many different types of testing devices [START_REF] Mazars | Dynamic behavior of concrete and seismic engineering[END_REF]. Using a Split Hopkinson Pressure Bars system or Kolsky bars is also a possible choice [START_REF] Zielinski | Stress behaviour of concrete and mortar at high rates of tensile loading[END_REF] [3] [START_REF] Gary | Behaviour of a quasi-brittle material at high strain rate. Experiment and modelling[END_REF]. The device consists mainly of three bars: an input (or incident) bar impacted by a striker (short bar) coming from an air gun and two parallels output (or transmission) bars. The beam specimen is placed between the incident and the transmission bars. « Three point » bending tests, using SHPB, are usually performed with notched specimen for dynamic fracture toughness tests [START_REF] Yokoyama | Determination of Dynamic Fracture-Initiation Toughness Using a Novel Impact Bend Test Procedure[END_REF][6] [START_REF] Richomme | Dynamic testing of concrete with SHPB[END_REF]. In this article, the specimens are unotched: the figure 1 shows pictures of specimen set onto the device. The strain waves in the bars are measured by resistive gauges glued on the input and the output bars. The velocity of the striker is measured just before the impact. The impact velocity and the impact force applied on the specimen are evaluated using classical formulas (1):

( ) ( ) ( ) ( ) t ε t ε C = t V r i B e - - ( ) ( ) ( ) ( ) t ε + t ε E S = t F r i B B e - (1)
The tests, presented in this paper, are carried out on micro concrete and brick specimens. The geometric and material parameters of both specimens are specified in table 1. It has been observed that these tests are highly reproducible. The load reactions at the supports are derived from the measurements on output bars.

MODELLING OF THE TEST

Analyzing the bending test is possible if and only if stresses and strains characterizing the material behaviour can be derived from the data (force and displacement at the impact point). A relevant model is therefore required to represent the specimen motion and to estimate the internal stresses. This modelling must be time and space adjusted. The characteristic time of the studied phenomenon is greater than the propagation time of a transverse wave within the specimen. The Euler Bernoulli's beam model ( 2) is used to represent the specimen mechanical response. The authors use the most A key feature of the bending tests performed on quasi-brittle materials is that the time needed for a specimen to fail is lower than the time required for a bending wave to reach the supports: the load reactions at supports remain very small and appear only after the specimen failure. Therefore the problem is equivalent to a problem involving a beam with an infinite length. It is linked to the fact that there is no feedback information from the supports before fracture. Boundary conditions, at the impact point, are given by [START_REF] Brara | Experimental characterization on concrete in dynamic tension[END_REF]. Functions U e (t) and F e (t) are derived from the measured incident and reflected waves using the relations (1). Functions ψ(t) and M (t) are unknown:

0 = t w α 4 + x w 2 2 4 4 4 ∂ ∂ ∂ ∂ EI ρS = α 4 4 [ + , 0 [ x ∞ ∈ (2) ( ) ( ) ( ) ( ) ( ) ( ) t 0, M EI 1 = t 0, x w t ψ = t 0, x w t U = t 0, w 2 2 e ∂ ∂ ∂ ∂ - ( ) ( ) t F EI 2 1 - = t 0, x w e 3 3 ∂ ∂ (3)
The transient dynamic elastic analytical response of the specimen can be derived from the measured boundary conditions (4):

( ) ( ) ( ) ( ) ( ) ( ) ( ) τ τ Ω τ - - τ - + τ τ Ω τ - - = ∫ ∫ d , x t G t G d , x t G t , x w 2 t 0 2 1 1 t 0 1 (4) ( ) ( ) ( ) ( ) ( ) τ α
τ = τ τ - π τ = ∫ ∫ d EI 4 F t G d t V t G 3 e t 0 2 t 0 e 1 ( )         α π = Ω t 2 x cos t 1 t , x 2 2 1 ( )         α π = Ω t 2 x sin t 1 t , x 2 2 2
Employing the analytical Euler-Bernoulli beam model, one gets many results from experimental data For example, we can deduce the bending moment along the beam. The figure 2.a represents the bending moment at the instant of fracture for a brick. We can also deduce the evolution of the bending moment in the section where the fracture occurs (x=0) and the opening of the crack (the angle ψ). The figure 2.b shows such a result (for the same brick as in figure 2.a). During the whole elastic phase the hypothesis ψ(t) = 0 remains valid. The computed results are very closed to this condition even if there exits a delay due to the imperfect initial contact between the input bar and the specimen. On this curve, it is easy to detect the failure and to estimate energy dissipation during fracture.

THE FAILURE MECHANISMS

Criterion for the quasi brittle fracture

The impact test on a beam made of a quasi brittle material can lead to several mechanisms of failure. The first one is the bending failure like in a quasi static test. The second one is the shear fracture mode which leads sometimes to multi fragmentation. Shear fracture mode and multi fragmentation can usually only occur in dynamic tests. For a quasi brittle material, the simplest failure criterion is the Rankine criterion. The failure occurs when the highest principal stress reaches a limit σ T . For a beam, the criterion can be expressed as a function of the generalized stresses like bending moment M and shear force V.

The stress tensor in a beam is given by.

        0 σ σ σ = σ 12 12 11 ( ) ( )        - - 2 T T 2 2 2 3 12 T T 2 11 x 1 σ V V = x 4 h bh V 2 3 = σ σ M M x = I M x = σ (5) h 2x = x 2 T M M M = k T V V V = k T T σ 3 bh 2 = V T 2 T σ 6 bh = M So the criterion is { } T III II I σ < σ , σ , σ Sup ⇔ [ ] ( ) ( ) 2 < x 1 k 4 + k x + k x 2 2 2 V 2 M M ax m 0,1 x       - ∈ (6)
This criterion is drawn in figure 3. We can distinguish two different parts. When k M = 1 and 0 < k V < 0.92, there is a bending failure and when 0 < k M < 1 and 0.92 < k V < 1, there is a shear failure. In quasi static tree point bending test, the shear failure mode is not possible. To obtain this failure mode, the following geometric condition should be verified

2 h < L bh FL 2 3 > bh F 4 3 2 ⇒
such a condition does not agree with the hypothesis of a beam structure.

The loading path

For a dynamic test, the loading path depends on the loading rate. We consider that the loading is an imposed force (this assumption is closed to the real loading).

( ) t β = t 0, V (7)
The beam model leads to ( )

t t π 3α β = t 0, M 3/2 T 2 T M t σ π α bh β 2 = M M = k T T V σ h b 2 t 3β = V V = k (8)
The relation [START_REF] Bacon | Dynamic fracture toughness determined from load-point displacement[END_REF] gives the loading path.

β h π 2bσ α 4 k = k T 2 / 3 V M (9)
According to the value of β, we can observe bending failure or shear failure. The loading rate corresponding to the transition (k M = 1 and k V = 0.92) is called β tr . The figure 3

MULTI FRAGMENTATION

With a high loading rate, we can observe a first failure followed by a second (or more) one. The figure 5a shows four bricks after impact tested with increasing loading rate. We can observe bending fractures (crack perpendicular to the axis of the beam) and shear fractures (crack inclined at 45° on the axis). The number of fractures increases with the loading rate. This effect may be explained by the elastic transient motion of a part of the specimen, after the first failure [START_REF] Audoly | Fragmentation of rods by cascading cracks[END_REF]. As an example the figure 5b presents, for the second brick, the bending moment at several times after the first fracture which occurs at x = 0.

CONCLUSION

This paper shows a procedure to study the three point bending response of quasi-brittle materials using SHPB. The analysis can always be performed even if the measure on the input bar is the only available one. For some brittle tested materials, the effect of strain rate on the behaviour is similar to the one observed in direct tension tests. This dynamic bending test is relatively simple because it does not require building specific specimens: experiments are directly performed on small building components (bricks…).

Figure 1 .

 1 Figure 1. (a) The SHPB device (aluminium bars φ 40), (b) a micro-concrete specimen, (c) a brick specimen

  according to the situation [8 -15].

Figure 2 .

 2 Figure 2. a) Bending moment M (x,T R ) at T R = 32 µs (first fracture time) for a brick specimen b) Bending moment M (0, t) versus rotation ψ (t)

Figure 3 .

 3 Figure 3. Criterion for a beam made of a quasi brittle material with several loading paths corresponding to several strain rates

Figure 4 .

 4 Figure 4. Specimens after the dynamic test. a) Bending failure for a relative loading rate d = 0.38. b) Shear failure for a relative loading rate d = 7.7.

Figure 5 .

 5 Figure 5. a) Some bricks after the test (strain rate 1: 57 s -1 , 2: 279 s -1 , 3: 341 s -1 , 4: 426 s -1 ) b) Specimen n°2, bending moment along the beam plotted at several times.

Table 1 . specimen geometries and material parameters

 1

		Concrete	Brick
	Length L	14 cm	20 cm
	Thickness a	4 cm	1.7 cm
	Width b	4 cm	6.5 cm
	Elastic modulus E Density ρ	12 GPa 2000 kg/m 3	7 GPa 1400 kg/m 3