Hassen Trabelsi

Pierre-Alain Yvars
email: payvars@supmeca.fr

Jamel Louati

Mohamed Haddar

P.-A Yvars

Evaluation of the effectiveness of the interval computation method to simulate the dynamic behavior of subdefinite system: application on an active suspension system

Keywords: Interval based simulation, Ordinary differential equation, MacPherson suspension system, Dynamic simulation

A new design approach based on methods by intervals adapted to the integration of the simulation step at the earliest stage of preliminary design for dynamic systems is proposed in this study. The main idea consists on using the interval computation method to make a simulation by intervals in order to minimize the number of simulations which allow obtaining a set of solutions instead of a single one. These intervals represent the domains of possible values for the design parameters of the subdefinite system. So the parameterized model of the system is solved by interval. This avoids launching n simulations with n values for each design parameter. The proposed method is evaluated by several tests on a scalable numerical example. It has been applied to solve parameterized differential equations of a Macpherson suspension system and to study its dynamic behavior in its passive and active form. The dynamic model of the active suspension is nonlinear but linearisable. It is transformed into a parameterized state equation by intervals. The solution to this state equation is given in the form of a matrix exponential. Three digital implementations of exponential have been tested to obtain convergent results. Simulations results are presented and discussed.

Introduction

Conventional design methodologies based on a loop "designsimulate-back at the initial stage in case of failure" appear to be increasingly obsolete. During the past years, many studies have been conducted to optimize design, more precisely at the pre-sizing step [START_REF]Modelica, modeling of complex physical systems[END_REF][START_REF]virtual design[END_REF][START_REF]Cosmosworks[END_REF]. This step is often expensive and requires higher computation times [START_REF]virtual design[END_REF][START_REF]Cosmosworks[END_REF][START_REF] Wang | Inflatable wing design parameter optimization using orthogonal testing and support vector machines[END_REF] in order to obtain the optimal [START_REF] Wang | Inflatable wing design parameter optimization using orthogonal testing and support vector machines[END_REF][START_REF] Croccolo | Design and optimization of shafthub hybrid joints for lightweight structures: Analytical definition of normalizing parameters[END_REF] values of design parameters. In this work we propose a new design procedure with a rigorous evaluation that allows designing a subdefinite system [START_REF] Telerman | Integrating symbolic mathematical computation and artificial intelligence[END_REF][START_REF] Shvetsov | Ushakov, object-oriented constraint programming environment based on subdefinite models[END_REF] with a less number of simulations. The new design approach is based on the Interval computation method [START_REF] Moore | Interval analysis[END_REF][START_REF] Jaulin | Applied interval analysis[END_REF][START_REF] Granvilliers | Symbolic-interval cooperation in constraint programming[END_REF][START_REF] Eldon | Global optimization using interval analysis: the one-dimensional case[END_REF][START_REF] Revol | Arithmétique par intervalles[END_REF][START_REF] Meslem | Guaranteed approximation of reachable sets for uncertain dynamical continuous-time systems[END_REF][START_REF] Croisard | Preliminary space mission design under uncertainty[END_REF][START_REF] Cabannes | An uncertaintybased approach to drive product preliminary design[END_REF] and technique of constraint satisfaction problem (CSP) [START_REF] Montanari | Networks of constraints: fundamental properties and application to picture processing[END_REF][START_REF] Brailsford | Constraint satisfaction problems: algorithms and applications[END_REF][START_REF] Yvars | Using constraint satisfaction for designing mechanical systems[END_REF][START_REF] Yvars | Optimization of mechanical system: contribution of constraint satisfaction method[END_REF][START_REF] Fischer | Strategie de conduite du calcul pour l'aide la decision en conception integree; application aux appareils a pression[END_REF][START_REF] Bensana | A generic approach for conceptual design based on object oriented and constraint logic programming[END_REF][START_REF] Scaravetti | Structuring of embodiment design problem based on the product lifecycle[END_REF][START_REF] Chenouard | Search heuristics for constraint-aided embodiment design[END_REF][START_REF] Vargas | An environment for mechnical design aid applications development[END_REF]. This study focus just on the evaluation of the effectiveness of the interval computation method to solve differential equations (type ODE) by intervals which represents an important step of the new design approach proposed. This study is limited to the design of linear or linearisable parametric systems, which can be reformulated in the form of state equation by intervals. The intervals represent the areas of possible values of design parameters of the system. In the new design approach, simulation step is integrated into the preliminary pre-sizing phase and the stages of ordinary design is rearranged.The paper is organized as follows: in Sect. 2, the various stages of the basic design and the new design p rocess were quoted and compared. Section 3 is dedicated to the presentation of interval computation method as a guaranteed method that can encompass all the solutions that satisfy requirements imposed. The idea is to apply this method to compute the dynamic simulation of a subdefinite system by solving its differential equations by intervals. In Sect. 4, the problem of exponentiation of a matrix interval and method of convergence of results is presented. Section 5 details the dynamic model of the MacPherson suspension system. A simple scal-able model of MacPherson suspension system is presented to check the validity of the interval computation method. The Sects. 5.1 and 6 of this paper are devoted to the application of the interval computation technique on a MacPherson suspension system in passive and active case. The interval computation technique was used to solve differential equations under uncertainties which have been applied to the values of the design parameters represented as intervals. Finally, conclusions are drawn in Sect. 6.1.2.

Design approach

Currently, the system design has to take into account most aspects of several disciplines (aerodynamics, thermal, crash, vibration ...) [START_REF] Pahl | Engineering design: a systematic approach[END_REF][START_REF] Roshanian | Effect of variable selection on multidisciplinary design optimization: a flight vehicle example[END_REF]. It is based on the art of compromise, i.e. accommodation of a large number of constraints. This aspect involves several multidisciplinary skills in many engineering areas. These imply the need for dialogue between all the different engineering areas during the design process, to reflect properly their specifications in the constraints of optimization, and to ensure that the generated solutions fulfil all the imposed constraints. In such cases, the human reasoning alone is not enough, because it becomes difficult to find arrangements between all requirements imposed from different disciplines, and to decide how to act to improve the chosen objectives. The complexity of the problem leads naturally to the use of robust optimization techniques to determine the best solutions. Optimization studies generate the production of a large number of simulations. So it is necessary to completely automate the launch of these simulations [START_REF] Eltimsahy | Computer simulation for a solar heating system[END_REF].

Stages of a usual system design process are described in Fig. 1. It is based on a loop "design-simulate-back to the ini-Fig. 1 General process of product design tial step in case of failure" [START_REF] Teorey | A logical design methodology for relational databases using the extended entity relationship Model[END_REF]. Generally, design process is split in three main steps: the first step is to provide a preliminary design of the system components that satisfy all the requirements from a given specification. During this step, the various performance specifications that the product has to satisfy are described by a set of analytical models. Indeed, the geometry plays an important role in the preliminary design phase of the system. So in this operation the design variables of the product are defined and the dimensions are fixed. The next step is the "test". It is devoted for evaluation and simulation, allowing the designer to reduce the space of decision where solutions exist, in order to make the best choice of design variables. This step is related to the number of simulations made. Many methodologies and software tools in different engineering areas [START_REF] Paredes | Methodology to build an assistance tool dedicated to preliminary design: application to compression springs[END_REF][START_REF] Yan | A methodology for creative mechanism design[END_REF][START_REF] Croccolo | Design and optimization of shafthub hybrid joints for lightweight structures: Analytical definition of normalizing parameters[END_REF] have been developed and applied to simulate the behaviour of a product. For example in electronics hardware description languages we find (VHDL, Verilog) which are intended to represent and simulate the behaviour and the architecture of a digital electronic system i.e. Modelica [START_REF]Modelica, modeling of complex physical systems[END_REF] for multi-physics simulation or CATIA [START_REF]virtual design[END_REF] for geometry, Cosmos Works [START_REF]Cosmosworks[END_REF] to analyse dynamic behaviour of the system using finite element method for dynamic simulations. The last step of the conventional design process consists in checking the requirements. In this step the designer makes a comparison between the responses of the product obtained by simulation of the system with parameters chosen by estimation and/or by expertise and the desired response. So the designer is faced to two situations. In the first case if the resulting behaviour of the system design has satisfied the constraints imposed in the specifications document then, the design parameters used in the simulation will be taken as a solution. In the second case if the system response does not satisfy the constraints imposed, then the designer has to change these parameters by taking into account the previous simulation and the same design steps should be repeated until obtaining the optimal solution.

The previous design cannot be retained as a good approach to ensure the optimization of the system design. The designer has to go through several simulations to determine an optimal solution without being sure if the obtained solution is the global optimum within solutions space. It comes from the fact that the number of simulations that can be done is limited by the time and cost constraints which leads to oversize the system. So in order to have a robust design, we proposed through this work a new design approach in which simulation step is integrated into the preliminary design phase as shown in Fig. 2.

In Fig. 2 we present the different steps of the new design process. In this method, dynamics requirements are integrated in the preliminary design before simulation step. And among the dynamics requirements we are interested in this study to integrate dynamics simulation by intervals. The new design process is made up of four main stages: the first stage is to express design variables by intervals. Here the choice of the design variables values does not require expertise. We can put a blind choice and use long intervals, but we cannot deny that expertise may reduce the computing time. In the next step, the designer identifies the requirements that must be satisfied and express these requirements as constraints on the design variables. Those constraints provide declarative descriptions of important requirements related to engineering objectives. Then all types of constraints are implemented in CSP (constraint satisfaction problem) model (in this study it is a dynamic model) in the preliminary design phase.Then the designer propagates those constraints in the intervals of design variables to frame the areas of parameter defining the system design. Here the role of CSP approach comes. In fact the CSP [START_REF] Montanari | Networks of constraints: fundamental properties and application to picture processing[END_REF][START_REF] Brailsford | Constraint satisfaction problems: algorithms and applications[END_REF] eliminates all values of design variables that do not respect the imposed requirements. The sets of values that remain represent the entire generated solution.

