Dr Alain Corfdir
email: corfdir@cermes.enpc.fr

Ph .D Jean Sulem

Comparison of extension and compression triaxial tests for dense sand and sandstone

Keywords: Compresssion test Extension test Friction angle Sand Compression test, extension test, friction angle, sand, sandstone

HAL is

After sending the first revised version, we noticed that reviewer # 3 had given an annoted version of our initial manuscript on the web site. So, this final revision takes his annotations into account. We apologize for this mistake.

The figures have been redrawn according to your remarks.

Introduction

Comparison of results from extension and compression tests on sands is an old dispute: is the friction angle at peak higher in extension or in compression ? This question is a particular point in the wider debate on the relationship between the friction angle and the Lode angle or the factor b=( 2 - 3 )/( 1 - 3 ) [START_REF] Lade Pv | Assessment of test data for selection of 3-D failure criterion for sand[END_REF].

Tests results collected by Wu and [START_REF] Wu W | On some issues in triaxial extension tests[END_REF] showed that for loose sand no general conclusion can be drawn and that for dense sand the friction angle is higher for extension tests than for compression tests performed at the same confining pressure. Special attention must be paid to the difficulties met when carrying out extension tests, the more specific one being the necking phenomenon which makes extension test unstable [START_REF] Lade Pv | Assessment of test data for selection of 3-D failure criterion for sand[END_REF][START_REF] Yamamuro | Strain localization in extension tests on granular materials[END_REF][START_REF] Wu W | On some issues in triaxial extension tests[END_REF]. Wu and [START_REF] Wu W | On some issues in triaxial extension tests[END_REF] listed the possible artefacts of extension tests such as the influence of the membrane, the effect of the weight of the sample, the effect of non-homogenous deformation. This last effect was corrected by measuring the local lateral deformation with a special device [START_REF] Kolymbas | A Device for Lateral Strain Measurement in Triaxial Tests with Unsaturated Specimens[END_REF] and evaluating the axial stress with a corrected area. Yamamuro andLade (1995, 1996) avoided necking of the sample by using membranes reinforced by metal plates as described in [START_REF] Lade Pv | Effects of shear band formation in triaxial extension tests[END_REF]. In this paper, the method of comparison of the measured friction angles is critically discussed in order to better understand the different values measured in the case of dense sands in dry or drained condition. High pressure tests where grain crushing can cause a change of behaviour [START_REF] Yamamuro | Drained Sand Behavior in Axisymmetric Tests at High Pressure[END_REF] are not considered in this study. The difference in measured friction angle appears to be due to the influence of stress level. Finally, the comparison is extended to cohesive granular media such as sandstones. For tests without brittle tension failure, the tangent friction angle is compared for compression and extension tests. We emphasize the fact that the comparison is made for the maximum mobilized friction angle so that only the failure state is discussed and not the complete constitutive behaviour of the material. This failure state may be the result of strain localisation. In particular, we are interested in discussing the appropriate stress quantity which controls the stress-dependent character of the friction angle at failure.

Comparison of experimental results for dense sand