The constraint programming [START_REF] Montanari | Networks of constraints: fundamental properties and application to picture processing[END_REF][START_REF] Brailsford | Constraint satisfaction problems: algorithms and applications[END_REF][START_REF] Yvars | Using constraint satisfaction for designing mechanical systems[END_REF][START_REF] Yvars | Optimization of mechanical system: contribution of constraint satisfaction method[END_REF] is a programming paradigm that emerged in the 1980s to solve combinatorial problems of large sizes such as problems with planning and scheduling. It is a technique intended to solve the mathematical problem that looks for states or objects satisfying a number of constraints. This technique is extensively used to treat problems manipulating intervals. Indeed it serves to minimize the areas of design variables. So the idea is to apply the CSP approach to solve design problems which involve design variables that must be expressed by intervals. The CSP calculation is essentially based on the interval compu- tation method, so the goal of this study is to test the capability of this method to solve differential equations type (ODE).

We notice that the CSP has been extensively used in different areas of engineering including systems control [START_REF] Jaulin | Applied interval analysis[END_REF][START_REF] Meslem | Guaranteed approximation of reachable sets for uncertain dynamical continuous-time systems[END_REF] and preliminary design [START_REF] Fischer | Strategie de conduite du calcul pour l'aide la decision en conception integree; application aux appareils a pression[END_REF][START_REF] Bensana | A generic approach for conceptual design based on object oriented and constraint logic programming[END_REF][START_REF] Scaravetti | Structuring of embodiment design problem based on the product lifecycle[END_REF][START_REF] Chenouard | Search heuristics for constraint-aided embodiment design[END_REF][START_REF] Vargas | An environment for mechnical design aid applications development[END_REF][START_REF] Sellini | Veri-KoMoD : Verification of knowledge models in the mechanical design field, knowledge acquisition, modeling and management[END_REF][START_REF] Yannou | A comparative study of constraint programming techniques over intervals in preliminary design[END_REF][START_REF] Aldanondo | Mass customization and configuration: requirement analysis and constraint based modeling propositions[END_REF][START_REF] Larroud | Constraint based approach for the steady-state simulation of complex systems : application to ship control[END_REF][START_REF] Tcherctchian | Benefits and limits of a csp/lca approach for the ecodesign of complex systems. a case applied to a hybrid passenger ferry[END_REF][START_REF] Canbaz | A multiplayer simulation of a set-based design framework with CSP techniques[END_REF]. In many excellent scientic works [START_REF] Fischer | Strategie de conduite du calcul pour l'aide la decision en conception integree; application aux appareils a pression[END_REF][START_REF] Bensana | A generic approach for conceptual design based on object oriented and constraint logic programming[END_REF][START_REF] Scaravetti | Structuring of embodiment design problem based on the product lifecycle[END_REF][START_REF] Chenouard | Search heuristics for constraint-aided embodiment design[END_REF][START_REF] Vargas | An environment for mechnical design aid applications development[END_REF][START_REF] Sellini | Veri-KoMoD : Verification of knowledge models in the mechanical design field, knowledge acquisition, modeling and management[END_REF][START_REF] Yannou | A comparative study of constraint programming techniques over intervals in preliminary design[END_REF][START_REF] Aldanondo | Mass customization and configuration: requirement analysis and constraint based modeling propositions[END_REF][START_REF] Larroud | Constraint based approach for the steady-state simulation of complex systems : application to ship control[END_REF][START_REF] Tcherctchian | Benefits and limits of a csp/lca approach for the ecodesign of complex systems. a case applied to a hybrid passenger ferry[END_REF][START_REF] Canbaz | A multiplayer simulation of a set-based design framework with CSP techniques[END_REF] the CSP was applied as an effective modeling language allowing the integration of requirements in the preliminary design phase. Most of requirements treated in these studies are static and/or ecological. So we propose in this work to integrate dynamics requirements such simulation by intervals of dynamic behavior of a technical system.

For the next step, the designer makes a simulation by intervals obtained by the CSP approach in order to determine the system response. Then he chooses the optimal solution from the set of solutions obtained by the interval simulation. As shown in Fig. 3 the interval based simulation method provides not only a single evaluation of the system behavior but also a set of performance bounds within the actual response system resides.

Theoretically and using what we presented before, with the new design process we can obtain the system response by interval and frame the solution space thanks to the interval computation. However the current problem, is that the evaluation of this method to solve a differential equation system and determine the dynamic response has not yet been tested on a scalable example. So the main objective of this work is to evaluate the effectiveness of the interval based simulation method to simulate the dynamic behavior of a parametric system (i.e. Macpherson suspension system) by intervals.

Interval computation

The interval computation [START_REF] Moore | Interval analysis[END_REF][START_REF] Jaulin | Applied interval analysis[END_REF][START_REF] Granvilliers | Symbolic-interval cooperation in constraint programming[END_REF][START_REF] Eldon | Global optimization using interval analysis: the one-dimensional case[END_REF][START_REF] Revol | Arithmétique par intervalles[END_REF][START_REF] Meslem | Guaranteed approximation of reachable sets for uncertain dynamical continuous-time systems[END_REF][START_REF] Croisard | Preliminary space mission design under uncertainty[END_REF][START_REF] Cabannes | An uncertaintybased approach to drive product preliminary design[END_REF] has been used to solve problems of uncertainty, it allows taking into account the measurement error by replacing x value measured with imprecision ε by the interval [x -ε, x + ε]. For example, replace a value not exactly representable, as π , by an interval containing [START_REF] Eltimsahy | Computer simulation for a solar heating system[END_REF][START_REF] Moore | Interval analysis[END_REF][START_REF] Eltimsahy | Computer simulation for a solar heating system[END_REF][START_REF] Jaulin | Applied interval analysis[END_REF]. Since the introduction of the interval computation, many applications have been developed such as linear optimization, solving ordinary differential equations, treatment of uncertainties,... The objective of the interval computation is to provide results that contain certain value or all searched values; so we talk about guaranteed or validated results. Indeed, in the interval computation, we do not handle numbers, but intervals. The result of an operation between two intervals [START_REF] Revol | Arithmétique par intervalles[END_REF][START_REF] Meslem | Guaranteed approximation of reachable sets for uncertain dynamical continuous-time systems[END_REF]: x y, of a function f (z) is the smallest interval (or vector interval) in the sense of inclusion containing:

{(x y | x ∈ X, y ∈ Y ) ⇒ ( f (z)| z ∈ Z )} Example: [x; x] [y; y]
= min x y, x y, x y, x y ; max x y, x y, x y, x y

The smallest interval means the smallest interval containing all possible results of the operation applied to all elements x of X and all elements y of Y, all the possible results of 'f' applied to all elements z of Z. The interval computation continues to grow, but with different objectives, for example, it represent an accurate tool to determine the global optimal or all zeros of a continuous function [START_REF] Granvilliers | Symbolic-interval cooperation in constraint programming[END_REF]. In this paper the objective was the application of an interval computation on a real system to get its dynamic behavior, through the resolution of its differential equations by interval.

Solving linear differential equations by interval

Solutions and methods for solving differential equations have been discussed in many papers, especially the methods that use generic intervals [START_REF] Ward | Numerical computation of the matrix exponential with accuracy estimate[END_REF][START_REF] Moler | Nineteen dubious ways to compute the exponential of a matrix[END_REF] instead of fixed values in their calculation algorithm. Unfortunately the evaluation of those methods in real cases provided disappointing results, that often diverge. In this paper we considered the linear case hoping to get a better convergence. A linear differential equation system is defined by:

{y (t) = A.y(t) and y(0) = y 0 |y ∈ R n and A ∈ R n×n }
where A is a matrix whose coefficients a ii that can be the design parameters of a system are expressed by intervals. The solution of this differential system is: y(t) = ex p(t.A).y 0 . So the exponentiation of a real matrix can solve problems of linear ordinary differential equations (ODE) with an initial value (IVPs). The numerical calculation of the exponential of the matrix has been extensively studied [START_REF] Moler | Nineteen dubious ways to compute the exponential of a matrix[END_REF][START_REF] Goldsztejn | On the exponentiation of interval matrices[END_REF] to get more accurate and realistic results. Indeed, if an approximate calculation from the last relationship led to an approximate solution for a linear ordinary differential system, our interval analysis provides more rigorous framework. Then, the strict enclosure of the solution would be obtained by interval computation of the matrix enclosing the exponentiation of the matrix of interval A. The most obvious way to get an interval enclosure of ex p(A) is to evaluate the truncated Taylor series using interval computation and to bind the rest. But the problem is that the results obtained by the exponential of Taylor and the exponential Horner with intervals computation often diverge especially the truncated Taylor series. For this reason, three implementations of exponentials (Exponential of Taylor, Exponential of Horner, and Exponential of Horner with reduction) given in [START_REF] Goldsztejn | On the exponentiation of interval matrices[END_REF][START_REF] Goldsztejn | Including ordinary difierential equations based constraints in the standard CP frame work[END_REF] are tested in our study. Its use in conjunction with interval analysis that allows us to control rounding errors to properly criticize and evaluate the results obtained by each exponential in the example of MacPherson suspension system will discussed in Sect. 6.1.2. The evaluation of those three implementations of exponential was made by programming in Matlab with the intervals library IntLab [START_REF] Rump | INTLAB -INTerval LABoratory[END_REF].

Taylor Series

The Taylor series is applicable only when the matrix has fixed values (like any diagonalizable matrix in C and the exponential of a complex is known. This method is very easy to implement). Unfortunately when it comes apart, this definition is no longer valid (we can not define an interval of eigenvalue due to its discontinuous nature). The easiest way is to apply the Taylor series:

ex p(A) = ∞ k=0 A k k! (1)
Let's give an example to evaluate its accuracy [START_REF] Vargas | An environment for mechnical design aid applications development[END_REF] and show that the truncation of the Taylor series is not well suited to the interval evaluation, even if no truncation of the series is made. We compute the formal expression of the exponential of the matrix A(t) for all t ∈ [-9 ; -7] with the eigenvalue method for its upper and lower limits. We find a misleading result because the lower and upper bounds can't be improved. This means [ex p(A), ex p(A)] is the best enclosure of ex p(A).

A = t 0 1 2 exp(A) = [0.0001234 ; 0.0009119] 1 2.7183 7.3891
Now we make the calculation with the Taylor series truncated to 18 terms with a rest practically negligible that can be estimated less than:

ρ( A , n)×[-9, -7] = | A n+1 (n + 1)! × (1 -|A n+2 ) ×[-9 ; -7] (2)
We get disappointing results as we can see below: As shown in the example above, even with a fairly high order for the expansion of the matrix [A], the influence of the rest is insignificant. So the evaluation of the Taylor exponential shows that it enlarges the limits of a matrix in an intervals calculation. The reason for this bad behavior of Taylor series in assessment by interval is the loss of dependency between the different expression occurrences of variables that occur during each assessment interval. Two well-known techniques can help to reduce the pessimism in the evaluation of the Taylor series by interval. First, a matrix whose intervals are centered forms which can provide results especially for very small input intervals. Such centered form for the exponential of a matrix which wass proposed in [START_REF] Oppenheimer | Application of interval analysis techniques to linear systems. I. Fundamental results[END_REF][START_REF] Oppenheimer | Application of interval analysis techniques to linear Sys II. The interval matrix exponential function[END_REF][START_REF] Oppenheimer | Application of interval analysis techniques to linear system[END_REF]. However, this centered assessment dedicated to exponentiation of interval matrix is very complex and difficult to follow or implement. The second technique is to rewrite the formal expression in order to obtain a more proportionate approach to the interval evaluation (usually by reducing the number of occurrences of variables). For example, the evaluation of a polynomial in Horner form is known to improve the assessment interval [START_REF] Tcherctchian | Benefits and limits of a csp/lca approach for the ecodesign of complex systems. a case applied to a hybrid passenger ferry[END_REF]. It can be naturally applied to the Taylor series in the exponential of a matrix, and was already used [START_REF] Oppenheimer | Application of interval analysis techniques to linear Sys II. The interval matrix exponential function[END_REF] to calculate the exponential of a matrix center as requested in the centered form. Proof of the best adaptation of the Horner exponential is provided in the following paragraph. We can conclude then that Taylor exponential must be used when there are a certain homogeneity in terms of the matrix and that its standard is not very important (the number of terms in the series must be chosen from the relation n + 2 > A ). Its bad behavior is attributed to a high number of operations in intervals evaluation.