Only a limited number of tests data allowing the comparison of extension and compression results are available in the literature. As stated above, extension tests are especially difficult to perform and we use here two data sets, one from (Kolymbas and[START_REF] Kolymbas | Recent results of Triaxial Tests with Granular Materials[END_REF][START_REF] Wu W | On some issues in triaxial extension tests[END_REF], and the other one from [START_REF] Yamamuro | Drained Sand Behavior in Axisymmetric Tests at High Pressure[END_REF].

For both data sets, the samples were built using a pluviation technique. Wu and [START_REF] Wu W | On some issues in triaxial extension tests[END_REF] tested a dry Karlsruhe medium sand with D 50 =0.395 mm and average relative density 100%. We also refer to the compression tests of [START_REF] Kolymbas | Recent results of Triaxial Tests with Granular Materials[END_REF] on the same material. [START_REF] Yamamuro | Drained Sand Behavior in Axisymmetric Tests at High Pressure[END_REF] tested a water saturated Cambria sand in drained conditions. This sand is uniformly graded with a grain size of 0.83 to 2 mm and an average relative density 89.5 %. Among the test data of Yamamuro and Lade, we selected those for which no important grain crushing was observed. This condition gives a limit of about 8 MPa for the shear stress ' a -' r where ' a is the (effective) axial stress and ' r is the confining stress. In the following we will consider total stresses for dry samples and effective stresses for saturated samples tested in drained conditions. For simplicity, the total stress (for dry samples) and the effective stress (for saturated samples) are both denoted by .

For cohesionless materials, the friction angle for each test has been simply evaluated in the Mohr plane as the secant angle at maximum shear stress sin /( )

a r a r        
. This definition does not imply a priori the validity of Mohr-Coulomb criterion, nor the independence from the intermediate principal value of the stress tensor.

The results of Wu and [START_REF] Kolymbas | Recent results of Triaxial Tests with Granular Materials[END_REF][START_REF] Wu W | On some issues in triaxial extension tests[END_REF] are presented in Figs. 1, and2 where the maximum friction angle in compression and in extension is plotted vs.

the confining stress  r and vs. the mean stress p = (2 r + a )/3 at peak. As shown by these plots, the difference between the friction angle in compression and in extension plotted against the confining pressure  r is reduced when the friction angle is plotted against the mean stress.

Figure 1

Figure 2 The difference is reduced even more by plotting the measured values of friction angle as a function of s = ( a + r )/2 (Fig. 3). The same remarks can be made for the test results by for a large range of pressures.

Figure 3

Figure 4 Figure 5 3

Comparison of experimental results on sandstone

The same approach is now applied to sandstone which can be viewed as a cohesive granular material. Two sets of tests results have been published by [START_REF] Bésuelle P | Déformation et rupture dans les roches tendres et les sols indurés : comportement homogène et localisation[END_REF] and [START_REF] Bésuelle | Experimental characterisation of the localisation phenomenon inside a Vosges sandstone in a triaxial cell[END_REF] for Woustwiller sandstone on one hand and by [START_REF] Sulem | Elasto-plastic modelling of Red Wildmoor sandstone[END_REF] and by [START_REF] Papamichos | Constitutive testing of Red Wildmoor sandstone[END_REF] on Red Wildmoor sandstone on the other hand. Their data are analysed in the following.

As it can be seen on Fig. 6 for Woustwiller sandstone, it is not possible to define a secant friction angle or a tangent friction angle for extension tests at relatively low confining pressure as for  r =10, 20 and 30 MPa, tensile failure occurs on a quasihorizontal plane [START_REF] Bésuelle | Experimental characterisation of the localisation phenomenon inside a Vosges sandstone in a triaxial cell[END_REF]. Thus, only the tests at higher confining pressure where shear failure occurs are considered hereafter. It is also observed on Fig. 6 that a unique envelope surface cannot be drawn for the Mohr circles of compression and extension tests with failure in the compressive regime.

Consequently, the discussion herein will refer to the tangent friction angle.

Figure 6

The tangent friction angle is evaluated by considering the tangent line to two adjacent Mohr circles. This quantity is plotted as function of the X-coordinate of the middle point of the centre of the two circles (Figs. 78). The plots of Figs. 7, and 8 show that the friction angles in compression and extension lay on a unique curve when plotted as a function of ( a + r )/2 (centre of the Mohr circle). The above results do not apply to loose sands. This has been checked for the tests of Wu and [START_REF] Kolymbas | Recent results of Triaxial Tests with Granular Materials[END_REF] and of [START_REF] Lade Pv | Relative density effects on drained sand behavior at high pressures[END_REF]. Fig. 9 shows the secant friction angle plotted versus ( a + r )/2 for tests by [START_REF] Lade Pv | Relative density effects on drained sand behavior at high pressures[END_REF] on Cambria sand with an initial relative density of 30%. Even at low confining pressure, extension tests and compression tests do not lie on a single curve. The gap between the two curves is even larger when the friction angle is plotted versus the radial stress  r or the mean stress p.

On Fig. 9, a clear change is observed for the slope of the curve at high confining pressure (>10MPa) for the compression tests. In this range of confining pressure the friction angle is larger in compression than in extension. The same phenomenon is also observed for dense sands at high confining pressure (Yamamuro and Lade, 2006). The change of response at high confining pressure can be attributed to particles breakage which is different in compression and in extension. (Yamamuro and Lade, 2006).

For these two cases, a simple stress quantity cannot be found for describing the stress dependency of the friction angle in compression and in extension. Haimson's observation corroborates our finding that sin

a r a r         is a unique function of  / 2 a r    .
/( ) ( b 3 1 3 2       
. In a quite similar way as for rocks, the octahedral shear stress seems to be better correlated with  

1 3 / 2    than with   1 2 3 / 3      (Fig 10).

Figure 10 4.3 Choice of ( a + r )/2 for describing stress-dependent friction

In triaxial condition, the study of the Mohr circles reduces to the study of the circle in any plane containing the axial direction. If the failure envelope can be drawn, for each value of secant angle there is a unique Mohr circle, the centre of which has its X coordinate equal to ( a + r )/2 for both extension and compression loading cases. The failure envelope is very simply described by the secant friction angle as a function of ( a + r )/2. It is the only natural choice which makes it possible to check immediately whether the failure criterion is the same in compression and extension.

Stress dependency of the friction angle is commonly expressed in terms of the mean stress. This invariant is chosen by [START_REF] Yamamuro | Drained Sand Behavior in Axisymmetric Tests at High Pressure[END_REF] for comparing compression and extension tests. For isotropic materials, a functional basis of independent invariants consists of 3 elements. The usual basis is kk 3

1 p   , ji ij s 2 s s 2 1 J  where ij ij ij p s    
is the deviatoric stress tensor and

kj jk ij s 3 s s s 3 1 J 
. The mean stress p is the only linear stress invariant for isotropic systems. For axisymmetric triaxial tests, it can be easily proved that 2 / ) ( r a    can be written in terms of stress invariants:

s 2 s 3 r a J J 4 3 p 2      (1)
Thus the choice of ( a + r )/2 is natural and is not in contradiction with invariant tensor theory.

Comparison of tests data and predicted values in the case of sands

               (3)
The axial stress can be expressed as:

)) s ( sin 1 ( s a     (compression); )) s ( sin 1 ( s a     (extension) (4)
Using equation ( 3), we get:

0.9355 1.24 (with , in MPa) 2 a r a r a r                (compression); (5a) 0.9355 1.24 (with , in MPa) 2 a r a r a r                (extension) (5b)
The predicted values of  a versus  r given by ( 5) are compared with tests results in Fig. 11.

Figure 11

Predicted values and experimental data appear to be in good accordance on a large range of confining pressure for both extension and compression tests, even for data corresponding to the two different sands (Karlsruhe sand and Cambria sand).

Failure mechanism and tests results

Different failure mechanisms induce different maximal shear values for triaxial tests (e.g. [START_REF] Suzuki | Double Strain Softening and Diagonally Crossing Shear Bands of Sand In Drained Triaxial Tests[END_REF] and for true triaxial tests [START_REF] Lade Pv | Assessment of test data for selection of 3-D failure criterion for sand[END_REF]). In the above related tests different failure modes (with no shear band, with an incomplete partially developed shear band, with a fully developed single shear band or with a system of several shear bands) have been obtained.

The tests data used for sand were not focussed on tracking shear band formation and evaluating its inclination. But Wu and [START_REF] Wu W | On some issues in triaxial extension tests[END_REF] mentioned shear banding in extension tests for dense sand. It is well known that dense sand triaxial tests involve shear banding [START_REF] Vardoulakis | Bifurcation Analysis in Geomechanics[END_REF]. Concerning tests on Woustwiller sandstone, for the tests that have been compared (i.e. compression and extension tests with similar value of ( a + r )/2), it has been found that the failure pattern is a unique failure plane in both cases. For Wildmore sandstone, all triaxial tests presented here have a shear band failure pattern [START_REF] Papamichos | Constitutive testing of Red Wildmoor sandstone[END_REF].

Thus for the tests considered in this paper, when the appropriate data are available, the failure mode is similar for compression and extension when the secant friction angles (for sand) or tangent friction angles (for sandstone) are identical in compression and extension.

In the pre-localisation regime, triaxial tests conditions impose fixed principal stress directions whereas strain localisation allows for principal stress rotations.

This phenomenon can be used for calibration out of plane shear moduli [START_REF] Vardoulakis I | Shear band inclination and shear modulus of sand in biaxial test[END_REF][START_REF] Sulem | Elasto-plastic modelling of Red Wildmoor sandstone[END_REF][START_REF] Desrues | Shear band analysis and shear moduli calibration[END_REF]. However an out of plane stress increment does not affect much the principal stresses. If we assume, as it is done here, that the failure criterion is isotropic then stress rotation induced by the strain localisation process does not affect the relation between the axial stress and the radial stress at failure. On the other hand, strain localisation can result in a stress-induced anisotropy of the material inside the shear band and could thus lead to an anisotropic failure criterion. A recent paper of [START_REF] Gajo | On certain critical material and testing characteristics affecting shear band development in sand[END_REF] shows that although stress induced anisotropy affects the conditions for onset of strain localisation, the global response of sample is not much affected.

Therefore, it is expected that stress induced anisotropy will not affect much the axial stress at failure.

The comparison of compression and extension triaxial tests on dense dry sands

shows that when special care has been taken to prevent or to correct the effects of necking, the secant friction angle at failure is a unique function of ( a + r )/2 for both loading cases. The comparison of compression and extension triaxial tests on sandstones shows also that the friction tangent angle is very similar in both cases (for the same value of ( a + r )/2), when failure occurs with shear banding. When strain localisation occurs, a transition between the imposed axisymmetric strain mode and a plane strain mode takes place and thus our results suggest that whether the failure plane is formed in extension or in compression does not affect the stress-dependency of the maximum mobilized friction angle. Therefore, constitutive models which incorporate stress-dependency of the friction angle should be compatible with this result. This result does not hold for loose sands or for sands at high pressures. In the case of loose sands, this can be explained by the fact that the failure mechanism is different from the one of dense sands as no peak is observed in the stress-strain curve and no failure plane is formed. For the case of sands tested at high pressures, one could relate the different response to the different mechanism of particle breakage in compression and in extension.

  Figure 7Figure8

Figure 9 4. 2

 92 Figure9

Fig. 1 Fig 3 Fig. 8

 138 Fig.1Secant friction angle as a function of radial stress; data from Wu and[START_REF] Wu W | On some issues in triaxial extension tests[END_REF],[START_REF] Kolymbas | Recent results of Triaxial Tests with Granular Materials[END_REF]

  Figure 1 Click here to download high resolution image

Ackowledgements

The authors wish to thank Prof. Dimitrios Kolymbas for fruitful discussions.


F Germain

P Y Rohan

G Rochcongar

P Rouch

P Thoreux

H Pillet

W Skalli

Role of Ligaments in the Knee Joint Kinematic Behavior: Development and Validation of a Finite Element Model

Introduction

The management of knee instability is a complex problem in orthopedic surgery. It usually involves the tear or the rupture of a cruciate or collateral ligament of the articulation and, in many cases, necessitates surgical operation and ligament reconstructions. The rupture of the ACL, in particular, is one of the most frequently occurring ligament injuries affecting about one person in 3000 every year and this trend has been constantly increasing with the rise of participation in sports in the general population [START_REF] Doroski | Techniques for biological characterization of tissueengineered tendon and ligament[END_REF]. Knee ligamentoplasty has become an issue of high clinical interest.

The current standard of care is based on ligament reconstruction by autografts from tendon tissues, allograft, or ligament substitutes. The assessment, however, of the long-term performance of the surgery is complicated by the complex interaction between the graft parameters (pretensioning, fixation method, etc.), knee anatomy, and the mechanical interaction between the graft and the passive anatomical structures. To this day, no diagnostic tool is available in clinics to quantitatively evaluate and predict the impact of surgery on knee kinematics. This highlights the need for the development of tools to investigate the contribution on knee joint mobility of the ligaments in the normal knee and to establish proper treatment strategies.

In that respect, many studies have been performed both in vitro [START_REF] Azmy | EOS ® orthopaedic imaging system to study patellofemoral kinematics: assessment of uncertainty[END_REF][START_REF] Li | Kinematics of the knee at high flexion angles: an in vitro investigation[END_REF][START_REF] Lorenz | The influence of asymmetric quadriceps loading on patellar tracking -an in vitro study[END_REF][START_REF] Nagamine | Effect of medial displacement of the tibial tubercle on patellar position after rotational malposition of the femoral component in total knee arthroplasty[END_REF] and in vivo [START_REF] Andersen | Do kinematic models reduce the effects of soft tissue artefacts in skin marker-based motion analysis? An in vivo study of knee kinematics[END_REF][START_REF] Lucchetti | Skin movement artefact assessment and compensation in the estimation of knee-joint kinematics[END_REF][START_REF] Ma | Analysis of three-dimensional in vivo knee kinematics using dynamic magnetic resonance imaging[END_REF][START_REF] Suzuki | In vivo kinematic analysis of posterior-stabilized total knee arthroplasty for the valgus knee operated by the gapbalancing technique[END_REF] to study knee kinematics, some of them focusing on influence of the main knee ligaments [START_REF] Blankevoort | Recruitment of knee joint ligaments[END_REF][START_REF] Feeley | Isometry of medial collateral ligament reconstruction[END_REF]. The results varied greatly across studies highlighting the high inter-individual variability.

Several finite element (FE) knee models have also been developed in an effort to comprehensively investigate various aspects of knee mechanics including contact pressure under various loads [START_REF] Adouni | Computational biodynamics of human knee joint in gait: from muscle forces to cartilage stresses[END_REF][START_REF] Hao | A Finite Element 3D Model of In Vivo Human Knee Joint Based on MRI for the Tibiofemoral Joint Contact Analysis[END_REF][START_REF] Kiapour | Finite element model of the knee for investigation of injury mechanisms: development and validation[END_REF][START_REF] Liu | Redistribution of knee stress using laterally wedged insole intervention: finite element analysis of knee-ankle-foot complex[END_REF][START_REF] Papaioannou | Patientspecific knee joint finite element model validation with high-accuracy kinematics from biplane dynamic Roentgen stereogrammetric analysis[END_REF][START_REF] Rao | A statistical finite element model of the knee accounting for shape and alignment variability[END_REF][START_REF] Wang | Comparison of stress on knee cartilage during kneeling and standing using finite element models[END_REF], ligament stress [START_REF] Limbert | Three-dimensional finite element modelling of the human ACL: simulation of passive knee flexion with a stressed and stress-free ACL[END_REF][START_REF] Ramaniraka | Biomechanical evaluation of intraarticular and extra-articular procedures in anterior cruciate ligament reconstruction: a finite element analysis[END_REF][START_REF] Xie | A study on construction three-dimensional nonlinear finite element model and stress distribution analysis of anterior cruciate ligament[END_REF], or the dynamic behavior [START_REF] Kiapour | Finite element model of the knee for investigation of injury mechanisms: development and validation[END_REF][START_REF] Abdel-Rahman | Three-dimensional dynamic behaviour of the human knee joint under impact loading[END_REF][START_REF] Beillas | A new method to investigate in vivo knee behavior using a finite element model of the lower limb[END_REF]. A few of them studied knee kinematics [START_REF] Adouni | Computational biodynamics of human knee joint in gait: from muscle forces to cartilage stresses[END_REF][START_REF] Kiapour | Finite element model of the knee for investigation of injury mechanisms: development and validation[END_REF][START_REF] Baldwin | Dynamic finite element knee simulation for evaluation of knee replacement mechanics[END_REF][START_REF] Mesfar | Biomechanics of the knee joint in flexion under various quadriceps forces[END_REF]. The high numerical costs associated with the computation of the biomechanical response of the knee joint and the full validation of these models against tibio-femoral and tibio-patellar kinematic data, however, are still a main issue [START_REF] Kiapour | Finite element model of the knee for investigation of injury mechanisms: development and validation[END_REF], particularly for patellofemoral motion because of uncertainty in patellar tracking.

A new experimental setup was recently developed combining 3D reconstruction imaging with the use of a motion capture system for accurate analysis of knee joint kinematics at our laboratory [START_REF] Azmy | EOS ® orthopaedic imaging system to study patellofemoral kinematics: assessment of uncertainty[END_REF]. The aim of the present study was to develop an FE model capable of reproducing the kinematic of the knee in flexion/extension and, building upon the work of Azmy, to validate it with in vitro experimentations for investigating the impact of ligament properties on the knee kinematics.

Materials and Methods

Finite Element Model

Geometry The model developed in this study is an adaptation of the one proposed by [START_REF] Noailles | Modelisation Elements Finis De L'articulation Du Genou[END_REF]. The 3D geometry of the bony structures was acquired using a sensor pen (Fastrak system, Polhemus, Olchester, USA), and was carried out on a left 63-yearold female lower limb.

Elements and Meshing

The FE mesh, represented in Fig. 2.1, consists of the distal and proximal part of the femur and tibia, respectively, the patella, the quadriceps tendon, the cartilage, and knee ligaments. Bones and cartilage were meshed with low-order (four-nodes) shell elements, the quadriceps tendon with loworder membrane elements, and the ligaments with tension-only cables.

The different ligament bundles were represented: (1) two bundles for each cruciate ligament (antéro-médial (AM) and postéro-lateral (PL) bundles for the Anterior cruciate ligament (ACL); antéro-latéral (AL) and postéro-médial (PM) bundles for the Posterior cruciate ligament (PCL)), (2) four bundles for the Collateral Ligaments (CoL) (three for the medial collateral ligament (MCL) and one for the lateral collateral ligament (LCL)), (3) five bundles for the patellar ligament, and (4) seven for the articular capsule.

Material Properties

The material properties used for each component are summarized in Tables 2.1 and 2.2. Bones and cartilage were defined as linear elastic isotropic material in accordance with the literature [START_REF] Choi | The elastic moduli of human subchondral, trabecular, and cortical bone tissue and the size-dependency of cortical bone modulus[END_REF][START_REF] Hori | Indentation tests of human articular cartilage[END_REF]. The regions covered with cartilage were modeled as a homogeneous bone cartilage material, with average properties. A thin strip of elements between the bones and cartilage regions were also defined with intermediate properties in order to avoid important mechanical discontinuity.

The stiffness of the different ligament bundles was chosen in accordance with the experimental data reported in literature [START_REF] Amis | Biomechanics of the PCL and related structures: posterolateral, posteromedial and meniscofemoral ligaments[END_REF][START_REF] Butler | Comparison of material properties in fascicle-bone units from human patellar tendon and knee ligaments[END_REF][START_REF] Chandrashekar | Sex-based differences in the tensile properties of the human anterior cruciate ligament[END_REF][START_REF] Race | The mechanical properties of the two bundles of the human posterior cruciate ligament[END_REF][START_REF] Robinson | Structural properties of the medial collateral ligament complex of the human knee[END_REF][START_REF] Sugita | Anatomic and biomechanical study of the lateral collateral and popliteofibular ligaments[END_REF]. Negative initial strain for PCL indicates a slack condition.

Contact Four frictionless surface-to-surface contact zones were created: femur (lateral)/tibia (lateral), femur (medial)/tibia (medial), femur/patella, and femur/quadriceps tendon. The interactions between the ligaments and bones represented in the model are not modeled in the present study.

Boundary Conditions

The experimental boundary conditions [START_REF] Azmy | EOS ® orthopaedic imaging system to study patellofemoral kinematics: assessment of uncertainty[END_REF] were imposed: the femur was fixed, the rope and pulley system was represented as two cables driven by temperature, and a set of 130 consecutive displacements of 0.5 mm were applied to the quadriceps tendon. This allowed to generate forces in the cables that were always directed toward the center of the femoral head.

Solution

The large deformation static response was computed using an implicit solver in ANSYS.

Model Evaluation

Sample Preparation Twenty-three fresh frozen lower limbs were used. The subjects from whom they were harvested were 47-97 years old. Absence of osteoarthritis and ligament laxity was checked. Specimens were disarticulated at the head of the femur and at the distal epiphysis of the tibia. All the soft tissues were removed except for the articular capsule of the knee, the quadriceps tendon, the CoL of the knee (medial and lateral), the patellar ligament, and the proximal and distal ligaments between the tibia and the fibula. The samples were frozen at 20 ı C then thawed at room temperature 12 h before experimentation.

Test Bench

The validated test bench, represented in Fig. 2.2, was adapted from [START_REF] Azmy | EOS ® orthopaedic imaging system to study patellofemoral kinematics: assessment of uncertainty[END_REF].

The femur was fixed and the tibia was free. Tension was applied to the quadriceps tendon using a weight of 9.8 N. A motor was connected to the tibial pilon by a rope and was used to perform the flexion movement. A progressive displacement was applied at the centroid of the malleolus and directed toward the center of the femoral head.

Movement Tracking Tripods were fixed on the femoral diaphysis, the proximal epiphysis, and the patella. During the cycles of flexion/extension, the position of the bony segments was tracked with a Polaris optoelectronic system (NDI, Waterloo, Ontario, Canada).

Anatomical Frames

The position of the tripods in their respective anatomical frames was calculated from 3D reconstruction made with the EOS ® bi-planar X-ray system. The experimental setup was put in the EOS cabin and bi-planar radiographies were acquired. A geometric reconstruction of the bony structures was performed and used for the calculation of the anatomical frames.

Movement Analysis

The variation of relative position tibia/femur and patella/ femur was extrapolated using a MATLAB ® routine. The rotations were calculated on mobile axis, with the following sequence (center of rotations): zy 0 x 00 (y 0 denoting the mobile y axis after the first rotation about z, and x 00 denoting the mobile x axis after the first two rotations).

Sensitivity Analysis

A sensitivity analysis was conducted to investigate the impact of ligaments on the knee kinematics. The stiffness and initial strain of one or several bundles of ligament were modified and the impact in the tibia position was computed. Ten cases were considered. These are summarized in Table 2.3. Stiffness values were defined based on data reported in the literature [START_REF] Amis | Biomechanics of the PCL and related structures: posterolateral, posteromedial and meniscofemoral ligaments[END_REF][START_REF] Butler | Comparison of material properties in fascicle-bone units from human patellar tendon and knee ligaments[END_REF][START_REF] Chandrashekar | Sex-based differences in the tensile properties of the human anterior cruciate ligament[END_REF][START_REF] Race | The mechanical properties of the two bundles of the human posterior cruciate ligament[END_REF][START_REF] Robinson | Structural properties of the medial collateral ligament complex of the human knee[END_REF][START_REF] Sugita | Anatomic and biomechanical study of the lateral collateral and popliteofibular ligaments[END_REF]. Initial strain values were chosen so that strains remained below the physiological limit whilst being sufficiently different from one another.

Results

Experimental/Numerical Comparison

Tibial Kinematics Both the experimental corridor of the tibial kinematics with respect to the femur and the numerical results computed with our model are shown in Fig. 2.3. Our results show that the angular position of the tibia of the finite element model is in the physiological corridor. During flexion, the tibial movement predicted by the FE model does an internal rotation (Fig. 2.3) of up to 11.9 ı at 70 ı flexion (while the in vitro experiment is in the corridor [ 18.3 ı ; 7.2 ı ]). Our results also

Influence on the Femoro-Tibial Kinematics

The influence of the ligament parameters on the femoro-tibial kinematics is reported in Table 2.4, mean value and standard deviation of tibia position at 60 ı flexion for each case. Concerning the Rx rotation (abduction/adduction), the ligament that had the most influence was the ACL, followed by the LLE, the PCL, and the MCL holding a smaller influence (SD < 0.35 ı ). The Ry rotation is mainly influenced by the PCL and the LCL, the PL bundle of ACL holding a smaller influence (SD < 1.01 ı ). The Tx translation was influenced mainly by the ACL and the PCL when both bundles are modified at the same time. The Ty translation was influenced mainly by the AM bundle of the MCL and the cruciate ligaments when the four ligaments are modified at the same time. The Tz translation is influenced mainly by the ACL, the AL bundle of the PCL, and the LCL. Despite the numerous contact elements in the model, the computation time was below 1 h. Model-predicted response was within the experimental corridors for all translations and rotations of tibia and patella with respect to the femur, except for a ˙2ı max discrepancy in the abduction/adduction rotation of the patella. The different ligament parameters appeared to have little effect on the patellar anterior/posterior and inferior/posterior shifts while the kinematic response of the tibial anterior/posterior shift was more sensitive.

Discussion

Very few FE models reported in the literature feature extensive validation against both tibio-femoral and tibio-patellar kinematic data. These are paramount to comprehensively assess the biomechanical role of the knee joint and to investigate various aspects of knee instability such as knee ligament injuries, ligament replacement, and ACL graft design.

In this contribution, an FE of the knee joint capable of reproducing the kinematics of the knee in flexion/extension for the investigation of the role of ligaments in the joint kinematic behavior was created and validated against the experimental data of 23 cadavers. Special attention was paid to keep the numerical cost as low as possible. Low-order shell elements and cable elements were therefore used to discretize the components of the geometric model.

For both the tibial and patellar kinematics, the results computed with our model are inside or at the limit of the experimental corridors. Both the experimental and computational results obtained for the femoro-tibial and the femoro-patellar kinematics are comparable to the literature. Tables 2.5 and 2.6 summarize and compare the kinematic values reported in this study with the main values reported in the literature.

As for model evaluation, Kiapour et al. [START_REF] Kiapour | Finite element model of the knee for investigation of injury mechanisms: development and validation[END_REF] presented extensive validation using 16 lower limbs and investigation of several parameters. However, patellofemoral kinematics was not evaluated.

[5] 8 1.5 ˙3 3 ˙5.5 - - - [37] 6 1 ˙5.5 1 ˙10 - - - [4] 7 0 ˙4 2 ˙5 - - -
The sensitivity analysis conducted on the main ligaments of the tibio-femoral articulation highlighted the importance of distinguishing the two bundles of the ACL and the PCL, since they do not impact the femoro-tibial kinematics in the same way: the AM bundle of the ACL exerts more than twice as much influence on Ry and Tx, but less than twice as much on Rx and Ty as compared with the PL bundle. Likewise, modifying both bundles of the PCL significantly increases the influence on Rx, Ry, and Tx, but doesn't have much effect on Ty and Tz.

Two main limitations of our model are the omission of (1) the meniscus, the popliteal tendon, and the popliteo-fibular ligament, in order to keep the computation time as low as possible. This might explain the small experimental-numerical kinematic discrepancies and could be taken into account to improve the realism of the model. The simplification of the ligament insertion and the absence of ligament/ligament contact could also explain the fact that, in the simulation results, the ligament effect remains small. (2) The geometry of the finite element model is not personalized. Yet the geometry has a strong impact on the femoro-tibial and the femoro-patellar kinematics for at least two reasons: it affects the morphology of the contact surface, and it can determine the position of anatomical landmarks from which the boundary conditions are defined (center of the femoral head, center of the malleoli).

Nonetheless, the model appears consistent with in vitro kinematics and the computational cost remains reasonable (1 h). From a clinical perspective, the developed knee model constitutes a valuable tool for exploring new treatment strategies, such as the development of ACL substitutes for ligament reconstructions where it can be used to test different designs and evaluate the restoration of knee kinematics right after implantation.

Conclusions

An FE of the knee joint capable of reproducing the kinematic of the knee in flexion/extension for the investigation of the role of ligaments in the joint kinematic behavior was created and validated against the experimental data of 23 cadavers. Special attention was paid to keep the numerical cost as low as possible. For both the tibial and patellar kinematics, the results computed with our model are inside the experimental corridor. Both the experimental and computational results obtained for the femoro-tibial and the femoro-patellar kinematics are also comparable to data reported in the literature. From a clinical perspective, the developed knee model constitutes a valuable tool for exploring new treatment strategies, such as the development of ACL substitutes for ligament reconstructions where it can be used to test different designs and evaluate the restoration of knee kinematics right after implantation.

Fig. 2 . 1

 21 Fig. 2.1 FE knee model. (a) Isometric view (front); (b) isomeric view (back). The femur and tibia were limited to their respective distal and proximal parts

Fig. 2 . 2

 22 Fig. 2.2 Validated experimental setup used for testing the lower limbs in flexion/extension. (a) Whole setup with the rope and pulley system for applying the flexion. (b) Position of the tripods

Fig. 2 . 3

 23 Fig. 2.3 Tibial kinematics (with respect to the femur) during the range of motion

Fig. 2 . 4 Table 2 . 4

 2424 Fig. 2.4 Patellar kinematics (with respect to the femur)Table 2.4 Mean value and standard deviation of tibia position at 60 ı flexion for each sensitivity case (the ligament column shows which ligament was modified) Ligament Rx ( ı ) Ry ( ı ) Tx (mm) Ty (mm) Tz (mm) Mean 1.2 13.9 17.3 5.9 3.0 ACL (AM) SD 1.0 0.5 0.5 0.5 0.2 Mean 0.5 15.5 16.8 5.0 3.0 ACL (PL) SD 0.6 1.0 0.9 0.2 0.2 Mean 1.3 14.2 17.2 5.9 2.9 ACL (all) SD 1.0 0.3 0.4 0.6 0.2 Mean 0.0 15.6 17.9 5.1 3.1 PCL (AL) SD 0.1 1.8 0.3 0.1 0.1 Mean 0.1 15.2 17.7 5.1 3.2 PCL (PM) SD 0.3 1.4 0.2 0.3 0.0 Mean 0.3 14.8 17.2 5.0 3.2 PCL (all) SD 0.3 3.0 0.8 0.4 0.1 Mean 0.2 14.2 17.7 5.2 3.2 MCL (all) SD 0.2 0.3 0.1 0.1 0.1 Mean 0.2 14.6 17.6 5.1 3.2 LCL SD 0.7 1.7 0.2 0.1 0.2 Mean 0.2 14.3 17. 5.1 3.0 CoL (all) SD 0.7 1.3 0.1 0.1 0.1

Table 2 .1 Material properties of bones and cartilage

 2

		E (MPa)
	Cortical bone and patella 12,000	0.3
	Bone/cartilage inter	2000	0.4
	Bone/cartilage	250	0.4
	Tendon	90	0.4

Table 2 . 2

 22 Material properties of the different ligament bundles

		K (N/mm) Initial strain (%)
	Anterior cruciate ligament (ACL) 75	5
	Posterior cruciate ligament (PCL) 75	3
	MCL	70	0
	LCL	20	0

Table 2 . 3

 23 Stiffness and initial strain chosen for the sensitivity cases on the ACL

		K (N/mm)	Initial strain (%)
		Min Mean Max Min (%) Mean (%) Max (%)
	ACL bundle 40	75	150 0	3	6
	PCL bundle 40	75	150 3	0	4
	MCL	70	100	120 0	2	4
	LCL	20	60	120 0	2	4

Table 2 . 5

 25 Comparison of the experimental and computational results with literature: kinematics of the femoro-tibial articulation at 60 ı flexion

					.9	0.2	0.2
	ACL (all)	Mean	1.3	14.2	17.2	5.9	2.9
		SD	1.0	0.3	0.4	0.6	0.2
	PCL (AL)	Mean	0.0	15.6	17.9	5.1	3.1
		SD	0.1	1.8	0.3	0.1	0.1
	PCL (PM)	Mean	0.1	15.2	17.7	5.1	3.2
		SD	0.3	1.4	0.2	0.3	0.0
	PCL (all)	Mean	0.3	14.8	17.2	5.0	3.2
		SD	0.3	3.0	0.8	0.4	0.1
	MCL (all)	Mean	0.2	14.2	17.7	5.2	3.2
		SD	0.2	0.3	0.1	0.1	0.1
	LCL	Mean	0.2	14.6	17.6	5.1	3.2
		SD	0.7	1.7	0.2	0.1	0.2
	CoL (all)	Mean	0.2	14.3	17.	5.1	3.0
		SD	0.7	1.3	0.1	0.1	0.1

Table 2 . 6

 26 Comparison of the experimental and computational results with literature: kinematics of the femoro-patellar articulation at 60 ı flexion

	Study	Specimens Rx ( ı )	Ry ( ı )	Tx (mm)	Ty (mm) Tz (mm)
	FE model	-	0	0.5	-41.5	-20	-34
	In vitro exp 23	3 ˙3	3.5 ˙4.5	41.5 ˙3.5	17 ˙3	41	˙4
	[2]	8	3 ˙3 0.5 ˙7.5	38 ˙4	14 ˙4	40	˙8


Jordan Gauvrit

Guillaume Ansanay-Alex

Comparative analysis of model reduction methods applied to building simulation benchmarks

Keywords: thermal model, model order reduction, computation time, spectral analysis

Many research and development studies use virtual buildings represented by building thermal models at a detailed level. Such models need relatively long computation times for one-year simulations: the use of these tools is then often too time-consuming when it comes to compare technical solutions during the design stage. Model reduction methods allow models with shorter computation times to be synthesized. In this paper, the linear state-space representation of a whole model is extracted and a balanced truncation method is applied to it. The detailed models are built from the SIMBAD library, a Simulink library of building modeling components developed at CSTB. Both linear state-space and reduced order models ensure shorter computation times than the full detailed model. However, the choice of the order of the reduced model has an impact on the final results. The main strengths and weaknesses of using the linear state-space and reduced models built from the same detailed model are investigated. To this end, the physical descriptions of idealized test buildings provided by ASHRAE standard 140 for building simulation tools assessment are used, and the results obtained through the different approaches for computation time reduction are compared.

INTRODUCTION

In 2012, commercial and residential buildings accounted for 44% of the final energy consumption in France and were responsible for 21% of Greenhouse Gas Emissions [START_REF] Ademe | Climat, air et énergie. Chiffres-clés[END_REF]. Among this consumption, heating and cooling systems stand for a large share (59%, [START_REF] Ademe | Climat, air et énergie. Chiffres-clés[END_REF]). Therefore, reducing energy consumption requires the improvement of HVAC (Heating, Ventilation, Air Conditioning) systems efficiency. There are two ways to reach this goal. First, replacing old systems with high efficiency systems increases energy performance (Visier, 2008). Secondly, control methods allow to adapt the system operation to match occupants' needs or renewable energy production ((Visier, 2008),("Energy future, Think Efficiency," 2010)) and reduce energy waste.

Modelling innovative systems, complex distribution networks or control strategies for design and optimization purposes require detailed models and often short time steps. This can lead to long computations, hardly affordable during the design stage or for embedded control models.

The SIMBAD (SIMulator of Building And Devices) library, developed at CSTB [START_REF] Husaunndee | SIMBAD: A simulation toolbox for the design and test of HVAC control systems[END_REF], is composed of a set of Matlab/Simulink components enabling user to build dynamic models of buildings, their systems and especially their control strategies. The components of this library will be used in this paper.

Simulating buildings with this level of detail leads to much longer computation times than static hourly models. In SIMBAD's case, an annual simulation of a mono-zone building takes approximately 130s with a fixed time step Δ 20s. Consequently, and considering even longer simulations when all production, distribution and emission systems are modelled, using such detailed models during the design stage may be difficult. It could then be useful to implement methods which reduce computation time while keeping sufficiently accurate results. The building envelope model considered in this paper will not take into account air transfer or non-linearized radiative transfers, but still shows a nonlinear behaviour regarding long wave radiations through windows.

The concern in this paper is to evaluate the performances of reduced models. On the one hand, the computation time saved with a reduced order model, and on the other hand, the accuracy difference between the results of detailed, linearized and reduced order models. Section 5 is dedicated to the problem of evaluating how well linearized or reduced order models perform with respect to the original detailed model.

In order to evaluate the detailed and reduced order building models on reference cases, the reference buildings used are those described by ASHRAE Standard 140-2001 ("ANSI/ASHRAE Standard 140-2001, Standard method of test for the evaluation of building energy analysis computer programs," 2001). It defines test procedures for evaluating capabilities of building energy programs. In fact, it gives a physical description of different buildings and the predictions from well-known programs (such as TRNSYS, DOE-2 or S3PAS) applied to them. Some of these test cases were implemented to compare synthesized models with the SIMBAD detailed model. Section 3 will introduce the SIMBAD models used as the detailed modelling level. In section 4, the linearization and common model order reduction methods will be presented. Section 5 clarifies the mobilized benchmarks indicators. Benchmark results will be presented and discussed in section 6. Finally, a variant of the order reduction method will be tested to improve its performance in section 7.

DETAILED SIMBAD MODEL

The SIMBAD multi-zone building model is a white box model where each component (zones, walls, windows …) is described by a set of physical parameters. It includes a dynamic model of each wall surface and zone air temperature. SIMBAD models are built using Simulink. Each component has inputs and outputs that are connected to create the whole model. Two cases are set up. The first one is the free floating case (indoor building temperature is not controlled, Figure 1). The second one is the regulated case (indoor building temperature is regulated using a PID controller, Figure 2). : Heating convective power (W). Positive for heating, negative for cooling.

Multizone SIMBAD model

Solar processor

The solar processor block projects solar flows on surfaces of the building envelope. This block uses many trigonometric functions and is as such strongly nonlinear. The multi-zone block calculates indoor temperature by evaluating surface temperatures taking radiative, convective and conductive heat exchanges into account.

The PID controller computes the heating convective power in order to reach the setpoint value. It is strongly nonlinear and won't be subject to model order reduction.

In the rest of this paper, the detailed model will be represented by the letter D. This technique works by nature without any loss of precision when blocks are linear.

COMPUTATION TIME AND MODEL ORDER REDUCTION TECHNIQUES

Linearization

There are two ways to apply this technique to the detailed building model:  The first one is the linearization of the multi-zone building block only. This block is then replaced with a SSM block but the solar processor is maintained. The resulting SSM block has 23 inputs.  The second one is the linearization of the assembly between the solar processor and the multi-zone block. As the solar processor is not linear, a simplification is made to make it linear: the sun position is set to a constant and becomes a model parameter.

Then, the assembly is replaced with a SSM block. Solar radiations are adjusted to fit the initial solar radiation projected on the window. This SSM block requires 8 inputs.

Model order reduction

After projecting (1) on a particular basis, the system is truncated (2).

𝑇 ̃̇( ) 𝐴 ̃𝑇 ̃( ) 𝐵 ̃ ( )

̃ 𝐶 ̃𝑇 ̃( ) (2)
In the following parts, the notation will stand for the outputs of linearized models as well as the outputs of the reduced order model.

The order of (2) has decreased from (1). Several methods performing this reduction are available . One of the first method was the balanced truncation and appeared in the early 80's. It was described by B.C. Moore in [START_REF] Moore | Principal component analysis in linear systems: Controllability, observability, and model reduction[END_REF]. Other methods followed in the same decade as the Hankel-norm reduction of Glover (Curtain & Glover, 1985) or whether the proper orthogonal decomposition described by Sirovich [START_REF] Sirovich | Turbulence and the dynamics of coherent structures[END_REF]. Finally, in the 90's, two other methods came out. The first one was explained by P. Feldmann in [START_REF] Feldmann | Efficient linear circuit analysis by Padé approximation via the Lanczos process[END_REF] and was related to Krylov subspaces and Padé approximation. The second one is called PRIMA and was the subject to the paper [START_REF] Odabasioglu | PRIMA: passive reduced-order interconnect macromodeling algorithm[END_REF]. In this study, the balanced truncation is used. It has been applied in thermal building modelling since the mid-1990s. Dautin's works (Dautin, 1997) or Déqué's works [START_REF] Déqué | Réduction d'un système linéaire et invariant par la technique de Moore. Application à la thermique du bâtiment[END_REF]) can be quoted. More recently, other studies like [START_REF] Kim | Reduced-order building modeling for application to modelbased predictive control[END_REF]) used this technique.

The main idea is to delete states that have a small sensitivity to inputs (controllability) and a small impact on outputs (observability). The controllability gramian (resp. observability gramian ) characterizes the state's controllability (resp. observability). It is shown [START_REF] Bruq | Théorie du signal: modélisation statistique, automatique et traitement[END_REF]) that they can be calculated by Lyapunov equations (3):

𝐴 𝐴 𝐵𝐵 𝐴 𝐴 𝐶 𝐶 (3)
and depend on the basis where (1) is represented. It should be noted that a state can be sensitive to inputs and have a small impact on outputs at the same time. States thus cannot be deleted by considering the gramians only. Using a balanced representation overcomes this obstacle.

It consists in applying a transform 𝑇( ) 𝑃 ( ) to (1) such that the gramians of system (4) are equal.

̇( ) 𝑃 𝐴𝑃 ( ) 𝑃 𝐵 ( )

( ) 𝐶𝑃 ( ) (4)
Palomo [START_REF] Palomo | Résolution de problèmes thermiques de grande dimension[END_REF] gives steps to establish 𝑃 :

-1 st step: Computation of and with (3) -2 nd step: Calculation of with Cholesky decomposition :

(5)

where is a lower triangular matrix.

-3 rd step: Diagonalization of (symmetric, positive definite) obtained by solving the following symmetric eigenvalues/eigenvectors problem:

(6)
where is a diagonal matrix and is such that .

-4 th step: Calculation of 𝑃 :

𝑃 / (7)
Then, (4) can be reduced by deleting the less influential states, leading to (2). The reduction order is the number of remaining states. This model will be called (R) for Reduced order model. A reduced order model of order 5 will be denoted R5.

BENCHMARKING

Computation time notation

The reference time is obtained by the simulation time of D. As regards L and Rx models, the time spent to reduce the detailed model is taken into account.

An a priori evaluation

The "quality" of R can be evaluated by analyzing its state space representation without the aid of a simulation. [START_REF] Palomo | Résolution de problèmes thermiques de grande dimension[END_REF], [START_REF] Palomo | Guidance for the selection of a reduction technique for thermal models[END_REF]) defines an index of the spectral quality called ( ). For each output , it gives the relative error (between L and R) according to input frequency .

To achieve this target, transfer function matrix of the linearized (resp. reduced order) model ( ) (resp. ̃( )) is calculated by :

( ) 𝐶( 𝐴) 𝐵 ̃( ) 𝐶 ̃( 𝐴 ̃) 𝐵 ̃ (8)
Then, ( ) is determined by:

( ) ∑ | ( ) ̃ ( )| ∑ | ( )| (9)
With:

: Number of inputs ( ) (resp. ̃ ( )) : Transfer function of the th output from the th input. It belongs to the matrix ( ) (resp. ̃( )).

Thus, the reduction order can be chosen by setting a risk tolerance threshold on the index.

An evaluation by simulation

Also, L and Rs have to be compared with the detailed model. A one year simulation is performed for each model and the results are studied. The gap between detailed model outputs and reduced model outputs, written ( ), is called the reduction error.

Statistical analysis

A first statistical analysis of ( ) by calculating the mean and the standard deviation will be made. Moreover, heating power and cooling power results will be integrated:

-over one hour time steps to compare maximum peaks (ASHRAE standard indicator) -over a year to compare annual consumptions

Spectral analysis

Secondly, a spectral analysis of ( ) will be performed. It is especially described for building thermal modeling in (Dautin, 1997) or [START_REF] Ramdani | How to improve building thermal simulation programs by use of spectral analysis[END_REF]. It aims at evaluating the dynamic behavior of reduced models in comparison with the detailed model.

The power spectral density (PSD) of ( ), noted ( ) and calculated by ( 10), shows how the energy of the signal is distributed over the different frequencies. It must not be confused with the energy of the physical quantity.

( ) | ( )| (10)
where ( ) is the Fourier transform of ( ).

The links between reduction error and model inputs are also investigated. Spectral coherence of ( ) and an input signal ( ) gives information about these. It is noted ( ) and calculated by (11). It expresses the linear correlation between both signals according to the frequency.

( ) | ( )| ( ) ( ) (11)
Where ( ) is the cross spectral density (CSD) of both signals and is determined by :

( ) ( ) ̂( ) (12)
Coherence allows comparing the influence of each input on reduction error of one reduced model. However, it is not a sufficient tool for comparing the influence of one input over two models. The partial coherence method (PCM) is an answer to this problem.

PCM offers a decomposition of ( ) such as:

( ) 0. Now, if the contribution of input 2 is calculated at first, then ( ) ( ) and ( ) 0. However, the unexplained portion stays the same regardless of the order of inputs. For the following parts, the order of inputs will be fixed for the comparison of models. Finally, results of the spectral analysis will be summed up by frequency bands. They are defined by the observation of the PSD of a signal. Here, (Figure 4), the signal is the error on the indoor temperature prediction of L. Four frequency bands are chosen:

-VLF band: Very low frequency band (between 0 to (48h) -1 ). It describes quasi-static behaviour -LF band: low frequency band (between (48h) -1 to (16h) -1 ). It describes slow dynamic behaviour.

-HF band: high frequency band (between (16h) -1 to (9h) -1 ). It describes medium dynamic behaviour.

-VHF band: very high frequency band (superior to (9h) -1 ). It describes fast dynamic behaviour.

The portion of the energy of the signal ( ) contained in the frequency band 𝐵 is noted

, and is calculated by:

, ∫ ( ) ∑ ∫ ( ) (16)
The part of the energy of the signal ( ) explained by input , contained in the frequency band 𝐵 , is noted , and is calculated by:

The portion which cannot be explained by the inputs is called the unexplained part. A graphic representation of these quantities is available in section 6.2.3.

RESULTS

Hypotheses

Test cases

In this paper, ASHRAE test cases used are case600FF, case900FF, case600 and case900. They reflect a mono-zone building with two south-oriented windows (Figure 5). Simulations are conducted over a one year period with a fixed time step Δ 20s. According to the case, envelope inertia and indoor temperature control change (Table 1). In particular:

-Free floating case: indoor temperature is not controlled -Regulation case: indoor temperature is regulated. The heating set-point is 20°C and the cooling set-point is 28°C.

Choice of reduced model orders

The spectral index, defined in subsection 5.2, is a useful tool to anticipate the quality of the results of Rs according to reduction order. It should be reminded that the spectral index is calculated by comparing L and Rs. Figure 6 and Figure 7 show indices estimated for case600 and case900 with SSM block including 23inputs or 8inputs. On each chart, a gap of the index may be observed between order5 and order6. Then, R6 and L make very close predictions. In the following parts, the results of R2, R4 and R6 will be exposed.

Computation times

Computation times are shown in Table 2. L with 23 inputs (resp. 8 inputs) roughly divides by two (resp. by three) the simulation time. In fact, for D, Matlab computes the outputs of each subcomponent at each time step: it considers the detailed model as a time variant system. In L's case, time invariance property is ensured through the linearization stage. That's why this computation is simplified and faster.

Secondly, it can be noticed that 8-input models have shorter computation times than 23-input models. Indeed, the solar processor block has been linearized and has been "included" in the SSM block. Moreover, a SSM with 8 inputs is smaller than a SSM with 23 inputs. Finally, Rs have results very close to L. It seems to be a paradox because matrices of Rs are smaller than L's ones. In fact, matrices of L are sparse (Figure 8) and Matlab optimizes the solving of SSM according to the number of null coefficient inside matrices. Then, even if sizes of L matrices are bigger, their sparsity must be taken into account to understand computation times.

Statistical analysis

Differences between indoor temperature prediction of D with L and Rs have been calculated. Means and standard deviation of the reduction errors are shown in Table 3 andTable 4.

As could be expected, L and R6 give very similar results. The mean of the difference between their results is close to 10 -11 °C.

Table 3: Mean of reduction error (°C) -Free floating cases

Case L R2 R4 R6 Case600FF/8inputs -3.17e-02 -2.77e+00 -2.66e-02 -3.17e-02 Case600FF/23inputs -8.00e-06 -3.82e+00 -2.27e-03 -8.00e-06 Case900FF/8inputs -2.90e-02 -9.91e-01 -7.51e-01 -2.90e-02 Case900FF/23inputs -5.12e-06 1.31e+00 -4.24e-01 -5.12e-06

Results of L with 8 and 23 inputs are close to D. However, R2s give results far from D. Models with 23 inputs give generally better results than model with 8 inputs. This probably stems from an exaggerated simplification of projected solar radiations made in the 8-input model.

Table 4

: Standard deviation of reduction error (°C) -Free floating cases

Case L R2 R4 R6 Case600FF/8inputs 5.09e-01 1.49e+00 5.09e-01 5.09-01 Case600FF/23inputs 7.27e-03 2.22e+00 7.24e-03 7.27-03 Case900FF/8inputs 3.60e-01 3.89e-01 3.34e-01 3.60-01 Case900FF/23inputs 6.03e-04 5.33e-01 1.45e-01 6.03-04

The same observations as before can be made about dispersion of reduction errors, in terms of difference between 8-input models and 23-input models, and the evolution of dispersion according to reduction order.

Spectral analysis

For each model, the decomposition of the reduction error was calculated thanks to (16) (energy of the residual in a frequency band) and ( 17) (energy explained by each specific input in a frequency band). The inputs are meteorological variables (outdoor temperature, sky temperature, diffuse and direct solar radiations). Table 5 (resp. Table 6) shows results for case600FF (resp. case 900FF). In these graphs "Inex" denotes the energy unexplained by inputs.

Less than 10% of the reduction error is not explained by the chosen inputs. Moreover, the reduction error is mainly connected with static and slow dynamic behaviors (more than 80% of the signal's energy, i.e. the error). Then, climatic data are relevant to explain the error. The decomposition of R2 looks similar for both models. The effect of radiative transfer modeling is then pretty close and they produce equivalent dynamic behaviors.

It can be noticed that the VLF part of the reduction error decreases with the reduction order being higher. Concerning case900FF, reduction error is more connected with static behavior than in case600FF. That shows the effect of building inertia on the reduction error.

Also, an important effect of inputs linked to radiative transfer on models with 8 inputs can be observed, which points towards the simplification made to the solar processing part.

Computation times

Computation times of regulation cases are shown in Table 7. The same observations as the previous case can be made about time savings and reduction methods.

Statistical analysis

Differences between output predictions of D and L and Rs have been calculated. Mean and standard deviation of the reduction error on indoor temperature are shown in Table 8 and Table 9.

Two behaviors can be distinguished in the regulation case. The first one is the period when indoor temperature is controlled (heating or cooling period). The temperature is quasi constant. The second one is a free floating period. When the indoor temperature is between 20°C and 28°C (the two set-points), indoor temperature is not controlled. In fact, periods when indoor temperature has a static behavior limit the value of reduction error. Then, the mean is often lower than free floating case. Mean reduction errors higher than in the free floating case are underlined in the following tables. 2.00e-02 -3.06e-01 2.07e-02 2.00e-02 Case600/23inputs 2.75e-04 -4.58e-01 -2.04e-05 2.75e-04 Case900/8inputs -1.80e-02 -2.62e-01 -1.78e-01 -1.80e-02 Case900/23inputs 4.89e-05 2.60e-01 -7.04e-02 4.89e-05 About the standard deviation, the same observations can be made about models with 8 inputs. However, in this case, the reduction error of models with 23 inputs is more scattered than free floating case.

Table 9: Standard deviation of reduction error on indoor temperature (°C) -Regulation cases

Case

L R2 R4 R6 Case600/8inputs 1.51e-01 5.76e-01 1.56e-01 1.51e-01 Case600/23inputs 4.38e-02 8.14e-01 6.20e-02 4.38e-02 Case900/8inputs 1.42e-01 4.00e-01 2.45e-01 1.42e-01 Case900/23inputs 2.69e-02 3.39e-01 1.13e-01 2.69e-02 Finally, these statistical indicators show predictions of L and Rs are pretty close to those of D.

For confirming this trend, annual consumptions and hourly peaks of heating and cooling system have calculated. Results are presented in Table 10 andTable 11. In the background of each chart, the banner represents the score interval of the set of tools presented in the ASHRAE standard. The dotted line is the mean result of the standard. R2, with 23 or 8 inputs, gives results outside the banner. It can be noticed for the cooling and the heating consumptions in Case600. However, all other models make coherent predictions with the standard.

Finally, R4 seems to be a sufficient model to predict consumptions and hourly peaks close to those of D. In fact, a maximum deviation of 3.2% can be found between R4 and D results. It concerns cooling peak of a model with 8 inputs.

Table 10 : Consumptions and hourly peaks -Case600

Heating Cooling

Consumption

Hourly peak

Table 11 : Consumptions and hourly peaks results -Case900

Spectral analysis

The decomposition of indoor temperature reduction error has been calculated for each model. It is shown in Table 12 for case600 and in Table 13 for case900. The effect of heating convective power has been added to the analysis.

In case600, reduction error of L, R4 and R6 is mainly in the very high frequency band: this could be expected considering that the VHF band spans frequencies higher than (9h). Decompositions are close between the model with 8 inputs and the model with 23 inputs but the simplification on the solar processor in the 8-input model creates low frequency error. The

Heating Cooling

Consumption

Hourly peak unexplained part is quite important (more than 10% here). It can be noticed all the inputs have the same contribution on the reduction error (around 16%).

In case900, a higher order model ( R6) is necessary to get rid of the reduction errors showing up in the VLF, LF and HF frequency bands. The errors on solar heat gain implied by the simplification of the solar processor have much more impact than in the lighter case600.

In the following graphs, Pcl denotes the convective power input.

Description

With a sufficient order, balanced truncation gave results close to L. However, computation time does not decrease because, contrary to L, matrices of reduced order models are not sparse and the number of inputs is large. Recent methods such as SVDMOR and RecMOR [START_REF] Feldmann | Sparse and efficient reduced order modeling of linear subcircuits with large number of terminals[END_REF] or TermMerg [START_REF] Liu | Termmerg: an efficient terminal-reduction method for interconnect circuits[END_REF] are able to reduce model with a large number of inputs. In this paper, the ESVDMOR method will be used, a generalised version of SVDMOR, which was introduced by Liu & Tan [START_REF] Liu | An efficient terminal and model order reduction algorithm[END_REF].

There are two steps to perform this method. First, the number of terminals is reduced by finding some correlations between the terminals thanks to Singular Value Decomposition (SVD). Secondly, the order of the SSM is reduced with a common method of Model Order Reduction. In this paper, the balanced truncation will be used. Step2: Reduction of the order of the identified state space with balanced truncation method.

In the current case study, there is only one output. ESVDMOR will then just reduce the number of inputs. It will be performed on the SSM with 23 inputs.

Choice of the number of inputs

Liu & Tan suggest choosing the number of inputs in order to preserve significant singular value computed in step 1. However, this method was not adapted to the problem at hand. Indeed, with this criterion, only one input was preserved and the predictions of the reduced model were imprecise. So, in order to determine the number of inputs, the spectral index previously defined (Figure 9) will be used. In this figure, the balanced truncation is performed at the order 6.When the 23 inputs are preserved, the spectral index is the same as Figure 6 (case600) and Figure 7 (case900). Spectral index is quickly downgraded when the number of inputs decreases. However, it is still satisfying. To confirm that, in the following part, the results of ESVDMOR will be compared to those of simple balanced truncation. The main results of simulation with ESVDMOR are shown in Table 15. The model with 5 inputs gives results with the best compromise between computation time and results quality. However, the saving is not significant (less than 0,2s on around 50s). In fact, the step which allows computing from U and from ′is added in this method. So, the gain obtained thanks to the reduced number of inputs is compensated by this additional step.

CONCLUSION

In this study, two methods to reduce the order of a detailed building model were used. The first was the balanced truncation. If the order is sufficient, the reduced order model will make predictions close to the detailed model. Two indicators were used to qualify the models. The spectral index gives an a priori evaluation of the adequate reduction order. The spectral analysis of the residual between the results of detailed and reduced models gives a good qualitative insight of the contribution of inputs on the residual and is as such an indication of the submodels responsible for the gaps between model results. It should, though, be used with caution when it comes to quantitative evaluations, since the order in which inputs are considered matters. The time saving of model order reduction is not significant on the test cases used, compared to the time saving obtained by the extraction of a time-invariant statespace model. The sparsity of the matrix corresponding to the linearized detailed model, and the large number of inputs of the linearized system, explain this low time saving.

Then, a second method called ESVDMOR was tested. It permits to create a state space model with a lower number of inputs and perform balanced truncation on it. However, an additional step is required. It has to be taken into account because it is computed at each time step. The results are quickly downgraded according to the number of inputs but still acceptable. Nevertheless, this method does not allow saving much time.

  Figure 4: PSD signal example and frequency bands

Figure 5 :

 5 Figure 5: Building geometry of cases

Figure 6 :Figure 7 :

 67 Figure 6 : Spectral indices according to model input numbers -Case600

Figure 8 :

 8 Figure 8: Sparsity of matrices for L and R4

Table 5 :

 5 Decomposition of reduction error -Free floating, case600FF the decomposition of reduction error of L with 8 inputs shows an important effect of meteorological inputs linked to radiative transfers (with regard to decomposition of the model with 23 inputs). The same observation can be made for R4 and R6. Here shows up the effect of the simplification of the solar processor for the 8 inputs model.

  on System Simulation in Buildings, Liege, December 10-12, 2014

Figure 9 :

 9 Figure 9: Spectral indices according to the number of inputs -ESVDMOR

Table 1 :

 1 Cases' names

		Low mass High inertia
	Free floating case600FF	case900FF
	Regulation	case600	case900

Table 2 :

 2 Computation times -Free floating cases

	Case	D	L	R2	R4	R6
	Case600FF/8inputs 110,45s Case600FF/23inputs	29,71s 28,29s 28,61s 30,45s 55,97s 53,90s 55,63s 60,33s
	Case900FF/8inputs 108,41s Case900FF/23inputs	28,65s 26,26s 27,06s 27,70s 55,59s 53,62s 55,03s 57,05s

Table 6 :

 6 Decomposition of reduction error -Free floating, case900FF

	Case900FF/8inputs	Case900FF/23inputs
	L
	R2
	R4
	R6
	6.

3 Regulation cases

Table 7 :

 7 Computation times -Regulation cases

	Case	D	L	R2	R4	R6
	Case600/8inputs 138,67s Case600/23inputs	55,85s 52,83s 53,88s 57,22s 83,78s 84,67s 85,37s 87,53s
	Case900/8inputs 147,94s Case900/23inputs	55,59s 53,62s 55,03s 57,12s 86,46s 83,83s 86,59s 88,06s

Table 8 :

 8 Mean of reduction error on indoor temperature (°C) -Regulation cases

	Case	L	R2	R4	R6
	Case600/8inputs

Table 12 :

 12 Decomposition of indoor temperature reduction error -Regulation, case600

	Case600/8inputs	Case600/23inputs
	L
	R2
	R4
	R6

9th International Conference on System Simulation in Buildings, Liege, December 10-12, 2014

Table 13 :

 13 Decomposition of indoor temperature reduction error -Regulation, case900

	Case900/8inputs	Case900/23inputs
	L
	R2
	R4
	R6
	7.

A VARIANT OF BALANCED TRUNCATION

  Table 14 summarizes these steps. Plain black rectangles symbolise dense matrices and hatched rectangles denote sparse matrices.

Table 14 :

 14 ESVDMOR process

	ESVDMOR with balanced truncation

Table 15 :

 15 Results of ESVDMOR -Free floating cases

	Method	Computation time	Case600FF ( ) mean/st. dev. (°C)	Case900FF ( ) mean/st. dev. (°C)
	Moore	54,33s	-8.00e-06/7.27e-03 -5.12e-06/6.03e-04
	ESVDMOR (1input)	51,53s	1.43e-01/8.04e-01 1.67e-01/7.80e-01
	ESVDMOR (5inputs)	54,13s	-1.57e-01/4.97e-02 3.18e-01/3.30e-01
	ESVDMOR (10inputs)	55,15s	-1.18e-01/3.21e-02 4.23e-01/3.37e-01
	ESVDMOR (15inputs)	56,45s	7.33e-02/9.49e-02 4.36e-01/5.78e-01

9th International Conference on System Simulation in Buildings, Liege, December 10-12, 2014

ACKNOWLEDGEMENT

This research was supported by an Agence Nationale de la Recherche grant (Project MAEVIA # ANR-12-VBDU-0005).

With:

( ) : Linear contribution of the input to ( ) ( ) : Unexplained part of ( ) [START_REF] Ramdani | How to improve building thermal simulation programs by use of spectral analysis[END_REF] provides equation ( 14). Knowing the contribution of the input requires to know the contribution of the ( ) input. The whole contribution computation is then an iterative process. The computation scheme, drawn on the Figure 3, is established in the case of a model with 3 inputs.

With:

This pattern is generally applicable with any number of inputs. In addition, it should be noticed that the order in which inputs are sorted affects the computed values of contributions. For instance, considering two inputs which are severely correlated with ( ) ( ( ) and ( ) ). If the contribution of input 1 is calculated at first, then ( ) ( ) and