exp(A) = [-

Horner exponential

To reduce the number of operations, Goldsztejn [START_REF] Goldsztejn | On the exponentiation of interval matrices[END_REF][START_REF] Rump | INTLAB -INTerval LABoratory[END_REF] applied the polynomial factorization of Horner [START_REF] Goldsztejn | On the exponentiation of interval matrices[END_REF] in Taylor series. This method consists in the evaluation of the polynomial image, at a point of the Euclidean space and performing several divisions. Specifically, instead of having a sum with power to assess, a product of simple terms that will reduce the number of occurrences of the terms of the matrix, is obtained. So the computation time is reduced and the cost of the computing operations is minimized. The exponential scheme of Horner is the following:

ex p(A) = I + A. I + A 2 I + A 3 ... I + A n (3)
However the rest is the same as that of the Taylor exponential (this was expected given that the writing comes from the expression of Taylor series). If we take the example of the previous section, we find for n=18 (with a rest below 10 -3 ): We have now a better framework as we move the component 2.1 from [486.589 ; 487.898] to [-67.172 ; 68.481]. This is an improvement in the calculation of the exponential. In addition, it is shown that the loss of dependency due to the significant improvement of the evaluation range, offering more accuracy and less expensive calculations. As explained previously, the center evaluation (centered matrix) provides more distinct enclosures but for the resolution of our problem is not yet clear. Therefore we turn to another digital trick.

exp(A) = [-

Horner exponentially with reduction

The objective here is to reduce the intervals to get a better precision in the calculations. The digital trick [START_REF] Goldsztejn | On the exponentiation of interval matrices[END_REF] consists in writing the exponential as follows:

ex p(A) = A 2 s 2 s (4)
When doing so, both the norm of the matrix and the intervals are reduced to the value we want. The number of terms to be chosen is determined by the following equation (n + 2).2s. A and can therefore be reduced depending on the value, without significant degradation of the results. But we have to be careful, when we choose a too big s because in the calculation of the power, the computer can interpret this value as infinite, which irretrievably distorts the result. If we take the example of the matrix A we obtain for n = 18 and s = 40 a very satisfying results as we can see below:

exp(A) = [-0.000241 ; 0.000912] 0 [0.656588 ; 0.821840] 7.389
We get for the same number of terms and a reduced number of operations, better results than the two previous exponentials. It therefore seems a priori that the Horner exponential with reduction is a suitable method to solve a linear differential equation. But to date, this approach, developed by Goldsztejn in [START_REF] Goldsztejn | On the exponentiation of interval matrices[END_REF], has been tested only on small matrix(si ze = 2, 2) as previously shown. One of our goals is to assess it on a real dynamic system. However, to confirm these results, we will test all the previous types of exponentials to simulate the behavior of a MacPherson suspension system.

MacPherson suspension system

Suspension systems cover several areas of applications. This involves different operating environments, ranging from domestic uses (office chairs ...) to applications with high reliability and requirement (vehicle). Depending on the application for which is dedicated to the suspension system, one or more objectives of pre-sizing are often specified and must be optimized.

The roles of a suspension system are to support the weight of the vehicle to isolate the vibrations from the road and to maintain the traction between the tire and the road.

In today's automotive industry, the optimization [49,50] of a vehicle in terms of comfort and security plays an active role. This involves performing a large number of simulations to find more efficient and innovative solutions, taking into account the time constraint. Yet, despite the diversity of existing tools and methods, the simulation step in most complex systems is relatively long. One of the main purposes of this study was to apply the interval computation method on a MacPherson suspension system Fig. 4 to determine several solutions, which offer the possibility to choose the ideal settings from existing solutions which satisfies all the imposed constraints.

Dynamic model of the arm of Macpherson suspension system

This study is based on the dynamic model of Sohn et al. [START_REF] Goldsztejn | On the exponentiation of interval matrices[END_REF] for MacPherson suspension system shown in Fig. 4. It is composed of a quarter vehicle body, an axle and a tire, a coil spring, a damper, an axle, a load disturbance and a control arm. The quarter vehicle body is assumed to have only a vertical motion. If the mass of the control arm is ignored and the joint between the control arm and the car body is control arm. The values of z s and θ will be measured from their static equilibrium points.The sprung and the unsprung masses are assumed to be particles. The mass and the stiffness of the control arm are neglected. The coil spring deflection, the tire deflection and the damping forces are in the linear regions of their operating ranges. Considering (y A , Z A ), (y B , Z B ) and (y C , Z C ) denote the coordinates of point A, B and C, respectively, when the suspension system is at an equilibrium point. Considering the sprung mass be translated upward by z s , and the unsprung mass be rotated by θ in the counter-clockwise direction. Then, the following equations hold.

y A = 0 ( 5 ) z A = z s ( 6
)
y B = l B (cos(θ -θ 0 ) -cos(-θ 0 )) ( 7
)
z B = z s + l B (sin(θ -θ 0 ) -sin(-θ 0 )) ( 8
)
y C = l C (cos(θ -θ 0 ) -cos(-θ 0 )) ( 9
)
z C = z s + l C (sin(θ -θ 0 ) -sin(-θ 0 )) ( 10
)
Where θ 0 is the initial angular displacement of the control arm at an equilibrium point. Considering α = α + θ 0 . Then, the following relations are obtained from the triangle OAB.

l = l 2 A + l 2 B + 2l A l B cosα 1 2 (11)
l = l 2 A + l 2 B + 2l A l B cos(α -θ) 1 2 (12)
Where l is the initial distance from A to B at an equilibrium state, and l' is the changed distance with the rotation of the control arm by θ. Therefore, the deflection of the spring, the relative velocity of the damper and the deflection of the tire are, respectively

( l) 2 = (l -l ) 2 = 2a 1 -b 1 cosα + cos(α -θ) -2 a 2 l -a l b l . cosα + cos(α -θ) + b 2 l cosα cos(α -θ) ( 13
) ( l) = (l -l ) = b l sin(α -θ)θ 2 (a l -b l cos(α -θ)) 1 2 (14) z C -z r = z s + l C (sin(θ -θ 0 ) -z r (15)
Where,

a l = l 2 A + l 2 B , b l = 2l A l B

Equations of motion

Lagrange's method is used to obtain the equations of motion. Considering T, V and D denote the kinetic energy, the potential energy and the damping energy of the system, respectively. Then, these are

T = 1 2 m s Z s 2 + 1 2 m u y C 2 + Z C 2 (16) V = 1 2 k s l 2 + 1 2 k t (z C + z r ) 2 (17) D = 1 2 c p l 2 (18)
Substituting the derivatives of ( 9), ( 10) and [START_REF] Shvetsov | Ushakov, object-oriented constraint programming environment based on subdefinite models[END_REF][START_REF] Moore | Interval analysis[END_REF][START_REF] Jaulin | Applied interval analysis[END_REF] into [START_REF] Granvilliers | Symbolic-interval cooperation in constraint programming[END_REF][START_REF] Eldon | Global optimization using interval analysis: the one-dimensional case[END_REF][START_REF] Revol | Arithmétique par intervalles[END_REF] yields

T = 1 2 (m s + m u ) Z s 2 + 1 2 m u l 2 C θ 2 + m u l C cosθθ z s (19) V = 1 2 k s [2a l -b l cosα + cos(α -θ) -2(a 2 l -a l b l . cosα + cos(α -θ) + b 2 l cos(α -θ)) 1 2 ]. 1 2 k t [z s + l C (sin(θ -θ 0 ) -sin(-θ 0 )) -z r ] 2 (20)
D = c p b 2 l sin 2 (α -θ)θ 8 (a l -b l cos(α -θ)) (21)
Finally, for the two generalized coordinates q 1 = z s and q 2 = θ , the equations of motion are obtained as follows:

(m s + m u )z s + m u l C cos(θ -θ 0 )θ -m u l C sin(θ -θ 0 )θ 2 + k t [z s -z r + l C (sin(θ -θ 0 ) -sin(-θ 0 ))] = -f d ( 22
)
m u l 2 C θ + m u l C cos(θ -θ 0 )z s + c p b 2 l sin 2 (α -θ 0 )θ 4 a l -b l cos(α -θ) -k t l C cos(θ -θ 0 )(z s + l C (sin(θ -θ 0 ) -sin(-θ 0 )) -z r ) - 1 2 k s cos(α -θ) ⎡ ⎢ ⎣bl + d l c l -d l cos(α -θ) 1 2 ⎤ ⎥ ⎦ = -l B f a (23)
Where, c l = a 2 la l b l cos(α + θ 0 ) and d l = a l b lb 2 l cos(α + θ 0 ) Now, introducing the state variables as [x 1 x 2 x 3 x 4 ]=[z s z s θ θ ]. Then, 22-23 can be written in the state equation as follows:

x 1 = x 2 x 2 = f 1 (x 1 , x 2 , x 3 , x 4 , f a , z r , f d ) x 3 = x 4 x 4 = f 2 (x 1 , x 2 , x 3 , x 4 , f a , z r , f d ) (24)
where,

f 1 = 1 D 1 [m u l 2 C sin(x 3 -θ 0 )x 2 4 + 1 2 k s sin(α -x 3 ). cos(x 3 -θ 0 )g(x 3 )+c p h(x 3 )θ . -k t l C sin 2 (x 3 -θ 0 ). z(.) + l B f a cos(x 3 -θ 0 ) -l C f d ] ( 25
)
f 2 = 1 D 2 [m 2 u l 2 C sin(x 3 -θ 0 )cos(x 3 -θ 0 )x 2 4 + (m s + m u ). c p h(x 3 )x 4 - 1 2 (m s +m u )k s sin(α -x 3 )g(x 3 )+m s k t . l C cos(x 3 -θ 0 )z(.) + (m s + m u )l B f a -m u l C cos(x 3 -θ 0 ) f d ] (26)
And

D 1 = m s l C + m u l C sin 2 (x 3 -θ 0 ) ( 27
)
D 2 = m s m u l 2 C + m 2 u l 2 C sin 2 (x 3 -θ 0 ) (28) g(x 3 ) = b l + d l c l -d l cos(α -x 3 ) 1 2
(29)

h(x 3 ) = b 2 l sin 2 (α -x 3 ) 4 (a l -b l cos(α -x 3 )) (30) z(.) = z(x 1 , x 2 , z r ) = x 1 + l c [sin(x 3 -θ 0 ) -sin(-θ 0 )] -z r (31)

Linearization

In order to simplify of the resolution the equations of the system, ( 20) is linearized at the equilibrium state where x e = (x 1e , x 2e , x 3e , x 1e ) =(0,0,0,0). Then, the resulting linear equation is

x = Ax(t) + B 1 f a (t) + B 2 z r (t) + B 3 f d (t), x (0) = x 0 ( 32
)
where,

A = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 1 0 0 ∂ f 1 ∂ x 1 ∂ f 1 ∂ x 2 ∂ f 1 ∂ x 3 ∂ f 1 ∂ x 4 0 0 0 1 ∂ f 2 ∂ x 1 ∂ f 2 ∂ x 2 ∂ f 2 ∂ x 3 ∂ f 2 ∂ x 4 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣ 0 1 0 0 a 21 0 a 23 a 24 0 0 0 1 a 41 0 a 43 a 44 ⎤ ⎥ ⎥ ⎥ ⎦ B 1 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 ∂ f 1 ∂ f a 0 ∂ f 2 ∂ f a ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 l B cos(-θ 0 ) m s l C +m u l C sin 2 (-θ 0 ) 0 (m s +m u )l B m s m u l 2 C +m 2 u l 2 C sin 2 (-θ 0 ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ B 2 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 ∂ f 1 ∂z r 0 ∂ f 2 ∂z r ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 k t l C sin 2 (-θ 0 ) m s l C +m u l C sin 2 (-θ 0 ) 0 m s k t l C cos(-θ 0 ) m s m u l 2 C +m 2 u l 2 C sin 2 (-θ 0 ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ B 3 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 ∂ f 1 ∂ f d 0 ∂ f 2 ∂ f d ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 l C m s l C +m u l C sin 2 (-θ 0 ) 0 m u l C cos(-θ 0 ) m s m u l 2 C +m 2 u l 2 C sin 2 (-θ 0 ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ with a 21 = -k t l C sin 2 (-θ 0 ) D 1 a 23 = 1 D 2 1 1 2 k s b l + d l (c l -d l cosα ) 1 2 cos(α + θ 0 ) - 1 2 k s sinα cos(-θ 0 ) . d 2 l sinα 2 (c l -d l cosα ) 3 2 -k t l 2 C sin 2 (-θ 0 )cos(-θ 0 ) .[m s l C + m u l C sin 2 (-θ 0 )] + m u k s l C sinα sin(-θ 0 )cos 2 (-θ 0 ) b l + d l (c l -d l cosα ) 1 2 a 24 = 1 D 1 . c p b 2 l sin 2 α 4 (a l -b l cosα ) a 41 = - m s k t l C cos(-θ 0 ) D 2 a 43 = 1 D 2 1 1 2 (m s + m u ) k s cosα b l + d l (c l -d l cosα ) 1 2 - 1 2 (m s + m u ) k s sinα . d 2 l sinα 2 (c l -d l cosα ) 3 2 -m s k t l 2 C cos(-θ 0 ) .[m s m u l 2 C + m 2 u l 2 C sin 2 (-θ 0 )] + 1 2 (m s + m u )m 2 u k s l 2 C sinα sin(-θ 0 ) b l + d l (c l -d l cosα ) 1 2 a 44 = - 1 D 2 . (m s + m u ) c p b 2 l sin 2 α 4 (a l -b l cosα )
6 Dynamic simulation of a passive Macpherson suspension system 6.1 Transient state

Fixed values

We handle the case where we impose a displacement of 10 cm on the wheel (the differential equations are with no second term), so the system is written as following: The system is shaped to apply the exponential matrix. Thus, without introducing uncertainty in the values of MacPherson suspension parameters (the system is solved with the method of exponential eigenvalues).

x = A.x(t) W ith : x = ⎡ ⎢ ⎢ ⎣ z s z s θ θ ⎤ ⎥ ⎥ ⎦ ( 33
)
The numerical values of the design variables of the MacPherson suspension system used for the simulation of its dynamic behavior are represented in Table 1.We obtain the reference curves shown in Fig. 6.

The dynamic behavior of a Macpherson suspension system when excited by a displacement of 10 cm on the wheel is represented in Fig. 6. The displacement z s converge to zero, which is normal as the role of a suspension is to mitigate the vibrations caused by road in a minimum time. We notice that the chassis takes 8s to stabilize and to attain its equilibrium position, which is not a long period considering the magnitude of excitation and the fact that the type of MacPherson suspension system is a passive suspension in this case. This is due to the fact that the design parameters are well chosen specifically on the chassis stiffness k s which quickly stiffens the structure. Also we notice the variation of the rotation angle θ is very small, it takes 8 s to attain its equilibrium position like the vertical displacement z s . The variation of the linear and angular velocity (z s , θ ) is very small compared to its equilibrium position, which implies that the damping coefficient c p is well chosen. Although the dynamic behavior of the MacPherson suspension responded well for the prospects expected but we are not sure if the design parameters chosen are the best. So in next paragraph we introduced uncertainty in these parameters.

Uncertainties

In this sub-section the evaluation of the dynamic behavior of the MacPherson suspension system was studied with uncertainty. The uncertainties have been applied to the values of design parameters which define the dynamic model of the MacPherson suspension system. The system of differential equations is solved with the interval Based Simulation method. The same calculation used in Sect. 6.1.1 is repeated but with uncertainty in reference values to determine the behavior of the suspension system by interval. Typically, for a first test, a 0.2 % uncertainty on the stiffness k s and 0.4 % on the stiffness k t , 0.3 % uncertainty on the mass m s and 0.1 % on the mass m u and 0.5 % on the damping coefficient c p were imposed. The areas of these design parameters are as follows: The calculation has been done in MATLAB with the intervals library IntLab [START_REF] Bochev | A self-validating numerical method for the matrix exponential[END_REF] for a number of terms equal to 14 (with Taylor, Horner or Horner with reduction) [START_REF] Croisard | Preliminary space mission design under uncertainty[END_REF] and a reduction factor equal to 45. The results are shown in Fig. 7.

In Fig. 7 the behavior of the Macpherson suspension system when excited by an instantaneous displacement is represented. With the Interval Based Simulation we determine an envelope which includes a set of performance for displacement z s and rotational movement θ of chassis 1. For the calculations of the dynamic response of the MacPherson suspension with uncertainty, two types of exponential have been applied. The calculations that have been made with Taylor (Fig. 7a,b) and Horner exponential (Fig. 7c,d).

We find that the results shown in the previous figures do not correspond to reality, since a significant difference (about 10 30 ) identified between the curves calculated with setting values and curves calculated with uncertainty. So as mentioned in the Sect. 4, the calculation by interval with the use of Taylor exponential or and Horner exponential diverges which confirms the divergence shape of curves of displacements z s and rotational movement θ . However, we notice that curves obtained by Taylor deviates in a increasing manner in a single direction while curves calculated by the Horner exponential deviates also in two directions in a symmetrical manner which is close to reality since the references curves In Fig. 8 the Horner exponential with reduction is tested to calculate the dynamic behaviour of MacPherson suspension. We notice that the set of solutions obtained are converging and have the same appearance like the curve calculated with fixed values which confirms the theories in Sect. 4 [START_REF] Goldsztejn | On the exponentiation of interval matrices[END_REF]. To ascertain the validity of the algorithm, the uncertainties were set to zero to check that the same curves were well received. However, the results obtained with the use of the mathematical formulation of Horner exponentially with reduction reflect that when the range of variables is divided into multiple sections, it will significantly reduce the error propagation. So we can conclude that the expression of Horner exponential with reduction is the most adaptable to the mathematical formulation to make exponential calculation by intervals. This asserts that the divergence of results is not due to the Interval based Simulation method but to the mathematical formulation of the problem.

Harmonic state

Setting values

Just like what has been achieved previously in Sect. 6.1.1, the same process and calculations are repeated in the next, except for the excitation of the dynamic model which has been done besides the displacement of 10 cm imposed on the wheel with a sinusoidal force. Repeating the previous matrices, the system to solve is:

x = Ax(t) + B 1 f a (t) + B 2 F(t) k t + B 3 f d (t) (34)
With:

F(t) = 0.02.k t .sin(10t)Even with zero initial conditions.

And:

f a = f d = 0
The solution of this system is the following form

x(t) = e A.t .x 0 + t 0 e A.λ .B 2 . F(t -λ). k t .dλ (35)
After an integration by parts we find the following solution: The Behavior of a Macpherson suspension system when excited by a sinusoidal force is represented in Fig. 9. The stars point curves ( -***-) represent the calculations which have been done with fixed values. So according to these results in Fig. 8 we notice that the curves z s and θ have the expected behaviours. They are sinusoids curves that follow the road as they have the same frequency and nearly the same amplitude. In addition, there is an angle shift between z s and the sinusoidal force of excitation F. It caused by the difference time between the instant when the peak occurred and the instant when users have felt it. Fortunately this is a part of the role of the suspension to delay and reduce the abrupt changes of the road. We can therefore use these new curves as reference curves.

x(t) = e A.t .x 0 + B 2 .K ω 1+ A ω 2 e A.t -cosω.t - A ω sinω.t (36)

Uncertainties

Uncertainties defined in Sect. 6.1.2 are now included but with the same percentages. According to results in Sect. 6.2.1. We apply just the Horner exponential with reduction since the results obtained by the others exponential expressions diverge. Curves in Fig. 9 show that the curves simulated con-sidering fixed values (reference curves) for displacements z s and θ is placed between the boundaries curves. So the validity of the calculation method is verified. We notice that the set of the system trajectories have similar look shape as references curves. The displacements z s and θ are sinusoids that follow the trajectory road.

Dynamic simulation of an active Macpherson suspension system

Comfort and road handling of a passenger car can be improved by replacing its passive suspension by a controlled suspension. A large number of studies [START_REF] Li | Linear and non-linear skyhook damping control laws for active railway suspensions[END_REF][START_REF] Priyandoko | Vehicle active suspension system using skyhook adaptive neuro active force control[END_REF][START_REF] Choi | Field test results of a semiactive ER suspension system associated with skyhook controller[END_REF] exist on the effectiveness of the skyhook control policy along with other optimal control techniques. Most of these studies indicate that skyhook control is the optimal control policy in terms of its ability to isolate the suspended mass from the base excitations. So we implement a skyhook controller in the Macpherson suspension system. One method of generating the skyhook damping force is to remove the passive suspension (i.e., both the damper and the spring) and replace it with an active force generator. This can be achieved by using a hydraulic actuator, however, the resulting system is rather complex and requires a significant amount of power [START_REF] Priyandoko | Vehicle active suspension system using skyhook adaptive neuro active force control[END_REF].

Another approach to achieving skyhook damping is to use semiactive dampers. Semiactive dampers allow the damping coefficient, and so the damping force, to be varied between high and low levels of damping. Early semiactive dampers were mechanically adjustable by opening or closing a bypass valve. The only power required for the damper is the relatively small power to actuate the valve. Since a semiactive damper is used in this study, we notice that the force applied by the skyhook damper to the suspended mass has the following form

F SK Y = -C SK Y .z s (37)
With C SK Y is the controller damping coefficient.

Transient state

Fixed values

Considering the previous matrices except the force f a is nonzero in this case as the suspension system is active. So the force f a is replaced by the relation of F SK Y and the system to solve becomes:

x = Ax(t) + B 1 (-C SK Y .z s ) + B 2 z r (t) + B 3 f d (t) (38)
After substitution of Eq. ( 37), we obtain the following system: After calculation which was done without uncertainty and with controller damping coefficient C SK Y =800N.s/m we obtain results in Fig. 10.

x = Ax(t) (39)
Figure 10a, b, c, d show the comparison of step response of MacPherson syspension system with skyhook and without skyhook controller for outputs z s ,z s , θ and θ . From Fig. 10 we notice that when we use controller overshot decrease and system reaches to steady state faster. It is noticeable from these figures that the magnitude for the body vehicle decreases even its velocity and the vehicle resumes its equilibrium state in 3 s instead of 8 s. This is an improvement of the system response which confirm the efficiency of the skyhook controller to the reduce of vibrations caused by the roads. The velocity variation (Fig. 10c,d) of active behaviour of MacPherson suspension is too small compared to the velocity variation of the passive behaviour which reflects that the acceleration of the body vehicle is reduced.

Uncertainties

The selection of the optimal controller damping coefficient is crucial, since it determines the optimal manner of the MacPherson suspension behaviour. So according to excitation kind, we need to determine if the controller damping coefficient is well chosen or not. So 2 % of uncertainties on the controller damping coefficient(C SK Y =[784 ; 816]) is introduced and the uncertainties of the values of design parameters is setted to zero.

From curve in Fig. 11 we notice that even if we take a higher percentage of uncertainty on the damping controller coefficient the results obtained converge. This is due to the fact that the uncertainty was implemented in one design variable C SK Y .

Harmonic state

With the uncertainties introduced on the damping controller coefficient , the designer can make a good idea on the effect of the increasing or the reducing of this coefficient on the behaviour of the MacPherson suspension system. So the simulation with interval can help to make a good decision about the choice of design parameters (Fig. 12).

With the interval based simulation method, several dynamic simulations result from the uncertainty imposed on the variables of the dynamic model. These performance bounds and the dynamic behavior of the system under nominal value can provide more complete information for a designer. So relying on these simulations, the designer can determine the solution that can fulfill the constraints imposed in the specifications document, which makes the sizing of the MacPherson suspension system easier.

Conclusion

From the previous study of dynamic behavior of the MacPherson suspension system, validity of interval based simulation to solve differential equations has been proved. With the use of Horner exponential with reduction, the exponential divergence problem when the Interval Based Simulation is applied and solved and the results of simulation are more realistic. However, this type of calculation is orientated to solve linear ordinary differential equations to make dynamic simulations. Indeed, the interval based simulation compared to a classical simulation allows multiple simulations simultaneously and recesses to reduce of the computational burden. The Interval Based Simulation method provides not only a single evaluation of the system behavior but also a set of performance bounds. Simulation results demonstrate that these performance bounds give a better description of the dynamic behavior of MacPherson suspension system with uncertain parameters.

Fig. 2

 2 Fig. 2 Interval based design process

Fig. 3

 3 Fig. 3 Classical simulation vs. interval based simulation

Fig. 5 Fig. 4

 54 Fig. 5 Dynamic quarter vehicle model of a Macpherson suspension

Fig. 6

 6 Fig. 6 MacPherson suspension system response when excited with snap-displacement. a Vertical displacement of mass m s . b Rotation of mass m s around X. c Linear velocity of mass m s . d Angular velocity of mass m s

  k t = [199200; 200800]m u = [24.975; 25.025] c p =[1990; 2010]k s =[19960; 20040]m s =[495.5; 501.5]

Fig. 7

 7 Fig. 7 The behavior of Macpherson suspension system when excited by an instantaneous displacement (comparison between the response calculated with the use of Taylor exponential and Horner exponential). a Vertical displacement of mass m s calculated by applying Taylor exponential. b Rotation of mass m s around X calculated by applying Taylor exponential. c Vertical displacement of mass m s calculated by applying Horner exponential. d Rotation of mass m s around X calculated by applying Horner exponential have a sinusoidal pseudo-periodic form, thus confirming the theory outlined in Sect. 4.In Fig.8the Horner exponential with reduction is tested to calculate the dynamic behaviour of MacPherson suspension. We notice that the set of solutions obtained are converging and have the same appearance like the curve calculated with fixed values which confirms the theories in Sect. 4[START_REF] Goldsztejn | On the exponentiation of interval matrices[END_REF]. To ascertain the validity of the algorithm, the uncertainties were set to zero to check that the same curves were well received. However, the results obtained with the use of the mathematical formulation of Horner exponentially with reduction reflect that when the range of variables is divided into multiple sections, it will significantly reduce the error propagation. So we can conclude that the expression of Horner exponential with reduction is the most adaptable to the mathematical formulation to make exponential calculation by intervals. This asserts that the divergence of results is not due to the Interval based Simulation method but to the mathematical formulation of the problem.

Fig. 8

 8 Fig. 8 The behavior of Macpherson suspension system when excited by an instantaneous displacement (response calculated with the use of exponential Horner with reduction). a Vertical displacement of mass m s calculated by applying Horner exponential with reduction. b Rotational movement of mass m s calculated by applying Horner exponential with reduction

Fig. 9

 9 Fig. 9 The behavior of Macpherson suspension system when excited by an instantaneous displacement (response calculated with the use of exponential Horner with reduction). a Vertical displacement of mass m s calculated by applying Horner exponential with reduction. b Rotational movement of mass m s calculated by applying Horner exponential with reduction

Fig. 10 a0 0 1 a+m 2 u l 2 C

 1012 Fig. 10 Comparison between active and passive MacPherson suspension system response when excited with snap-displacement. a Vertical displacement of mass m s . b Rotation of mass m s around X. c Linear velocity of mass m s . d Angular velocity of mass m s

Fig. 11

 11 Fig. 11 MacPherson suspension system response when excited with snap-displacement. a Vertical displacement of mass ms calculated by applying Horner exponential with reduction.

Fig. 12

 12 Fig. 12 MacPherson suspension system response when excited with snap-displacement. a Vertical displacement of mass ms calculated by applying Horner exponential with reduction.

Table 1

 1 Design parameters of the MacPherson suspension system


Lucrezia Veggi
email: lucrezia.veggi@ltf.mw.tum.de

Julian D Pauw

Bernd Wagner

Oskar J Haidn

S Y Mpo

O Tatin

A Study on the Design of LOx Turbopump Inducers

Keywords:

In this paper, the design methodology for turbopump components is presented and the design process of the inducer for a LOX turbopump is discussed in detail. e design methodology is based on classical methods and results of more recent studies. With the design process presented here, three inducers have been designed and numerically investigated. e aim of the study is to analyze the influence of the incidence angle on the performance by linking cavitation behavior and internal flow phenomena of the inducer. Results show that with increasing incidence angle, the cavitation formation increases causing flow deflection in the blade passage. ese flow phenomena provoke flow separation on the blade pressure side and increase the tendency of the flow to separate near the hub at the trailing edge.

INTRODUCTION

e need of high thrust levels for future space transportation systems results in high combustion chamber pressures in liquid rocket engines. e advantage of pump-fed over pressure-fed engines consists in the requirement of relatively low pump inlet pressures, and thus propellant tank pressures.

e major part of the pressure required in the chamber is supplied by the pumps. e use of turbopumps results in a decrease of system weight and an increase of performance as compared to pressurized-gas fed systems. Both dimensions and performance of the turbopump depend on the engine cycle and the requirements of the combustion chamber. Due to the high rotational speeds in the pump, the fluid pressure might drop below the vapor pressure resulting in cavitation.

is flow phenomenon can lead to unstable operating conditions and pump failure. Rocket propellant fed turbopumps have inducers in order to avoid cavitation, improve the suction performance and reduce the propellant tank pressure and weight.

e main task of the inducer is to provide a modest increase in pressure upstream of the main pump, typically a radial impeller, which in turns prevents cavitation and allows sufficiently reasonable operating conditions in the main pump.

First documents describing the design of inducers for liquid rocket engines came from the United States in the s-s [ , ] where recommended ranges of relevant geometric parameters were given. Based on studies that have followed in the years, Acosta [ ] summarized in the early s results concerning inducer design, performance, flow phenomena as well as system stability. In the last two decades, numerous experimental and numerical studies have been conducted on the impact of geometric variations on the performance. Acosta [ ], Bakir [ ] and Cervone [ ] made some efforts to understand the influence of the leading edge sweep on the performance in cavitating regimes, whereas Coutier-Delgosha [ ] published results on the impact of the leading edge shape on the cavitation behavior. e effect of tip clearance on inducer performance has been investigated e.g. by Hong [ ] and Kim [ ]. e common characteristic of these studies is the focus on the effect of geometric variations on cavitation behavior and performance. However, a detailed analysis of the flow phenomena inside the blade passage, which is influenced by cavitation and has an impact on the inflow of the main pump, is rarely covered in the literature. A detailed study of flow phenomena inside rocket pump inducers has been conducted by Lakshminarayana [ ], however operating the inducer with air in non-cavitating regime in the experimental setup.

Within the frame of the research project KonRAT at the Technical University of Munich (TUM), where a LOX turbopump for a LH / LOX expander cycle engine is investigated [ ], a design methodology for turbopump components has been developed. e design process for the pump and especially the design of inducers are discussed in detail. e design methodology is based on classical methods derived from the NASA documents [ , ], however taking into account more recent studies for the choice of some design parameters.

e intention of the authors is to initiate a process for a be er comprehension of the link between the cavitation behavior and the complex fluid mechanics inside the blade passage. e approach applied to achieve this goal is to re-duce the inducer geometry to its simplest form in order to limit the parameters responsible for flow phenomena, isolate individual effects and dissolve their dependencies. e study presented in this paper is the first step taken towards this direction. ree flat-plate inducers have been designed with a variation of the blade angle at the inlet and have been numerically investigated. During the analysis particular attention has been paid to the influence of the incidence angle on the fluid mechanics inside the blade passage along with the occurence of blade cavitation.

. INDUCER DESIGN METHODOLOGY

. Pump Specification e inducers investigated in the frame of this study are designed to meet the requirements of a LOX turbopump for a LOX/LH expander cycle engine with a thrust in the order of magnitude of the VINCI upper stage engine. e operational parameters of the turbopump at the design point were extrapolated from other works, e.g. [ , , ], dealing with cryogenic expander cycle engines in the thrust class of VINCI.

e operational conditions on the pump side at design point are listed in Table . Depending on the engine requirements listed above, the specific speed (Eq. ) can be calculated and the pump efficiency estimated by means of empirical correlations taken from Gülich [ ]. On the base of the N s -D s -diagram, the specific diameter is obtained giving a first guess of the outer pump diameter. Additionally, the sha diameter and the impeller hub diameter at inlet are derived from similar published pump designs and kept fixed. According to Huzel [ ], if the suction specific speed at inlet exceeds , or in SI-units, an inducer is needed in order to improve the suction performance. Similar values are given also by Lakshminarayana [ ] decreasing the limit to , or in SI-units. Evaluating the suction specific speed of the pump based on Jakobsen [ ], see Equation , indeed confirms the need of an inducer for this particular pump. Having a first estimation of the main geometrical dimensions and the configuration of the pump, a more detailed design follows. Table summarizes the pump characteristics necessary for the design process.

N s = N • √ Q H 3/4
( ) . Design Process e preliminary design of the inducer is part of an iterative design process involving all relevant components of the pump such as inducer, impeller, diffuser and volute. e preliminary design of the pump matches the requirements and characteristics presented in Table and and starts with the sizing of the impeller by assuming some parameters that are iteratively adjusted.

N ss = N • √ Q N PSH 3/4 ( )
In the first loop of the pump design, a swirl free flow upstream of the impeller is assumed considering a pump configuration without inducer. Once the main geometrical dimensions of the impeller are calculated, these are used as input for the design of the inducer and the diffuser. From the second loop in the design process, once the inducer dimensions are known, velocities and flow angles at inducer exit can be calculated. Simplifyingly assuming that in the axial clearance distance between inducer and impeller the flow does not redistribute, the velocities and the flow angles at inducer outlet are used to correct the previous assumption of a swirl free flow upstream of the impeller. Figure shows schematically the design process.

e sizing of the inducer is strictly dependent on the geometrical dimensions resulting from the preliminary design of the impeller. Depending on the head rise of the inducer, one differentiates between low head (ψ ≤ 0.15) and high head (ψ > 0.15) inducers. In the presented design methodology the head rise capability of the inducer is an assumption made at the very beginning of the inducer design. Su on provides a guideline value for the pressure increase in the inducer of typically % of the total pump pressure rise [ ]. According to Logan, simple inducers with noncambered blading, so called flat-plate inducers, where the blade exit angle is equal to the blade inlet angle, have head coefficients limited to values around . [ ]. When the geometry becomes more complex by introducing e.g. cambered blading or a tapered hub contour, the head rise capability of the inducer increases.

Since the goal of this study is to be er understand the influence of single design parameters on the internal flow

Outputs:

-main geometrical dimensions -velocities and flow angles at inlet and outlet -airfoil design

Inducer

Outputs:

main geometrical dimensions -pressure and velocities at slected positions -3D geometry

Impeller

Assumption 1 st loop:

Outputs:

-velocities and flow angles an inlet and outlet -main geometrical dimensions for the meridional contour and the airfoil design

Volute

Figure .

Design process of the pump components and therefore on the performance of the inducer, the first design of the inducer has to be as simple as possible with complexity being introduced step by step. Consequently, the geometries used for this study are flat-plate inducers with cylindrical hub and tip, which means that the hub and tip diameters are constant from inlet to outlet. e head rise capability of the inducer must be sufficient to deliver the demanded % of the total pump pressure rise.

ese assumptions for the pump requirements, shown in Table , result in a head coefficient of . . For lowhead inducers the inducer discharge dimensions must match those of the impeller eye [ ]. e hub-to-tip diameter ratio at the inducer outlet is estimated by means of the d t1 /d t2 -N s -diagram given in the NASA document SP-[ ]. e tip diameter is derived at the impeller inlet and at inducer outlet with the hub diameter kept fixed e sizing of the inducer blade profiles is based on the assumption that the flow can be described by potential flow models.

ese models assume a two-dimensional, irrotational, steady flow of an incompressible inviscid flow through a two-dimensional blade cascade [ ]. e tridimensional inducer blade is derived from the design of blade cascades at different heights between hub and shroud. Once the geometrical dimensions and characteristics of the inducer are known, the conditions of the entering and leaving flows necessary for the cascade geometry design are easily calculated by means of velocity triangles. Figure shows the blade cascade geometry in a certain radial position with inlet and outlet conditions as well as characteristic flow parameters. Characteristic geometric parameters of the cascade are blade chord length C, blade spacing along the cascade axis S, blade angle β, incidence angle i and deviation angle δ. In order to obtain these parameters a number of assumptions is necessary.

First, the number of blades Z has to be fixed. Pumping high-density fluids requires: • Large blade chord lengths to provide adequate bending strength to withstand the high hydraulic loads. Moreover, the larger the blade chord length the higher is the residence time for the collapse of cavitation bubbles and for the gradual addition of energy before entering the impeller.

• Large blade spacing in order to obtain a small ratio of blade thickness to blade spacing to accommodate the blade thickness inside the cavity and to minimize the blockage.

ese requirements are equivalent to a requirement for low blade numbers. In the literature, different recommendations can be found: the document of the NASA [ ] suggests a number between two and five, Huzel [ ] and Lakshminarayana [ ] recommend three to four for be er suction performance. For this study a two-bladed inducer has been chosen.

In the next step of the design process, the inducer blade angle has to be defined. For a flat-plate inducer, once the blade angle at inlet is derived, it is simply kept constant along the blade chord. Using a flat-plate inducer, to generate an increasing pressure rise results in increasing inlet incidence angles with the consequence of high losses and large-scale cavitation. Gülich [ ] suggests inlet incidence angles between °and °. Acosta [ ] and Lakshminarayana [ ] just give blade angle ranges between °and °(axial definition). In the NASA document [ ], the ratio of incidence angle to blade angle (circumferential definition) ν B should lie in the range of . and . . For thin blades the ratio should be . and for thick blades . . e recommended design value is the mean of . .

In the design process two possibilities have been implemented to fix the blade angle:

. e ratio ν B is assumed and kept fixed along the entire radius . Incidence angles at hub and tip are assumed. e incidence angle at the hub must be greater than at the tip and the gradient in between is linear.

A further assumption in the design process is related to the solidity. e solidity is the ratio of blade chord length to blade spacing. For flat-plate inducers, the solidity is supposed to stay constant over the radius. According to Gülich [ ] and

Jakobsen [ ], for be er suction performance current practice requires solidity values between and . With the solidity kept constant over the radius, further geometry parameters such as blade chord, blade axial length and angular wrap angle θ, can be derived. However, the parameter values calculated at this point of the process could change depending on upcoming decisions related to the inducer geometry. e solidity and therefore the blade chord and axial length can be affected by the leading edge and the trailing edge shapes. Very o en, inducers present a swept back contour at the leading edge, which positively affects the suction performance and the cavitation behavior. When introducing a sweep back at the leading edge, the loss in solidity should be compensated with a corresponding increase in axial length in such way that the solidity is maintained at its full value. Similarly, the trailing edge geometry has an impact on solidity, blade chord and blade axial length.

To complete the set of geometrical parameters necessary for the tridimensional design of the blade, the blade thickness has to be defined. Both hydrodynamic and structural design considerations determine the blade thickness resulting in a tapered form with the maximum at the juncture with the hub and the minimum at the tip [ ]. Gülich [ ] suggests a ratio of the blade thickness at the tip to the tip diameter between .

and . , keeping in mind that these values are related to industrial pumps where manufacturing differs from that applied in aerospace. However, evaluating available data from the literature [ , , , ], this ratio seems to be observed also in aerospace applications.

. INDUCER GEOMETRIES

In this study, three flat-plate inducers with cylindrical hub and tip have been investigated. e primary goal is to investigate the influence of the blade angle on the flow and on the suction performance. In order to be er recognize the influence of a variation of the blade angle on the internal flow, the geometries were designed to have as many as possible geometrical parameters in common.

e inducers do not have a sweep back at the leading edge. Nevertheless, the blade geometry does include an increase in axial length, since it is intended to introduce a sweep back at the leading edge in future studies and use the geometries presented here as a reference for the comparison of suction performance. All geometries have been designed to present the same angular wrap angle θ T E at the trailing edge at any blade height. A thickness ratio of .

has been applied resulting in a blade thickness at the tip of . mm. For all inducer geometries, the same value of blade chord length is applied at the tip.

e tip clearance of all three inducers has been set to a constant value of . mm which corresponds to . % of the blade height. e decision was taken considering the work of Jakobsen [ ] who recommends a maximal clearance area of % of the flow area. Here, the clearance area results in . % of the flow area.

As mentioned above, the only parameter that has been actively varied is the blade angle at inlet. Nevertheless, other parameters such as blade chord, blade axial length and angular wrap angle, are functions of the blade angle. A variation of its value inevitably leads to a variation of the related parameters. Moreover, the constraint of constant angular wrap angle over the radius at the trailing edge (see previous paragraph) results in different trailing edge shapes.

For the design of the first geometry (inducer A), the blade angle has been defined by assuming a ratio ν B of . , fixed throughout the entire blade height. e blade chord length at the tip of inducer A is the reference value for the design of the other geometries. For inducer B the ratio ν B has been decreased to . , which, according to Jakobsen [ ], is the minimum recommended value. In this way the incidence angles are decreased and the blade angles increased. For inducer A and B, the parameter variations result in a slightly backward swept trailing edge geometry as shown in the meridional view in Figure .  e incidence angles of inducer A and B are far from the recommended range of °… °. erefore, the geometry of inducer C has been designed with the second approach for the blade angle definition, namely by assuming an incidence angle of °at hub and of °at tip. . .

It is also clear, that in a flat-plate inducer an increase of the blade angle at inlet while keeping the blade chord length at the tip constant inevitably results in a decrease of blade axial length. All inducers present an elliptic geometry at the leading edge as well as at the trailing edge.

e main constructive parameters are summarized in Table .   . COMPUTATIONAL METHODOLOGY e three-dimensional incompressible flow calculation of the inducer was performed in ANSYS ® CFX . . For all computations a block structured mesh of around , elements has been used. e mesh was created with the mesh generator ANSYS ® TurboGrid. Figure shows the applied computational mesh.

For all numerical simulations, the tip clearance was not implemented. Although being aware of the impact of tip clearances on the pump performance, the neglection of the tip clearance was motivated by the intention to suppress the tip vortex due to the following reasons. e tip vortex interacts with the main flow and this interaction would have interfered with the detection of small changes in the flow due to blade angle variations. Furthermore, tip clearance leads to tip clearance and tip vortex cavitation in the peripheral zone, e.g. at the casing. e presented inducers have a radial leading edge design with the consequence that the blade cavitation is limited to the periphery of the inducer [ ]. e superposition of blade cavitation and tip vortex cavitation makes the detection of changes of the blade cavitation due to blade angle variation very difficult.

e geometrical domain has been modeled as a rotating fluid domain.

e casing is specifed as a counter-rotating wall so that the wall boundary is assumed to be stationary with respect to the stationary frame and in relative rotation e inducers have been simulated in steady state at the design point. Since alternate blade cavitation is not expected, one single passage has been calculated using rotational periodicity for the connection between the periodic boundaries.

is approach is widely used for CFD in turbomachinery and has shown good agreement with experimental results in prediciting the correct trends of performance and the size and location of blade cavitation [ , ].

e inlet was positioned mm upstream of the leading edge and the outlet mm downstream of the leading edge. At the inducer inlet, total pressure and total temperature were set. e total pressure was set to a constant value of . bar. At the inducer outlet the mass flow rate was used as boundary condition with a value of kg/s. All simulations were made with a fluid model of liquid oxygen at K. e standard SST model was used for modelling turbulence. For all the numerical simulations cavitation has been taken into account by using the cavitation model implemented in ANSYS ® CFX, which is a multiphase model based on the use of the Rayleigh-Plesset equation to estimate the rate of vapor production. In order to achieve convergence while the cavitation model is turned on, the following steps have been undertaken: First, an incompressible solution without cavitation with a total pressure at the inlet much higher than the desired one is computed. Here, the initial total pressure used was bar. en, from this non-cavitating initial solution, the cavitation model is turned on and the total pressure at the inlet is decreased by steps of . bar until the desired value is reached.

. RESULTS AND DISCUSSION e inducer geometries have been evaluated in terms of pressure rise capability, cavitation behavior and flow distribution in the blade passage.

Streamwise distribution of velocities and pressures

e absolute and axial velocities distributions as well as the total pressure distribution from inlet to outlet at midspan are summarized in Figure .  As expected, the pressure rise decreases with decreasing incidence angle, so that inducer C presents the lowest total pressure at outlet. However, comparing the resulting total pressures of the geometries with the design value, inducer C approximates the design value the most. It seems that for the design of a flat-plate inducer with the requirement of % of the total pump pressure rise, the recommended value for the ratio ν B of . results in incidence angles larger than necessary. Inducer B, which has been designed with a ratio ν B of . , also shows too high pressure rise. e choice of incidence angles in the suggested range of °to °is the most suitable option here. Looking at the development of the total pressure in axial direction in Figure , for all the geometries the most of the pressure rise occurs within the first / of the chord length, which is consistent with the measurements of Lakshminarayana [ ]. All inducers present a peak with a following decrease of total pressure. Absolute and axial velocity exhibit the same behavior, most pronounced for inducer A. is effect is related to the occurrence of cavitation on the suction side of the inducer blade.

Figure shows in the blade to blade view the distribution of the vapor volume fraction near the hub, at midspan and near the tip for the three inducers. A further analysis at the same spanwise positions has shown that in the region where the vapor volume fraction is near or equal to one the static pressure drops below the vapor pressure, confirming the presence of cavitation. Firstly, it can be observed that a variation of the blade angle has a much higher impact at the hub than at the tip. Comparing the distribution of vapor volume fraction of inducer A with that of inducer C at tip, the decrease of cavitation is visible but limited, whereas at the hub the variation of blade angle of inducer C results in a blade suction side almost free of cavitation.

Secondly, Figure shows that, independently from the geometry, cavitation covers the blade suction side up to the circumferential position where the adjacent blade begins. At midspan and below cavitation even extents into the blade passage. Consequently, the actual cross-sectional area at the entrance of the blade passage is reduced, increasing suddenly where the cavitation region ends. is causes a recirculation zone directly downstream of the cavitation region on the suction side and a flow deflection towards the axial direction which explains the strong increase of axial velocity and consequently of absolute velocity in Figure . Furthermore, the flow deflection from pressure side to suction side provokes downstream flow separation on the blade pressure side followed by rea achment of the flow.

e distribution of relative velocity at midspan in Figure (b) visualizes the flow phenomena described above. For all geometries, the axial velocity shows the biggest discrepancy in slope between numerical results and design value. In the hub region velocity shows a tendency of the flow to separate, while its maximum values develop near the tip. is observation is consistent with the experimental results of Lakshminarayana [ ]. e three geometries show similar behavior, however near the hub the distribution of axial velocity of inducer C is gradual, whereas for inducer A and B the axial velocity further strongly decreases. A redistribution of the axial velocity from the trailing edge towards the outlet can be observed for all the inducers especially in the near wall regions.

Spanwise distribution of velocities

According to Lakshminarayana [ ], the tendency to separate near the hub is caused primarily by radial flows in the wake and secondary flows inside the blade passage. Here, a clear correlation between the radial and the axial velocities could not be found. A closer look at the distribution of axial velocity on a surface directly downstream of the trailing edge, given in Figure , reveals localized zones of low axial velocities.

ese regions are located in the middle of the blade passage in the lower half of the blade height. e tangential velocity shows a slow increase towards the tip with a gradual distribution in the tip region for all three geometries. Generally, the tangential velocity is much higher than the design value and the discrepancy between design value and numerical results increases when the blade angle is reduced. Lakshminarayana [ ] justifies the high tangential velocity values with real fluid effects combined with radial flows that exist inside the blade passage.

Despite the different trailing edge shape, the inducers show very li le differences in the radial distribution of the velocity components and the velocity profiles are just shi ed towards lower or greater values.

. CONCLUSION

In the frame of the research project KonRAT at TUM a design methodology for turbopump components has been developed. In this paper the design process for the pump and especially for the inducer is described in detail. A numerical study on the influence of the incidence angle on the internal flow and on the performance was conducted on three flat-plate inducers with cylindrical hub and tip operating in liquid oxygen. In the numerical model, cavitation has been taken into account on the base of the Rayleigh-Plesset equation. e variation of incidence angle has been realized with two different approaches: e first two inducers were designed by keeping a fixed value of the incidence angle to blade angle ratio ν B over the entire radius, whereas for the third inducer incidence angles at hub and tip have been assumed. Results

show that an inducer designed with the recommended value of .

for ν B considerably exceeds the requirement of % total head rise capability. e bigger the incidence angle, the larger is the discrepancy between design value and numerical results in terms of total pressure at outlet and velocitiy distribution at exit. Downstream of the trailing edge, the flow shows a tendency to separate near the hub. is is due to flow deflections in the blade passage caused by cavitation on the suction side on the first / of blade chord length. With increasing incidence angle, the amount of cavitation rises, the flow deflections become stronger causing considerable decrease of axial velocity at outlet and flow separation on the blade pressure side. Big discrepancies between effective velocity profiles and design values at inducer exit, on which the impeller is designed, also increase the number of iteration loops in the design process.

e objective of this paper is to be er link the study of fluid mechanics of the internal flow with the cavitation behavior of an inducer. e first step in this direction is accomplished by reducing the inducer to its simplest form. Although the study is conducted on simple flat-plate inducer geometries, tendencies regarding the occurence of blade cavitation on the suction side as well as the flow behavior directly at the leading edge can be applied to more complex geometries such as cambered inducers which are more representative of the current technological standards. Surely, in order to confirm this supposition, the results obtained in this study and those obtained on similar geometries with cambered blading should be compared. It will be a further step to extend the design methodology by introducing the possibility to generate cambered inducer blades. e intention of the authors in starting from a very simple geometry and gradually introducing complexity is to be able to recognize the influence of single geometrical parameters in the complex flow behavior of state-of-the-art inducers so that targeted modifications of the geometry can be undertaken.

Further studies will cover the comparison between radial and swept back leading edges as well as the comparison between inducers with constant and variable hub contours. Furthermore, numerical simulations will be conducted by operating the inducers also in hot water in order to validate the numerical models with experimental data coming from the test facility at TUM [ ].

Figure .

 . Figure . Blade cascade geometry and flow parameters

Figure .

 . Figure . Comparison of the inducers geometries in the meridional view

Figure .

 . Figure . Computational grid of inducer A

  Figure shows a comparison of the three inducers in terms of tangential velocity (absolute), axial velocity and relative flow angle along the radius. e velocities and the flow angle are circumferentially mass flow averaged and evaluated directly downstream of the trailing edge and directly upstream of the outlet surface.

Figure . Figure . Figure . Figure . Figure .

 ..... Figure .Velocities and total pressure distribution from inlet to outlet at midspan

Table .

 . LOX pump operating conditions

	Pump requirements
	Rotational speed	,	rpm
	Nominal mass flow rate		kg/s
	Temperature at pump inlet		K
	Total pressure at pump inlet	.	bar
	Total pressure at pump outlet	∼	bar

Table .

 . Main constructive parameters of the inducers

	Flow coefficient [-]	.
	Head coefficient [-]	.
	Blade number [-]
	Hub diameter [mm]
	Tip diameter [mm]		.
	Hub-to-tip ratio [-]	.
	Blade chord length at tip [mm]
	Blade thickness at tip [mm]		.
		Inducer A Inducer B Inducer C
	Incidence-blade angle ratio ν B [-]	.	.	-
	Blade angle at tip [°]	.		.
	Blade angle at hub [°]	.	.	.
	Incidence at tip [°]	.	.
	Incidence at hub [°]	.
	Solidity at tip [-]	.	.	.
	Solidity at hub [-]	.	.	.
	Blade axial length at hub [mm]	.	.
	Angular wrap angle [°]

ACKNOWLEDGMENTS

is project is supported by Ludwig Bölkow Campus, funded by the Bavarian government. e good cooperation with the consortium partners is appreciated.


Anastasia Datsogianni
email: anastasia.datsogianni@psy.lmu.de

Stefan Ufer
email: ufer@math.lmu.de

Beate Sodian
email: sodian@psy.lmu.de

A pilot study on elementary pupils' conditional reasoning skills and alternatives generation skills in mathematics

Keywords: Conditional reasoning, alternative generation, mathematics, primary students

Background and rationale

One of the main aspects of logical reasoning is conditional reasoning, i.e. reasoning with ifthen statements. Conditional reasoning tasks usually present a rule of the form "if p then q" as a major premise, and a minor premise. Four minor premises differentiate four possible logical forms of inference: p is true (Modus Ponens, MP), p is false (Denial of the Antecedent, DA), q is true (Affirmation of the Consequent, AC), and q is false (Modus Tollens, MT). Definite conclusions can be drawn for MP ("q is true") and MT ("p is false"), while AC resp. DA do not allow definite conclusions. Even very young children can show conditional reasoning skills in familiar everyday contexts (Markovits & Thomson, 2008). However, a link between conditional reasoning and mathematics has been found only in adolescents and adults in the context of proof [START_REF] Stylianides | The mental models theory of deductive reasoning: Implications for proof instruction[END_REF].

According to Mental Model Theories, inferences are drawn by constructing mental models that encode information about specific situations in which the conditional is valid [START_REF] Johnson-Laird | Conditionals: a theory of meaning, pragmatics, and inference[END_REF]. The reasoners' ability to generate alternative models (beyond a model representing "p and q") for the given conditional is considered a crucial prerequisite to draw valid inferences. Studies on conditional reasoning in the everyday context [START_REF] De Chantal | The capacity to generate alternative ideas is more important than inhibition for logical reasoning in preschool-age children[END_REF] have shown that the alternative generation skills do predict early development of conditional reasoning. This raises the question, if the ability to generate multiple alternative models for a given mathematical premise, has an influence on students' conditional reasoning with these concepts.

Methods

Our study aimed at examining the relation between primary students' skills in conditional reasoning with mathematical concepts, and their performance on the corresponding, alternative generation tasks. Participants were 55 elementary students (4th graders n=13: M=9.5 years, 6th graders n=42: M=11.5 years) from a public school in Cyprus. Their conditional reasoning skills were assessed in two conditional reasoning tasks including mathematical concepts [START_REF] Datsogianni | Elementary students' conditional reasoning skills in mathematical and everyday contexts[END_REF]. Each task contained one item for each logical form (MP, MT, DA, AC). For example, one conditional reasoning task focused on the circumference of rectangles in a context of dwarf houses. It was explained that dwarf houses always consist of several aligned rows, all with the same number of same-sized quadratic rooms. The major rule in this task was "If a dwarf house has exactly 2 rows of 4 rooms each, then it has 12 windows". The alternative generation tasks asked students to generate many examples satisfying the conclusion of the major premise (i.e. "Draw as many dwarf houses as possible, that have 20 windows!").

Results and discussion

Students' solved 62.7% of the mathematics conditional reasoning items correctly, illustrating early conditional reasoning skills under specific conditions. Regarding the first alternative generation task 62% of students gave 4 to 6 correct alternative solutions (M= 3.42, Mdn= 4.00). The second alternative generation task seemed to be more difficult only 20% of students gave 4 to 5 correct solutions (M= 1.81, Mdn= 1.00). This indicates that the alternative generation tasks are feasible, but not too easy for our target sample. The number of generated alternatives correlated with students' reasoning performance (task 1: ρ=0.358, Sig.=0.001, task 2: ρ=0.343, Sig.=0.003), supporting our assumption that mathematical conditional reasoning is based on alternatives generation. Some first insights of our main study, which followed the above one, confirm that alternative generation in mathematics predicts correct logical reasoning and especially with AC form.

The results showed that the applied instrument is accessible to students while they replicate early conditional reasoning skills reported in studies from developmental psychology (e.g. [START_REF] Markovits | Different developmental patterns of simple deductive and probabilistic inferential reasoning[END_REF]. This study replicates also previous results (De [START_REF] De Chantal | The capacity to generate alternative ideas is more important than inhibition for logical reasoning in preschool-age children[END_REF] regarding the significant relation between alternative generation skills and students' conditional reasoning skills, extending these results to tasks that involve mathematical concepts. The small sample size forbids taking far-reaching conclusions. However, our future research will be based on these developed tasks to describe the role of mathematical knowledge in conditional reasoning with mathematical concepts in more details investigating also whether alternative generation training could be one approach to support students' conditional reasoning in mathematics.