Théo Baumard
email: tbaumard01@qub.ac.uk

Gary Menary
email: g.menary@qub.ac.ukc

Olivier De Almeida
email: olivier.dealmeida@mines-albi.frd

Fabrice Schmidt
email: fabrice.schmidt@mines-albi.fre

Peter Martin
email: p.j.martin@qub.ac.uk

Jérôme Bikard

Towards a Coupled Heating-Forming Simulation of the Thermoforming of Thermoplastic Composites

A strategy for the simulation of the whole thermoforming process, from the infrared heating to the stamping, is presented here. Two loosely coupled simulation tools are being developed: the first one computes a realistic 3D, transient temperature field of the composite stack inside an infrared oven, considering the radiative, conductive and convective heat transfers; the temperature distribution is used as an input for the second that aims at simulating the thermomechanical behaviour of the composite during the forming step via a non-orthogonal constitutive model. The steps for the identification of the model parameters are introduced. Initial validation tests show realistic results in term of shear angle distribution.

INTRODUCTION

The use of thermoplastic composites in the automobile sector was up until recently limited to non-structural parts or low volume, high-end vehicles. With new CO 2 emissions regulations coming into effect by 2020 [START_REF]EU CO2 Emission standards for passenger cars and light commercial vehicles[END_REF], their lightweightness and shorter manufacturing time achievable compared to thermosetting composites makes them interesting as possible replacement for metal in high volume applications. In order to efficiently reduce the costs associated with the new tooling design and ensure the quality of the parts produced, virtual manufacturing tools must be used to predict the optimal process parameters [START_REF] Lee | Applications of virtual manufacturing in materials processing[END_REF].

(a) (b) (c) (d) Figure 1 -Thermoforming process steps: (a) infrared heating, (b)-(c) forming and consolidation, (d) demoulding

The thermoforming process consists in heating a laminate above the melting temperature of the thermoplastic matrix in an oven; it is then transferred to a press where it is formed and cooled down before demoulding. The sequential steps are presented in Figure 1. Originally developed for thermoplastic sheets, it is now used for manufacturing thermoplastic composites parts as high production volumes are in demand.

The field of forming simulation for woven composites has seen a lot of activity over the last years. While most of the effort is focused on developing constitutive models to represent the behaviour of the dry fabric during draping [START_REF] Boisse | Composites Forming[END_REF], it has been evidenced in recent studies that the forming cannot be considered adiabatic: when the laminate contacts the tool, its temperature rapidly decreases which can lead to local rigidification and formation of wrinkles. The thermomechanical models proposed thus far are either isothermal [START_REF] Machado | A rate-dependent non-orthogonal constitutive model for describing shear behaviour of woven reinforced thermoplastic composites[END_REF] or assume an homogeneous temperature field in the blank at the beginning of the forming step [START_REF] Guzman-Maldonado | Simulation of thermoplastic prepreg thermoforming based on a visco-hyperelastic model and a thermal homogenization[END_REF], which is complex to achieve in practice given the low conductivity of the polymers.

In this paper, we present a more comprehensive approach for the simulation of the thermoforming of thermoplastic preimpregnated composites. Two numerical tools are being developed to simulate the preheating step and the forming step respectively. The temperature field resulting from the infrared heating simulation will be used as an input for the forming simulation.

HEATING SIMULATION

The first step of the thermoforming process aims at bringing the composite blank to a temperature above the melting point of the matrix suitable for forming (295°C in the case of the PA66 matrix) in an oven. Infrared oven are often favoured in industrial composite applications for their fast heat-up time and the good absorptivity of polymers in the infrared wavelength range. However, the low thermal conductivity of the thermoplastic matrix (between 0.2 and 0.5W.m -1 .K -1 [START_REF] Yang | Thermal conductivity[END_REF]) coupled with the spatial configuration of the oven in terms of number of heating elements, distance to the blank etc., can lead to temperature gradients both in plane and through the thickness [START_REF] Cunningham | Modelling of pre-heating of flat panels prior to press forming[END_REF].

The aim of the heating simulation developed is then to predict a realistic three dimensional temperature distribution in the composite stack during the preheating stage, to be later used as an input of a forming simulation; this tool can also be used to optimize an oven configuration for a given processing temperature.

The simulation is developed with Comsol Multiphysics finite elements software, and solves the 3D heat equation in the composite stack:

𝜌𝐶 𝑝 𝜕𝑇 𝜕𝑡 = ∇. (𝜅∇𝑇) + 𝑞 𝑅 + 𝑞 𝑐 (1)
Where q r and q v are the radiative source term and convective losses respectively. The lamps used in this study emit in the near infrared (with 95% of the energy radiated in the range [0.6µm -6µm] at their nominal temperature of 2400K). Transmission experiments showed that both the prepreg material and consolidated plates can be considered opaque in this wavelength range; the radiative term is thus computed using a view factor method via hemicube projection, rather than a more computationally expensive ray tracing approach. The entire environment of the oven (lamps, reflectors and preform) is explicitly modelled. The surfaces that make up the enclosure are assumed to be opaque, diffuse and grey. The convective losses are considered via the determination of a constant convective transfer coefficient.

The modelisation of the lamps is based on the work of Monteix [START_REF] Monteix | Modélisation du chauffage convecto-radiatif de préformes en PET pour la réalisation de corps creux[END_REF]: only the tungsten filament is modelled as a cylinder of equivalent diameter taking into account the emissive surface of the spires. The temperature dependent emissivity and input power-temperature relation are considered for the tungsten.

Solving this equation requires the determination of the thermal properties of the lamps and composite, both in its non-consolidated and consolidated state . The considered material is a glass fabric/polyamide 6.6 thermoplastic composite developed by Solvay. It is supplied as a twill 2/2 woven fabric, powder-impregnated on both sides. Consolidated plates were manufactured using a in-house thermocompression machine.

This study was conducted at Institut Clément Ader, using a novel method based on infrared thermography and inverse analysis to identify the temperatures of the lamps and thermal conductivity tensor of the material [START_REF] Baumard | Determination of thermal contact conductance in thermoplastic prepreg stacks using infrared thermography[END_REF]. This experimental setup also allowed to characterize the inter-ply thermal contact resistances in a non-consolidated composite stack.

(a) (b) Figure 2 -(a) Temperature field on the backface of a [0/90]2s plate heated by a single lamp , (b) Comparison between experimental and simulated temperature

The results of a validation test case with a single lamp heating are presented in Figure 2. The temperature evolution is accurately predicted, with an underestimation of the temperature of the points lying in the convective plume, where the use of a constant heat transfer coefficient is not sufficient. Then, the ability of the model to optimize an infrared oven to minimize the temperature variation on a composite stack is investigated. An oven with 10 lamps disposed in two banks shifted by 4cm, with a center-to-center distance of 8cm between each lamp in the same bank. The oven configuration, viewed from the top, is shown in Figure 3. On each bank (left and right), the lamps are numbered from the closest to the furthest of the plate symmetry axis.

Figure 3 -Oven configuration

The composite stack is heated from the top only. In a first step, only the input power of the lamps is considered as a variable. The cost function to be minimized is therefore:

min 20≤𝑃 𝑖 ≤100 𝑓(𝑷) , 𝑓(𝑷) = ∫(𝑇(𝒙) -𝑇 𝑜𝑏𝑗 )²𝑑𝑉 𝑉 (2)
With P a vector containing the individual input power values of the lamps. Given the symmetry of the oven, the dimension of P is reduced to 5 (i.e., P(L i )= P(R i ) ). The objective temperature is in this case 295°C, which is approximately 30° above the melting temperature of the thermoplastic matrix. Optimization is conducted in steadystate, using a Levenberg-Marquardt algorithm. The lamp input power is constrained between 20% and 100% of the nominal power to avoid under voltage. The temperature on the top of the composite stack along the red dotted line is shown in Figure 4 before and after the optimization procedure. The dashed blue line shows the temperature distribution along the profile with all lamps at the same input power (50% of nominal power). The solid red line corresponds to the temperature field after optimization, and exhibits a much flatter profile: the central hot spot is avoided.

The transfer phase between the infrared oven and the mold and die in the actual process allows considering a loose coupling between the two simulation tools: the three dimensional temperature field that results from the heating simulation will be used as an input for the forming model.

FORMING SIMULATION

The forming simulation is based on a non-orthogonal constitutive model based on the work of Khan [START_REF] Khan | Numerical and experimental analyses of woven composite reinforcement forming using a hypoelastic behaviour[END_REF] and Pierce [START_REF] Pierce | Implementation of a non-orthogonal constitutive model for the finite element simulation of textile composite draping[END_REF]. This model considers the composite as a continuous, anisotropic media, and tracks the orientation of the yarns, thus allowing the prediction of possible defect locations by visualizing the shear angle distribution in the part. This is done via a material user subroutine VUMAT implemented in Abaqus/Explicit. The subroutine works by calculating the fibre directions from the deformation gradient; the incremental strain is rotated from the Green-Naghdi frame used by Abaqus to the fibre frame; the incremental stress is computed in this frame from the constitutive law and converted back to the Green-Naghdi frame. The constitutive tensor assumes the following form in the fiber frame:

𝐶 𝑓 1 = [ 𝐸 1 0 0 0 0 0 0 0 𝐺 12 ] , 𝐶 𝑓 2 = [ 0 0 0 0 𝐸 2 0 0 0 𝐺 12 ] (3)
The superscripts 1, 2 represents the fiber directions, E the Young's modulus and G 12 the shear modulus of the material. Tension-shear coupling and bending stiffness of the yarn are neglected. The advantages of this approach are the relative ease of implementation in commercial finite element codes, and the simplicity of the constitutive material tensor as the computation is conducted in the material frame.

A phenomenological, 3 rd order polynomial expression of the shear stress as a function of shear angle is derived from bias extension tests conducted on the dry fabric at ambient temperature, the resulting expression for the shear modulus being:

𝐺 12 = 𝑎 + 2𝑏. 𝛾 12 + 3𝑐. 𝛾 12 2 (4)
A bias extension test is simulated to validate this formulation. The resulting shear stress as a function of shear angle is plotted in Figure 5 against the experimental results from optical shear angle measurements. The model, originally developed for draping simulations of dry fabric, is currently being extended to take into account the temperature dependent behaviour of the resin through its contribution to the shear modulus of the material. A bias extension campaign at high temperature and various strain rates is planned and experimental hemispherical punch apparatus is developed at Queen's University Belfast to validate the model formulation with stamping experiments. The fabric is held into place by the mean of a spring-actuated clamp to easily vary the clamping force, while the punch itself is mounted on a 60kN Promess press unit.

An example of the shape of the fabric after deformation by the punch is shown in Figure 6, with the camera direction aligned with the warp yarns. The first validation simulations (see Figure 7) show a good qualitative agreement with the experiments in terms of shear angle distribution and draw-in profiles, with the zones of maximum shear appearing at a 45° angle from the principal fiber directions. The hemispherical punch will be used coupled with an infrared oven as a complete, instrumented lab-scale thermoforming setup to investigate the effect of temperature on the type and onset of defects generated.

CONCLUSION

Understanding the thermomechanical behaviour of thermoplastic composites during the thermoforming process is critical to accurately predict the defects that can appear during the manufacturing. In this article, the strategy for the development of a predictive, virtual manufacturing tool that relies on two coupled simulations of the infrared heating and forming steps has been presented. The material characterization procedure has been introduced for the thermophysical and mechanical properties, and is currently underway for high temperatures. The first validation test cases used for the simulations show the pertinence of the approach with promising qualitative results; these will be confirmed with the complete thermoforming setup being finalized.

Figure 4 -

 4 Figure 4 -Temperature along horizontal profile, dashed line: before optimization, solid line: optimized input power

Figure 5 -

 5 Figure 5 -Comparison between experimental and simulated bias extension tests on twill 2/2 dry fabric

Figure 6 -Figure 7 -

 67 Figure 6-Hemispherical forming experiment on dry fabric

ACKNOWLEDGEMENTS

The work was conducted at Institut Clément Ader, Ecole des Mines d'Albi and at the School of Mechanical and Aerospace Engineering, Queen's University Belfast. The authors would like to thank Solvay for supporting this study and providing the test materials.


Jiuchun Gao
email: jiuchun.gao@irccyn.ec-nantes.fr

Anatol Pashkevich
email: anatol.pashkevich@mines-nantes.fr

Stéphane Caro
email: stephane.caro@irccyn.ec-nantes.fr

⋅ A Pashkevich

Manipulator Motion Planning in Redundant Robotic System for Fiber Placement Process

Keywords:

This paper proposes a new methodology for motion planning in redundant robotic system for the high-speed fiber placement technology. The considered system is composed of a 6-axis industrial robot and a one-axis actuated positioner. Compared to the previous works assuming constant rate, the particularity of this work lies in the time-optimal motion coordination based on the combinatorial optimization and the polynomial approximation techniques. Advantages of the developed technique are confirmed by an application example that deals with a planar robotic system for fiber placement.

Introduction

Composites have been increasingly used in aerospace and automotive industries because of their good strength-to-weight ratio and durability [START_REF] Kaw | Mechanics of Composite Materials[END_REF]. For fabricating complex composite parts, fiber placement is an important technique [START_REF] Gay | Composite Materials: Design and Applications[END_REF]. The fiber placement process can be implemented by using either specifically designed machines or robotic systems, which are redundant in this application. The machines have no limitations on the component size, but they usually require large work floor areas and are expensive [START_REF] Gallet-Hamlyn | Multiple-use robots for composite part manufacturing[END_REF]. Compared to the process-dedicated machines, the robotic system (usually includes a 6-axis robot and a one-axis positioner) is relatively cheap and flexible allowing changing the product type easily.

In robotic fiber placement, planning of manipulator motions is an important issue. The main difficulty here arises due to the redundancy of the robotic system with respect to the manufacturing task. In literature, there are a number of works that deal with the redundancy resolution. Relevant techniques are usually based on the pseudo inverse of the kinematic Jacobian. However, they can be hardly applied to the considered problem because they do not allow generating optimal trajectories satisfying real-life industrial requirements [START_REF] Kazerounian | Redundancy resolution of serial manipulators based on robot dynamics[END_REF]. There are a few researches dealing with this problem in particular technological process. For laser-cutting and arc-welding applications, relevant numerical techniques were proposed in [START_REF] Dolgui | Manipulator motion planning for high-speed robotic laser cutting[END_REF], but the main assumption (constant speed of the tool) is not in agreement with the nature of the fiber placement. For fiber placement process, the motion planning problem was studied by Martinec et al. [START_REF] Martinec | Calculation of the robot trajectory for the optimum directional orientation of fibre placement in the manufacture of composite profile frames[END_REF] and Mlynek and Martinec [START_REF] Mlynek | Mathematical model of composite manufacture and calculation of robot trajectory[END_REF], who also assumed the end-effector velocity is constant. For this reason, the problem of optimal manipulator motion planning for the fiber placement has not been solved completely yet.

This paper proposes a new methodology that aims at obtaining the time-optimal trajectory taking into account real industrial constraints. The proposed approach is based on the discretization of the redundant variable and dynamic programming technique. Further, a smoothing technique is also applied.

Robotic Fiber Placement Problem

A typical robotic fiber placement system (a 6-axis robot and a one-axis positioner) is presented in Fig. 1. The workpiece is mounted on the positioner that is able to change its orientation in order to improve accessibility of certain desired zones by the robot end-effector. The latter is attached to the manipulator flange and ensures placement of the fiber tows in the desired locations. To describe the fiber placement task, let us introduce several coordinate systems attached to the robotic manipulator, the positioner, the workpiece and the robot end-effector (see Fig. 1). It is assumed that the fiber placement task is defined with respect to the frame F W attached to the workpiece. The task is presented as a set of n discrete task locations (position and orientation) sampled on the given placement path. For each of these locations, it is also defined a frame F TLi , i = 1, 2, . . . n, in such a way that the axis X TLi is directed along the path direction and the axis Z TLi is normal to the workpiece surface pointing outside of it. Using these notations, the fiber placement task can be presented as a sequence of homogenous 4 × 4 matrices describing the mapping from the frame

F W to F TLi W T TLi = n i s i a i p i 0 0 0 1 4 × 4 = Rðφ i Þ p i 0 0 0 1 ; i = 1, 2, . . . , n ð1Þ
where p i = ðx i , y i , z i Þ T is the position vector of the frame

F TL i origin, φ i = ðα i , β i , γ i Þ T
is the orientation vector in accordance with the Z-Y-X Euler angles definition, and n i , s i and a i are the unit vectors giving directions of the axes X, Y, and Z. To execute the given fiber placement task, the robot end-effector must visit the frames (1) sequentially and as fast as possible.

Robotic System Model

The robotic system started in this paper is composed of two components: an industrial robot and a positioner. The positioner model describes the mapping from the positioner base to the positioner flange that depends on the positioner joint coordinate q p . This transformation can also be presented as a 4 × 4 homogenous transformation matrix posðq p Þ = R X ðq p Þ. Further, the matrix function posðq p Þ will be referred to as the positioner direct kinematics that depends on the positioner parameters a and d (see Fig. 1). The robot model describes the mapping from the robot base to the robot flange. For serial manipulators, it can be obtained using the D-H technique and presented as the product of the 4 × 4 homogenous matrices

rob ðq r Þ = RB T 1 ðq 1 Þ ⋅ 1 T 2 ðq 2 Þ ⋅ 2 T 3 ðq 3 Þ ⋅ 3 T 4 ðq 4 Þ ⋅ 4 T 5 ðq 5 Þ ⋅ 5 T 6 ðq 6 Þ ð2Þ
where the matrix function robðq r Þ is also called the robot direct kinematics; and q r is the vector of the robot joint coordinates.

To derive the designed model of the robotic system, let us present the task frames in two ways, using the robot and positioner kinematics as follows

0 T TLi ðq p Þ = 0 T PB ⋅ posðq p Þ ⋅ PF T W ⋅ W T TLi ð4Þ
where T is the matrix which defines the mapping from the frame referred to the left superscript to the frame referred to the right subscript.

After equating these two presentations, the robotic model can be expressed as

0 T RB ⋅ robðq r Þ ⋅ RF T T ⋅ T T TLi = 0 T PB ⋅ posðq p Þ ⋅ PF T W ⋅ W T TLi ð5Þ
It should be stressed that this model describes the geometry of a closed kinematic loop whose configuration depends on the joint coordinates q r and q p of the robotic manipulator and the positioner, respectively. Besides, it should be mentioned that this model is redundant with respect to the joint variables since it gives six independent scalar equations for seven unknowns. The latter does not allow the users to obtain the unique solution for the robot and positioner motions corresponding to the given fiber placement path. On the other hand, it provides some rooms for the optimization of robot /positioner motions when tracking a specified curve.

Optimization Problem

To utilize the redundancy in the best way, it is reasonable to partition the required movement between the robot and the positioner ensuring that the end-effector executes the desired task with smooth motion in the shortest time.

Motion Generation Problem

To present the above mentioned problem in a more formal way, let us introduce the functions q r ðtÞ and q p ðtÞ that describe the robot and positioner motions on the time interval t ∈ ½0, T. In addition, let us define the time instances t 1 , t 2 , . . . t n f g corresponding to the cases where the robot end-effector visits the task nodes

F TLi i = 1, 2, . . . n j f
g defined by Eq. ( 1), where t 1 = 0, t n = T. Using this notation, the problem can be presented as minimization of the total travelling time T → min q r ðtÞ, q p ðtÞ ð6Þ over q r ðtÞ and q p ðtÞ subjected to the equality/inequality constraints, which are imposed by the prescribed task

0 T RB ⋅ robðq r ðt i ÞÞ ⋅ RF T T ⋅ T T TLi = 0 T PB ⋅ posðq p ðt i ÞÞ ⋅ PF T W ⋅ W T TLi ; i = 1, 2, . . . n ð7Þ
and which describe the capacities of the robot/positioner actuators

q min j ≤ q j ðt i Þ ≤ q max j ð8aÞ q̇m in j ≤ qjðt i Þ ≤ q̇m ax j ð8bÞ q̈m in j ≤ q̈jðt i Þ ≤ q̈m ax j ð8cÞ
where j = 0, 1, . . . 6 is the common index for the joint variables, and j = 0 corresponds to q p . Besides, the collision constraints are defined as follows: cols ðq p ðtÞ, q r ðtÞÞ = 0; ∀t ∈ ½0, T ð 9Þ must be taken into account, where the binary function cols ð.Þ verifies the intersections between the system components (robot and positioner links, fixture, etc.). For this optimization problem, which aims at finding continuous functions q r ðtÞ and q p ðtÞ, there is no standard technique that can be straightforwardly applied to. The main difficulty here is related to the equality constraints that are written for the unknown time instances t 1 , t 2 , . . . t n . Besides, this problem is highly nonlinear and includes redundant variables (for the non-redundant case without the collision constraints the problem was solved by Bobrow et al. [START_REF] Bobrow | Time-optimal control of robotic manipulators along specified paths[END_REF]). For these reasons, this paper proposes a discrete optimization based approach that is able to take into account both the system redundancy and the collision constraints.

Search Space Discretization

For the considered robotic system, there is a single redundant variable with respect to the task. It is convenient here to consider q p as the redundant one. This allows us, after solving the positioner direct kinematics for any given q p , straightforwardly apply the manipulator inverse kinematics and find corresponding configurations q r of the robot. This approach permits us to take into account explicitly the equality constraint in Eq. [START_REF] Mlynek | Mathematical model of composite manufacture and calculation of robot trajectory[END_REF].

To present the problem in a discrete way, let us sample the allowable domain of the redundant variable q p ∈ ½q min p , q max p with the step Δq p q ðkÞ p = q min p + Δq p ⋅ k; k = 0, 1, . . . m ð10Þ

where m = ðq max p q min p Þ ̸Δq p . Then, applying sequentially the positioner direct kinematics and the manipulator inverse kinematics in accordance with Eq. ( 7), one can get a set of possible configuration states for the robotic system. For the unique mapping from the task space to joint space, let us define the configuration index vector μ = ðμ 1 , μ 2 , μ 3 Þ T that corresponds to the manipulator posture, where μ 1 , μ 2 and μ 3 specify the shoulder, elbow and wrist indices respectively. Then, the manipulator configuration corresponding to the given q ðkÞ p can be obtained as follows:

q ðkÞ r ðt i Þ = g -1 r ðg p ðq ðkÞ p ðt i ÞÞ, μÞ; k = 0, 1, . . . m; i = 1, 2, . . . n ð11Þ
where t i specifies the unknown time instant corresponding to the task point TL i , the functions g p ð.Þ and g -1 r ð.Þ denote the positioner direct kinematics and the manipulator inverse kinematics, respectively. Therefore, for each task location we can generate a number of configuration states, i.e. TL i → L ðk, iÞ c ; ∀k, i, where L ðk, iÞ c = ðq ðkÞ r ðt i Þ, q ðkÞ p ðt i ÞÞ will be further referred to as the location cell. While considering the task locations are ordered in time, the original sequence of TL i described by Eq. ( 1) may be converted into a directed graph presented in Fig. 2. It should be noted that some of the configurations generated by Eq. ( 11) should be excluded from the graph because of violation of the collision constraints (9) or the actuator joint limits (8a). Besides, from an engineering point of view, it is prudent to avoid some configurations that are close to the manipulator singular postures. These cases correspond to the "inadmissible" in Fig. 2, which are not connected to any neighbor. It is clear that due to time-irreversibility, the allowable connections between the graph nodes are limited to the subsequent configuration states

L ðk i , iÞ c → L ðk i + 1 , i + 1Þ c
, and the edge weights correspond to the minimum travelling time that are restricted by the actuator capacities expressed by Eq. (8a-8c).

Using the discrete search space, the considered problem can be transformed to the searching of the shortest path on the above presented graph. In accordance with the physical sense, the initial and final nodes must belong to the sets L ðk 1 , 1Þ c , ∀k 1 È É and L ðk n , nÞ c , ∀k n È É respectively. In the frame of this notation, the desired solution can be represented as the sequence

L ðk 1 , 1Þ c È É → L ðk 2 , 2Þ c È É → . . . L ðk n , nÞ c È É .
In accordance with the actuator constraints (8b), the distance between subsequent nodes can be evaluated as the displacement time for the slowest joint

dist ðL ðk i , iÞ c , L ðk i + 1 , i + 1Þ c Þ = max j = 0, ..6 q ðk i Þ j, i -q ðk i + 1 Þ j, i + 1 ̸ q̇m ax j ð12Þ
The latter allows us to present the objective function (travelling time) as follows

T = ∑ n -1 i = 1 dist ðL ðk i , iÞ c , L ðk i + 1 , i + 1Þ c Þ ð 13Þ
that depends on the indices k 1 , k 2 , . . . k n . It should be noted that the applied method of edge weights computing automatically takes into account the velocity constraints (8b), but the acceleration constraints (8c) must be examined as follows

2 Δt i ðq ðk i + 1 Þ j, i + 1 -q ðk i Þ j, i Þ -Δt i + 1 ðq ðk i Þ j, i -q ðk i -1 Þ j, i -1 Þ ̸ Δt i + 1 Δt i ðΔt i + 1 + Δt i Þ ≤ q̈m ax j ð14Þ
which is based on the second order approximation of the corresponding functions q r ðtÞ and q p ðtÞ. Here,

Δt i + 1 = dist ðL ðk i , iÞ c , L ðk i + 1 , i + 1Þ c Þ and Δt i = dist ðL ðk i , iÞ c , L ðk i -1 , i -1Þ c Þ.

Path Planning Algorithm

After discretization, the original optimization problem ( 6)-( 9) is converted to the combinatorial one, which can generally be transformed to the classical shortest path search on the graph (by simply adding virtual nodes at the beginning and the end). However, this straightforward approach is extremely time-consuming. For instance, it takes more than 20 h to find desired solution for 100 task locations and the positioner discretization step 1 • (Intel ® i5 2.67 GHz). Besides, standard techniques are not able to respect acceleration constraints (8c) that are very important here. For these reasons, a more efficient optimization algorithm was proposed that takes into account particularities of the graph describing the search space, which can be also treated as a multi-layer network (see Fig. 2).

The developed algorithm is based on the dynamic programming principle, which breaks down the full-size problem into a set of sub-problems [9], aiming at finding the shortest path from the initial node set L ðk

d k, i + 1 = min k ′ d k ′ , i + dist ðL ðk, i + 1Þ c , L ðk ′ , iÞ c Þ n o ; ð15Þ
that is applied sequentially starting from the second layer, i.e. i = 1, 2, . . . n -1. Finally, after selection of the minimum length d k, i + 1 corresponding to the final layer and applying the backtracking, one can get the desired optimal path. Therefore, the desired solution for n sets

L ðk 1 , 1Þ c È É → L ðk 2 , 2Þ c È É → . . . L ðk n , nÞ c È É is obtained.
The desired path is described by the recorded indices k 1 , k 2 , . . . k n f g . The developed algorithm is rather time-efficient and perfectly suits industrial requirements. In particular, it takes about 30 s to find the optimal solution for the above mentioned problem (100 task locations, discretization step 1 • ).

Application Example

To show the efficiency of the proposed technique, let us apply it to a fiber placement system (a two-axis robot and a one-axis positioner). This system possesses a 1-dof redundancy with respect to the task. The geometrical parameters of the system and the actuator constraints are given in Fig. 3. The number of task locations is equal to n = 100 and the discretization step is equal to Δq p = 0.5 • . Fig. 3 Fiber placement system: two-axis robot and one-axis positioner A desired solution was obtained by implementing the developed algorithm described in Sect. 4. The result is presented in Fig. 4, which contains the displacement and velocity profiles for all three actuators. It should be noted that the trajectories obtained by direct application of the dynamic programming technique include oscillations in derivatives, which are undesirable from engineering point of view. For this reason, a specially developed smoothing technique was applied to improve locally the velocity and acceleration profiles. This technique is based on the polynomial approximation of the redundant variable profile and relevant correction of the displacement intervals, to ensure that the velocities and accelerations are within their limits. It is clear that this modification slightly increases the travelling time, but it yields important benefits to the actuator working conditions.

Conclusions

This paper contributes to a new methodology to optimize the manipulator motion in redundant robotic system for the fiber placement process. In contrast to the previous works, the proposed approach generates the time-optimal smooth trajectories and allows using the full capacities of the robotic system expressed as the velocity, acceleration and collision constraints. The developed algorithm is based on dynamic programming technique, and leads to significant computing time reduction compared to the standard shortest-path algorithms. Advantages of the developed technique are highlighted with an application example that deals with a planar robotic system for fiber placement. Future research works will concentrate on its integration into industrial robot control system software.

Fig. 1

 1 Fig. 1 Robotic fiber placement system

Fig. 2

 2 Fig. 2 Graph-based presentation of the discrete search space

  basic idea, let us denote d k, i as the length of the shortest path connecting one of the initial nodes L ðk 1 , 1Þ c , ∀k 1 È É to the current node L ðk ′ , iÞ c n o . Then, taking into account the additivity of the objective (13), the shortest path for the nodes belonging to the next layer L ðk, i + 1Þ c , ∀k È É can be found by combining the optimal solutions for the previous layer L ðk ′ , iÞ c , ∀k ′ n o and the distances between the nodes with the indices i and i + 1. The latter corresponds to the formula

T TLi ðq r Þ = 0 T RB ⋅ robðq r Þ ⋅ RF T T ⋅ T T TLið3Þ

Acknowledgments This work is supported by the China Scholarship Council (Grant N°2 01404490018). The authors also acknowledge CETIM for the motivation of this research work.

Fig. 4 Optimal trajectories for the robotic fiber placement system


Sylvain Deville

Eduardo Saiz

Ravi K Nalla

Antoni P Tomsia
email: aptomsia@lbl.gov

Freezing as a Path to Build Complex Composites

published or not. The documents may come

L'archive ouverte

Text:

While the potential of layered materials has long been recognized (1), their creation requires solving a two-fold problem, namely the design of optimum microstructures, and development of fabrication procedures to implement these designs.

Natural materials such as nacre offer a wealth of information to guide such a design process (2,3). The unique properties of natural layered materials are achieved through a fine control of the layer thickness, selection of the right components, and manipulation of roughness and adhesion at the organic-inorganic interface (4,5). The ideal fabrication process has to be not only simple, but also adaptable enough to fabricate layers with a large number of material combinations and a wide range of layer dimensions. Previous techniques for mimicking nacre are bottom-up chemical approaches (6,7) that are intrinsically limited to a narrow range of materials exhibiting the proper interfacial reactions and compatibility. Other techniques offer only a coarse control of the layer thickness or have practical limitations regarding the number of layers that can be fabricated (7,8).

In sea ice, pure hexagonal ice platelets with randomly oriented horizontal c crystallographic axes are formed, and the various impurities originally present in sea water (salt, biological organisms, etc.) are expelled from the forming ice and entrapped within channels between the ice crystals (9). We apply this natural principle to ceramic particles dispersed in water to build sophisticated, nacre-like architectures in a simple two-step approach. Ice-templated (IT), porous, layered materials-with layers as small as 1 m-are first fabricated through a freeze-casting process, which involves the controlled unidirectional freezing of ceramic suspensions. These porous scaffolds are then filled with a selected second phase (organic or inorganic) to fabricate dense composites. Using a natural, self-organizing phenomenon, we allow nature to guide the design and processing.

The physics of water freezing has drawn the attention of scientists for a long time.

With few exceptions (10), much of this work has concentrated on the freezing of pure water or very dilute suspensions (9,11). This phenomenon is critical for various applications, such as cryo-preservation of biological cell suspensions (12) and the purification of pollutants (13). An important observation in these studies is that, during the freezing of such suspensions, there is a critical particle size (11) above which the suspended particles will be trapped by the moving water-ice front. Another important observation is that the hexagonal ice crystals exhibit strong anisotropic growth kinetics, varying over about two orders of magnitude with crystallographic orientation. Under steady-state conditions, it is possible to grow ice crystals in the form of platelets, with a very high aspect ratio. The ice thus formed will have a lamellar microstructure, with the lamellae thickness depending mainly on the speed of the freezing front. We designed a simple experimental setup (Fig. S1) that allowed us to precisely control the freezing kinetics. By freezing concentrated suspensions containing ceramic particles with suitable granulometry, we are able to build homogeneous, layered, porous scaffolds.

In the method proposed here (supporting online text), directional freezing of the ceramic slurries is achieved by pouring them into Teflon molds placed between two copper cold fingers (Fig. S1) whose temperature is regulated to control the speed of the solidification-front (Fig. S2). As in nature, during the freezing of sea water, the ceramic particles concentrate in the space between the ice crystals (Fig. 1A). When the freezing rate increases, the magnitude of supercooling ahead of the solidifying interface is increased and, as a result, the tip radius of the ice crystals decreases. A finer microstructure is thus obtained without affecting the long-range order of the entire structure. Afterwards, the ice is sublimated by freeze drying, such that a ceramic scaffold whose microstructure is a negative replica of the ice is produced (Fig. 1F). We controlled the growth of lamellar ice by adjusting the freezing kinetics. In this way, we achieved a layered microstructure with relevant dimensions that vary over two decades PAGE 4 (Fig. 2A), from 1 μm (almost the same as nacre, typically ~0.5 m) (14) to 200 μm, without affecting the ordered architecture. To a large extent, the mesostructure of natural materials determines their mechanical response (15), and it has been difficult to replicate this synthetically. Our results indicate that, by controlling the freezing kinetics and patterns of the cold finger, it is also possible to build mesostructural features and gradients (Fig. 2F) that could optimize the mechanical response of the final materialsfor example, by stiffening the structure and limiting torsion (15).

The IT porous scaffolds obtained by this process exhibit striking similarities to the meso-and micro-structure of the inorganic component of nacre (Fig. 1). The inorganic layers are parallel to each other and very homogeneous throughout the entire sample (Fig. 1E). Particles trapped in between the ice dendrites (Fig. 1A) lead to a dendritic surface roughness of the walls (Fig. 1F), just as in nacre (Fig. 1D) (16). Finally, some dendrites span the channels between the lamellae (Fig. 2B), mimicking the tiny inorganic bridges linking the inorganic platelets of nacre, which are believed to increase the fracture resistance (17).

The inorganic portion represents 95% of the volume of nacre, but its highly specific properties (in particular its great toughness) are due to the interaction of this inorganic component with the organic phase (protein) that is found between the calcium carbonate platelets (14). To obtain similar synthetic materials, the next step is to fill the porous ceramic scaffolds with a second phase. For example, we filled the IT scaffolds with a simple organic phase (epoxy) or with an inorganic component (metal, Fig. 1E).

Nature shows that the optimum fracture properties are encountered not only when the organic/inorganic interface is strong, but also when delamination at the organic/inorganic interface occurs before the crack goes across the stiff, brittle layer. In the IT composites, extensive crack deflection at the interface between layers was observed (Figs 3B and3D). As in nacre (Fig 3C ), this delamination creates tortuous PAGE 4 cracks that propagate in a stable manner (Fig 3A ) and increases the toughness of the materials. It is believed that nature manipulates adhesion in two ways-mechanical and chemical. In nacre, this is done by controlling the roughness and the highly specific properties of the polymer adhesive phase (4). Our process allows us to control the morphology of the inorganic layers and the chemistry of the interface. For example, the mechanical response of alumina/Al-Si (45/55 vol.%) layered composites (Fig. 1E) can be manipulated by controlling the interfacial bonding. By adding as little as 0.5 wt.% Ti (known to segregate at the metal/ceramic interfaces ( 18)) to the aluminum eutectic, the strength increases from 400 to 600 MPa and fracture toughness from 5.5 to 10 MPam (Fig. 3E).

Our technique shows promise for a large number of applications that require tailored composite materials. One such scientific challenge that could be solved is the development of new biomaterials for orthopaedic applications (19). Despite extensive efforts in the development of porous scaffolds for bone regeneration, all porous materials have an inherent lack of strength associated with porosity. By applying freezing to commercial hydroxyapatite (HAP, the mineral component of bone) powder suspensions, we processed IT highly porous lamellar scaffolds that are four times stronger in compression than conventional porous HAP (Fig. 4). These IT scaffolds exhibit well-defined pore connectivity along with directional and completely open porosity of an adequate size to allow bone ingrowth (20). Hence, most of the current shortcomings (low strength, random organization, multiple pore size, uncontrolled pore connectivity) that plague such bone substitutes are eliminated by this innovative approach.

Current ceramic and metallic implant materials have serious shortcomings due to the mismatch of physical properties with those of bone. In bone, intrinsically weak materials, such as calcium phosphates and collagen, are combined into composites exhibiting intermediate modulus (10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20), fairly high strength (30 to 200 MPa), and high work of fracture (100-1000 J/m 2 ) [START_REF] An | Mechanical testing of bone and the bone-implant interface[END_REF]. The unique properties of bone arise from the controlled integration of the organic (collagen) and inorganic (apatite) components ( 5) with a sophisticated architecture from the nano-to meso-levels. Our approach to the problem is to infiltrate the IT porous HAP scaffolds with a second organic phase with tailored biodegradability. Because the biodegradation rates of the scaffold and the infiltrated compound can be designed to be different, porosity will be created in situ to allow bone ingrowth. Using this approach, we have been able to fabricate HAP-based composites, with stiffness (10 GPa), strength (150 MPa), and work of fracture (220 J/m 2 ) that match that of compact bone for an equivalent mineral/organic content (around 60/40 vol.%).

Figure 1. Processing principles and materials. While the ceramic slurry is freezing, the growing ice crystals expel the ceramic particles, creating a lamellar microstructure oriented in a direction parallel to the movement of the freezing front (A). For highly concentrated slurries, the interaction between particles becomes critical: a small of particles are entrapped within the ice crystals by tip-splitting and subsequent healing (A), leading to the formation of inorganic bridges between adjacent walls. Dense composites are obtained by infiltrating the porous lamellar ceramic with a second phase (e.g., a polymer or a liquid metal). Natural nacre has a brick-mortar-bridges microstructure where inorganic calcium carbonate layers are held together by organic protein "glue" (B and C), the roughness of the inorganic walls (D) is a key contributor to the final mechanical properties of nacre. The layered microstructure of the IT dense composites resembles that of nacre (for example the alumina/Al-Si composite in E). The particles entrapped between the ice dendrites generate a characteristic roughness on the lamella surface (F) that mimics that of the inorganic component of nacre. Scale bars (B) 5 m (C) 0.5 m (D) 0.3 m (E) 300 m (F) 10 m. [START_REF] Almirall | [END_REF](23)(24)[START_REF] Bignon | Optimization of the porous structure of calcium phosphate implants for bone substitutes and in situ release of active principles[END_REF][START_REF] Ramay | [END_REF](27)(28)(29)(30) vs. IT porous HAP scaffolds. The typical pore sizes of conventional porous HAP scaffolds are of the order of 100 m to 800m, in order to allow bone ingrowth. In the IT materials exhibiting the greatest strength, the pores are typically ~20 by ~200 microns wide and several millimetres long-previous studies have indicated that these dimensions are adequate for bone tissue engineering (20). For the IT porous materials, compression is applied in the direction parallel to the ceramic layers. Each style of points corresponds to a different source in the literature. The presence of inorganic bridges between the ceramic layers (a feature that parallels the microstructure of nacre) prevents Euler buckling of the ceramic layers and contributes to the high strength. The insert shows typical compression load-displacement curves for materials with 56% porosity (three different samples shown here). The samples fail gradually, and due to the large degree of control of hierarchical architecture, the mechanical behavior is very consistent from one sample to another.

Figure 2 .

 2 Figure 2. Microstructural control at several levels. (A) Effect of the speed of the

Figure 3 .

 3 Figure 3. Mechanical response of natural and synthetic IT composites. The 3-points

Figure 4 .

 4 Figure 4. Compressive strength of porous hydroxyapatite scaffolds. Results from literature

Acknowledgments

This work was supported by the National Institute of Health (NIH/NIDCR) under grant No. 5R01 DE015633 (Novel Scaffolds for Tissue Engineering and Bone-Like Composites) and by the Director, Office of Science, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering of the Department of Energy under Contract No. DE-AC03-76SF00098 (Metal/Ceramic Composites). The authors wish to thank Prof. Robert O. Ritchie for useful discussions and James Wu for help with the synthesis of the aluminum-infiltrated composites.


Sylvain Deville
email: sdeville@lbl.gov

Eduardo Saiz

Antoni Tomsia

Freeze casting of hydroxyapatite scaffolds for bone tissue engineering

Although extensive efforts have been put into the development of porous scaffolds for bone regeneration, with encouraging results, all porous materials have a common limitation: the inherent lack of strength associated with porosity. Hence, the development of porous hydroxyapatite scaffolds has been hindered to non-load bearing applications. We report here how freeze-casting can be applied to synthesize porous hydroxyapatite scaffolds exhibiting unusually high compressive strength, e.g. up to 145 MPa for 47% porosity and 65 MPa for 56% porosity. The materials are characterized by well-defined pore connectivity along with directional and completely open porosity. Various parameters affecting the porosity and compressive strength have been investigated, including initial slurry concentration, freezing rate, and sintering conditions. The implications and potential application as bone substitute are discussed. These results might open the way for hydroxyapatite-based materials designed for load-bearing applications. The biological response of these materials is yet to be tested.

Introduction

The requirements for a synthetic bone substitute appear deceptively simple, that is, to supply a porous matrix with interconnecting porosity that promotes rapid bone ingrowth, yet possesses sufficient strength to prevent crushing under physiological loads during integration and healing.

The ideal bone substitute is not a material that interacts as little as possible with the surrounding tissues, but one that will form a secure bond with the tissues by allowing, and even encouraging new cells to grow and penetrate. One way to achieve this is to use a material that is osteophilic and porous, so that new tissue, and ultimately new bone, can be induced to grow into the pores and help prevent loosening and movement of the implant. Resorbable bone replacements have been developed from inorganic materials that are very similar to the apatite composition of natural bone [START_REF] Suchanek | Processing and properties of hydroxyapatite-based biomaterials for use as hard tissue replacement implants[END_REF].

In recent years, considerable attention has been given to the development of fabrication methods to prepare porous ceramic scaffolds for osseous tissue regeneration [START_REF] Deville | Freezing as a path to build complex composites[END_REF][START_REF] Almirall | Fabrication of low temperature macroporous hydroxyapatite scaffolds by foaming and hydrolysis of an alpha-TCP paste[END_REF][START_REF] Charriere | Hydroxyapatite cement scaffolds with controlled macroporosity: fabrication protocol and mechanical properties[END_REF][START_REF] Del Real | A new method to produce macropores in calcium phosphate cements[END_REF][START_REF] Kawata | Development of porous ceramics with well-controlled porosities and pore sizes from apatite fibers and their evaluations[END_REF][START_REF] Milosevski | Preparation and properties of dense and porous calcium phosphate[END_REF][START_REF] Ramay | Preparation of porous hydroxyapatite scaffolds by combination of the gel-casting and polymer sponge methods[END_REF][START_REF] Tian | Preparation of porous hydroxyapatite[END_REF]. The ideal fabrication technique should produce complex-shaped scaffolds with controlled pore size, shape and orientation in a reliable and economical way. However, all porous materials have a common limitation: the inherent lack of strength associated with porosity. Hence, their application tends to be limited to low-stress locations, such as broken jaws or fractured skulls. Therefore, the unresolved dilemma is how to design and create a scaffold that is both porous and strong.

Freeze casting is a simple technique to produce porous complex-shaped ceramic or polymeric parts [START_REF] Fukasawa | Pore structure of porous ceramics synthesized from water-based slurry by freeze-dry process[END_REF]. In freeze casting, a ceramic slurry is poured into a mold and then frozen. The frozen solvent acts temporarily as a binder to hold the part together for demolding. Subsequently, the part is subject to freeze-drying to sublimate the solvent under vacuum, avoiding the drying stresses and shrinkage that may lead to cracks and warping during normal drying. After drying, the compacts are sintered in order to fabricate a porous material with improved strength, stiffness and desired porosity. The result is a scaffold with a complex and often anisotropic porous microstructure generated during freezing. By controlling the growth direction of the ice crystals, it is possible to impose a preferential orientation for the porosity in the final material [START_REF] Kang | Fabrication of porous gelatin scaffolds for tissue engineering[END_REF].

The technique was applied more specifically to polymeric materials, for tissue engineering. A wide variety of materials have already been investigated, including chitin [START_REF] Madihally | Porous chitosan scaffolds for tissue engineering[END_REF], gelatin [START_REF] Kang | Fabrication of porous gelatin scaffolds for tissue engineering[END_REF][START_REF] Huang | In vitro characterization of chitosan-gelatin scaffolds for tissue engineering[END_REF], collagen [START_REF] Schoof | Dendritic ice morphology in unidirectionally solidified collagen suspensions[END_REF], PLA [START_REF] Zhang | Poly([alpha]-hydroxyl acids)/hydroxyapatite porous composites for bone-tissue engineering. I. Preparation and morphology[END_REF][START_REF] Ho | Preparation of porous scaffolds by using freeze-extraction and freeze-gelation methods[END_REF], PDLLA [START_REF] Zhang | Poly([alpha]-hydroxyl acids)/hydroxyapatite porous composites for bone-tissue engineering. I. Preparation and morphology[END_REF][START_REF] Maquet | Porous poly([alpha]hydroxyacid)/Bioglass(R) composite scaffolds for bone tissue engineering. I: preparation and in vitro characterisation[END_REF], PLGA [START_REF] Zhang | Poly([alpha]-hydroxyl acids)/hydroxyapatite porous composites for bone-tissue engineering. I. Preparation and morphology[END_REF][START_REF] Maquet | Porous poly([alpha]hydroxyacid)/Bioglass(R) composite scaffolds for bone tissue engineering. I: preparation and in vitro characterisation[END_REF], poly(HEMA) [START_REF] Oxley | Macroporous hydrogels for biomedical applications: methodology and morphology[END_REF], agarose [START_REF] Stokols | The fabrication and characterization of linearly oriented nerve guidance scaffolds for spinal cord injury[END_REF], sericin [START_REF] Tao | Preparation and structure of porous silk sericin materials[END_REF] and alginate [START_REF] Ho | Preparation of porous scaffolds by using freeze-extraction and freeze-gelation methods[END_REF][START_REF] Shapiro | Novel alginate sponges for cell culture and transplantation[END_REF][START_REF] Zmora | Tailoring the pore architecture in 3-D alginate scaffolds by controlling the freezing regime during fabrication[END_REF][START_REF] Chung | Preparation of alginate/galactosylated chitosan scaffold for hepatocyte attachment[END_REF][START_REF] Tampieri | HA/alginate hybrid composites prepared through bio-inspired nucleation[END_REF]. Although not stiff and strong enough for load-bearing applications, all these materials have in common a homogeneous structure with open porosity, favorable for rapid cell proliferation. In particular, pore size and its structure can be controlled by heat transfer rate.

We show here how freeze-casting can be applied to hydroxyapatite, an osteophilic ceramic related to the inorganic component of bone, to process bone substitute materials with suitable physical and mechanical properties. In particular, we describe here how the processing conditions (concentration, freezing rate, sintering) affects the scaffold characteristics (size and amount of porosity, compressive strength) and discuss the limits of the technique.

Experimental techniques

The porous inorganic scaffolds were produced by controlled freezing of hydroxyapatite (HAP) slurries. Slurries were prepared by mixing distilled water with a small amount (typically 1 wt.%) of ammonium polymethacrylate anionic dispersant (Darvan 811, R. T. Vanderbilt Co., Norwalk, CT), an organic binder (1 wt.%, polyvinyl alcohol) and the hydroxyapatite powder in various content (Hydroxyapatite#30, Trans-Tech, Adamstown, MD), depending on the desired total porosity. Slurries were ball-milled for 20 hrs with alumina balls and de-aired by stirring in a vacuum desiccator.

Freezing of the slurries was done by pouring them into a PTFE mold, placed on a Cu cold finger whose temperature is controlled using liquid nitrogen and a ring heater. Typical cooling rates range from 0.5 to 5°C/min. Frozen samples were freeze-dried (Freeze-dryer 8, Labconco, Kansas City, MI) at low temperature (-100 °C, temperature of the cold finger of the freeze-dryer) and low pressure (1300 Pa), for 24 hrs. Sintering of the green bodies was done in an air furnace (1216BL, CM Furnaces Inc., Bloomfield, NJ), in the temperature range of 1250-1350°C, with heating and cooling rates of 2°C/min. Samples of 4x5x5 mm 3 were cut for compression tests. Compression tests were carried out on a testing machine (Instron 1122, Instron, Norwood, MA), with a crosshead speed of 0.2 mm/min. Apparent porosity was derived from the density, measured by Archimede's method. The microstructure of the samples was analyzed by scanning electron microscopy (SEM), and the phases present were identified by X-ray diffraction. SEM work was performed using an environmental scanning electron microscope (S-4300SE/N, Hitachi, Pleasanton, CA). X-ray diffraction was performed using a Siemens D500 diffractometer (Siemens AG, Munich, Germany).

Results

By directional freezing of the slurry, we force the particles in suspension to be rejected from the moving ice front and piled up between the growing columnar ice crystals [START_REF] Deville | Freezing as a path to build complex composites[END_REF]. Afterwards, the ice is sublimated by freeze-drying, such that a ceramic scaffold whose microstructure is a negative replica of the ice is produced. The porosity of the sintered materials is a replica of the ice structure before drying.

Sample general features

According to their microstructure, the samples can clearly be divided into three distinctive zones (Fig. 1), each of them characterized by a different shape and dimensions of the pores. In zone 1, the closest to the initial cold finger, no porosity at all is observed and the material is dense. In the second zone, the material is characterized by a cellular morphology. Finally, in the upper zone (zone 3), the ceramic is lamellar, with long parallel pores aligned in the movement direction of the ice front. The underlying reason for this architecture will be described later in the discussion.

Unless explicitly stated otherwise, this paper will focus on the last zone (zone 3).

The microstructure of the lamellar zone is to some extent very similar to that obtained on polymeric materials for tissue engineering processed with the same technique. It consists of hydroxyapatite plates with flat interconnected macropores between them, aligned along the ice growth direction (Fig. 2). On the internal walls of the lamellae, a dendritic, branch-like structure can be observed following the microscopic ice formation. Some ceramic bridges, linking adjacent plates, are also observed. The sintered scaffolds had relative densities between 30 and 70% depending on the solids concentration in the initial slurries and the width of the open interconnected macropores, which typically ranges between 20 and 100 µm in their smallest dimension and 50 to 500 µm in the largest one.

Phase analysis

The XRD pattern (Fig. 3) of sintered samples can be completely indexed with hydroxyapatite (JCPDS 09-0432), the only phase found present. No processing residue or secondary phases were found in the materials.

In the following sections, we describe the effect of the processing conditions, namely, the slurry concentration, freezing rate, and sintering, on the main properties of interest for tissue engineering: shape and amount of porosity, and compressive strength.

Influence of slurry concentration -Total porosity

By changing the initial slurry concentration, a broad control of the final porosity can be achieved (Fig. 4). The relationship between the final porosity in the lamellar zone and the initial slurry concentration is linear, provided sintering conditions are constant. Not only does the technique offer very reproducible results, but also the final porosity can easily be tuned by the initial slurry concentration.

-Porosity size

For highly porous samples (Fig 5a), the microstructure is almost cellular (bridges and walls have similar dimensions). Walls can be made very thin (<5m), especially for fast cooling rates (>5°C/min), and their thickness is very homogeneous throughout the whole sample. For lower porosity content (47%, Fig 5b), some scattering is observed in the distribution of wall thickness, as repelling particles from the ice front becomes more difficult with increasing particle content.

-Compressive strength

Compressive strength vs. total porosity is plotted in Fig. 6. Although for high porosity content (typically >60 vol%) the strength (16 MPa) is comparable to that reported in the literature, it increases rapidly when the porosity decreases, reaching 65 MPa at 56% porosity and 145 MPa at 47% porosity. Values obtained for these samples are well above those reported so far in the literature. In fact, the strength of the porous lamellar HAP is close to that of compact bone [START_REF] An | Mechanical properties of bone. Mechanical testing of bone and the boneimplant interface[END_REF].

Such high values allow considering the potential of these materials for some load-bearing biological applications.

Influence of cooling rate

-Porosity size Being highly anisotropic, the pore channels have been characterized here by two dimensional parameters: the long axis and the short axis. The pore size can be modified mostly by increasing or decreasing the cooling rate of the cold finger during freezing. For low cooling rate (<1°C/min), the final microstructure will be characterized by large lamellae thickness (>50 m) and pore width (>500 m in the long axis and 40 in the short axis, Fig. 7). Conversely, for fast cooling rate (>5°C/min), finer lamellae and pores are obtained.

-Compressive strength

A series of samples with identical slurry concentration (50 wt.%) have been processed using similar conditions, the only difference being the cooling rate during freezing. The total porosity is independent of the cooling rate and was the same for all the samples, as expected. It has previously been shown how the cooling rate strongly affects the porosity features (size and morphology). These changes in cooling rate have, consequently, a direct repercussion on the compressive strength (Fig. 8). For low cooling rate (<2°C/min), the compressive strength remains moderate, at 10 to 20 MPa. However, for fast cooling rates (>5°C/min), the compressive strength increases up to 60 MPa as the cooling rate increases. This effect can be interpreted by considering the effect on the microstructure previously illustrated in Fig. 7. Porosity size (i.e. defect size) decreases as the cooling rate increases, so that the strength also increases as the defects (pores) become smaller. This also explains why the rate of increase in strength diminishes for fast cooling rates (>5°C/min). Indeed, the relationship between porosity size and cooling rate is not linear (Fig. 7), so that it becomes more and more difficult to reduce the pore size by accelerating the freezing. Samples cooled at 7 and 10°C/min exhibit similar porosity dimensions, and therefore have similar compressive strength. The technical limits of the experimental setup do not allow cooling rates greater than 10°C/min, and biological response should be taken into consideration before trying to size down the porosity as much as possible.

Influence of sintering conditions

The sintering conditions are of prime importance for optimizing the physical and mechanical properties. Porosity and compressive strength were measured as a function of the sintering temperature (sintering of 3 hrs). Results are given in Fig. 9. When sintering temperature increases from 1250°C to 1325°C, porosity goes down and compressive strength goes up. Both then reach a plateau value, which means that densification is achieved at 1325°C. Further increase in sintering temperature will only increase the grain size. Therefore, the optimum sintering temperature was determined to be 1325°C.

It is also worth noting the increase of grain size with temperature is very limited, as compared to that observed with materials processed by conventional techniques (e.g. [START_REF] Benaqqa | Slow crack growth behaviour of hydroxyapatite ceramics[END_REF]). This particular behavior of freeze-dried materials was already observed for freeze-dried ceramics [START_REF] Lu | Preparation and mechanical properties of dense polycrystalline hydroxyapatite through freeze-drying[END_REF]. Although the underlying reasons of this behavior are still to be elucidated, it is probably partly related to the pinning of grain boundary at the surface during sintering of thin films, since the grain size lies in the same range of order as the lamellae thickness.

Discussion

Control of the microstructure and the physics of ice formation

The final microstructure is a replica of the ice. Modifying the amount of water in the slurry as well as the shape of the ice crystals will modify the final amount, shape and size of porosity of the porous ceramic scaffold. In particular, during the steady freezing regime, the ice crystals exhibit a homogeneous morphology throughout the whole sample; this explains why the lamellae thickness is very homogeneous throughout the whole sample. This also means that the ratio porosity thickness/layer thickness is primarily determined by the initial slurry concentration. For instance, if the volume fraction of ice and ceramic particles is the same, layer thickness (e.g. 20 microns) will be the same as the porosity (20 microns in their smallest dimension). For low slurry concentration, the porosity becomes predominant and layer thickness decreases.

The microstructural features can be controlled by applying the physics of ice formation. In particular, the ice-tip radius (and as a consequence the thickness of the ice crystals), which is physically determined by the magnitude of supercooling ahead of the freezing front [START_REF] Hunt | Pattern formation in solidification[END_REF], can be modified by increasing or decreasing the cooling rate during freezing. For fast cooling rates, supercooling becomes larger and the microstructure can be scaled down. On the other hand, under a very slow cooling regime the layer (or pore) thickness can be noticeably increased (Fig. 7).

Physical limits of the technique

Fundamentals of the physics of ice will determine the limits of the technique. In particular, the sequence leading to a steady lamellar ice front is well known and involves a progressive transition of the freezing front morphology from planar to cellular and cellular to lamellar (also called dendritic), due to the progressive build up of supercooling ahead of the freezing front. This evolution is replicated in the microstructure close to the initial freezing front, i.e. the one initially in contact with the cold finger, as seen in Fig. 10. If the zone corresponding to the planar to cellular zone (Fig. 1) is difficult to ascertain, the cellular-to-dendritic transition zone is clearly observed in the middle of the micrograph. Once the transition has occurred, a steady regime is reached, and the freezing front retains its morphology. Hence, two very different architectures are obtained in the final part. The lower part, formed during the cellular regime, exhibits a porous structure, but with poor pore interconnectivity. The average pore size is small; from a biomedical point of view, this will be of little interest and this part will probably have to be removed. The upper part exhibits the porous lamellar architecture, with open and interconnected porosity.

Although no further experiments were performed at this stage, the size of the transition zone seems highly dependant on cooling rate, slurry concentration and ceramic particle size.

Although already observed in freeze-dried polymeric materials (e.g. see fig. 2a in [START_REF] Madihally | Porous chitosan scaffolds for tissue engineering[END_REF] and fig. 1a in [START_REF] Zmora | Tailoring the pore architecture in 3-D alginate scaffolds by controlling the freezing regime during fabrication[END_REF]) the presence of such transition zones has not yet been explained. Combining the physics of ice and the interaction of the freezing front and inert particles, the present and previous experimental results can be interpreted in terms of modification of the freezing regime and transition of the freezing front morphology. The first zone in contact with the cold source exhibits faster freezing kinetics, and the equilibrium in terms of supercooling and heat diffusion has not been reached yet. Hence, the freezing interface exhibits an unstable cellular morphology, resulting in poor pore interconnectivity.

The influence of slurry concentration on the final porosity has been illustrated previously. The range in which porosity can be adjusted is physically limited by two phenomena. At low concentrations, even in the presence of a binder, collapse of the green body occurs during the sublimation of the ice. However, above a certain concentration (65wt.%), the interaction between particles in the slurry becomes more important and eventually particles cannot be repelled from the ice front anymore. As a result, the interconnected and open porous structure is lost (Fig. 11). Such architecture is not suitable anymore for tissue engineering. This critical slurry concentration is strongly dependent on the particle size. Earlier theoretical and experimental results have established the existence of a critical velocity for particles engulfment [START_REF] Asthana | Review -The engulfment of foreign particles by a freezing interface[END_REF][START_REF] Korber | Interaction of particles and a moving ice-liquid interface[END_REF]. If the velocity of the front is above the critical value, particles will be engulfed by the interface. In addition, it was found that the critical velocity is inversely proportional to the particle radius [START_REF] Korber | Interaction of particles and a moving ice-liquid interface[END_REF]. Repelling small particles from the ice front is easier than larger ones. Hence, by using a powder with a smaller granulometry (the particle size of the powder used here was 2.4 m, data provided by manufacturer), reaching lower porosity while retaining the interconnected open porosity becomes feasible.

Suitability as bone substitute

Klawitter and Hulbert [START_REF] Klawitter | Application of porous ceramics for the attachment of load bearing internal orthopedic applications[END_REF] established a minimum pore size of ~100 µm for bone growth into ceramic structures, and a similar conclusion was reached by Simske et al. [START_REF] Simske | Porous materials for bone engineering, in Porous Materials for Tissue Engineering[END_REF]. More recently, Itala claimed that bone ingrowth occurred in pores as small as 50 m [START_REF] Itala | Pore diameter of more than 100 mu m is not requisite for bone ingrowth in rabbits[END_REF], and similar results were obtained for porosity engineered in the range 15-40 m [START_REF] Dutta | Performance of degradable composite bone repair products made via threedimensional fabrication techniques[END_REF]. The porosity must be interconnected to allow the ingrowth of cells, vascularization and diffusion of nutrients. Typical porosities for HA scaffolds described in the literature range from 35% [START_REF] Pilliar | Porous calcium polyphosphate scaffolds for bone substitute applications -in vitro characterization[END_REF][START_REF] Simske | Porous materials for bone engineering, in Porous Materials for Tissue Engineering[END_REF] to 75% [START_REF] Tamai | Novel hydroxyapatite ceramics with an interconnective porous structure exhibit superior osteoconduction in vivo[END_REF] with pore sizes between 50 m [START_REF] Itala | Pore diameter of more than 100 mu m is not requisite for bone ingrowth in rabbits[END_REF][START_REF] Chang | Osteoconduction at porous hydroxyapatite with various pore configurations[END_REF] and 400 m [START_REF] Tamai | Novel hydroxyapatite ceramics with an interconnective porous structure exhibit superior osteoconduction in vivo[END_REF][START_REF] Tsuruga | Pore size of porous hydroxyapatite as the cell-substratum controls BMP-induced osteogenesis[END_REF][START_REF] Pilliar | Porous calcium polyphosphate scaffolds for bone substitute applications -in vitro characterization[END_REF][START_REF] Flautre | Porous HA ceramic for bone replacement: Role of the pores and interconnections -experimental study in the rabbit[END_REF]. A preponderance of pores this large can make the implant material rather weak. Moreover, even small movements of the implant (which are hard to avoid) can cause complications by cutting off blood supply to tissue in the pores, which can lead to inflammation. In addition, bone response reflects meso-scale structure, and direct bone attachment and conduction along the struts of porous scaffolds has been reported [START_REF] Zmora | Tailoring the pore architecture in 3-D alginate scaffolds by controlling the freezing regime during fabrication[END_REF] even for small porosity [START_REF] Dutta | Performance of degradable composite bone repair products made via threedimensional fabrication techniques[END_REF], revealing that pore connectivity and orientation are at least as important as pore size. Hence, the porosity range reported in this study might be suitable for osseous regeneration.

An additional might be taken into consideration, namely the dendrites surface relief covering the surface of the lamellar walls (Fig. 12). Considerable attention has been paid recently to the relationships between cellular response and morphological features such as mesostructure, patterns, roughness, etc. Particularly, this was made possible with the recent development of computer-assisted processing techniques. Results from a wide range of experimental conditions pointed out the beneficial effects of roughness, patterns and architectures, which seem to shorten the time necessary for new bone ingrowth and tissue formation [START_REF] Stokols | The fabrication and characterization of linearly oriented nerve guidance scaffolds for spinal cord injury[END_REF][START_REF] Zhu | Nanotopographical guidance of C6 glioma cell alignment and oriented growth[END_REF][START_REF] Yu | Guided cell adhesion and outgrowth in peptide-modified channels for neural tissue engineering[END_REF][START_REF] Recknor | Oriented astroglial cell growth on micropatterned polystyrene substrates[END_REF][START_REF] Patist | Freeze-dried poly(D,L-lactic acid) macroporous guidance scaffolds impregnated with brain-derived neurotrophic factor in the transected adult rat thoracic spinal cord[END_REF][START_REF] Miller | Oriented Schwann cell growth on micropatterned biodegradable polymer substrates[END_REF][START_REF] Manwaring | Contact guidance induced organization of extracellular matrix[END_REF][START_REF] Mahoney | The influence of microchannels on neurite growth and architecture[END_REF][START_REF] Anselme | Osteoblast adhesion on biomaterials[END_REF]. The presence of such patterns allows directional growth of the cells and faster vascularization of the implants. In our case, these features, originating from the dendritic shape of the growing ice crystals during freezing, are also oriented along the ice growth direction and hence the porosity. By adjusting the slurry concentration and freezing conditions, their size might be controlled to some extent, so as to match more favorable dimensions of guidance patterns [START_REF] Mahoney | The influence of microchannels on neurite growth and architecture[END_REF], improving the osteoconduction and osteoinduction characteristics. Experiments are under way to assess the validity of these hypotheses and the biological response of these materials.

It is finally worth mentioning that, because the process relies on physical and not chemical interactions, it can be easily extended to any calcium phosphate or suitable ceramic for biomedical applications. The unique features of the porosity obtained by this technique have already attracted attention for non-load bearing applications using polymeric materials. The technique can now be extended to load-bearing applications.

Conclusions

Based on an experimental study of freeze-drying of hydroxyapatite powders with various slurry concentrations and sintering conditions, the following conclusions can be made:

1. Porous scaffolds with total porosity ranging from at least 40% to 65% can be obtained by freezing of hydroxyapatite aqueous suspensions and subsequent ice sublimation and sintering.

The resultant porosity is open and unidirectional, exhibiting a lamellar morphology. Size of the porosity can be controlled by modifying the freezing rate of the slurries and the slurry concentration.

2. Formation of the porosity can be understood by applying basic principles of the physics of ice and the interaction of a moving solidifying interface (the freezing front) and inert particles in suspension. In particular, the influence of cooling rate, particle size and slurry concentration can be interpreted.

3. Due to their lamellar architecture and the pore shape anisotropy, the processed scaffolds exhibit unusually high compressive strength for such materials, i.e. up to 145 MPa for 47% porosity and 65 MPa for 56% porosity, allowing them to be considered for some load-bearing applications. [START_REF] Currey | Handbook of composites[END_REF][START_REF] Fung | Biomechanics. Mechanical properties of living tissues[END_REF] and literature values [START_REF] Almirall | Fabrication of low temperature macroporous hydroxyapatite scaffolds by foaming and hydrolysis of an alpha-TCP paste[END_REF][START_REF] Charriere | Hydroxyapatite cement scaffolds with controlled macroporosity: fabrication protocol and mechanical properties[END_REF][START_REF] Del Real | A new method to produce macropores in calcium phosphate cements[END_REF][START_REF] Kawata | Development of porous ceramics with well-controlled porosities and pore sizes from apatite fibers and their evaluations[END_REF][START_REF] Milosevski | Preparation and properties of dense and porous calcium phosphate[END_REF][START_REF] Ramay | Preparation of porous hydroxyapatite scaffolds by combination of the gel-casting and polymer sponge methods[END_REF][START_REF] Tian | Preparation of porous hydroxyapatite[END_REF][START_REF] Ioku | Preparation and characterization of porous apatite ceramics coated with beta-tricalcium phosphate[END_REF][START_REF] Landi | Carbonated hydroxyapatite as bone substitute[END_REF][START_REF] Liu | Influence of porosity and pore size on the compressive strength of porous hydroxyapatite ceramic[END_REF][START_REF] Ota | Preparation and compressive strength behavior of porous ceramics with beta-Ca(PO3)(2) fiber skeletons[END_REF][START_REF] Ramay | Biphasic calcium phosphate nanocomposite porous scaffolds for load-bearing bone tissue engineering[END_REF][START_REF] Sous | Cellular biocompatibility and resistance to compression of macroporous beta-tricalcium phosphate ceramics[END_REF][START_REF] Bignon | Optimization of the porous structure of calcium phosphate implants for bone substitutes and in situ release of active principles[END_REF][START_REF] Pilliar | Porous calcium polyphosphate scaffolds for bone substitute applications -in vitro characterization[END_REF][START_REF] Barralet | Preparation of macroporous calcium phosphate cement tissue engineering scaffold[END_REF][START_REF] Tancred | A synthetic bone implant macroscopically identical to cancellous bone[END_REF][START_REF] Dong | In vivo evaluation of a novel porous hydroxyapatite to sustain osteogenesis of transplanted bone marrow-derived osteoblastic cells[END_REF][START_REF] Chu | Mechanical and in vivo performance of hydroxyapatite implants with controlled architectures[END_REF] of porous HAP for tissue engineering, conventional and computer-assisted processing techniques. Each style of point corresponds to a different literature source. Fig. 7: Influence of the freezing rate over the pore width (56% total porosity).

Figures

  Figures

Fig. 1 :

 1 Fig. 1: Typical sample (left) and general microstructure of the porosity (right). Three distinctive zones are found, their dimensions depending largely upon processing conditions.

Fig. 2 :

 2 Fig. 2: Microstructure of porous HA samples with 64% porosity. Cross-sections parallel to the ice front (a) and detail (b), and perpendicular to the ice front (c) and detail (d). The slight difference in wall thickness (b and d) arises from the presence of a small gradient of density at the upper surface of the sample (a and b), due to particle segregation at the top of the mold during freezing.

Fig. 3 :

 3 Fig. 3: XRD spectrum, sintered materials.

Fig. 4 :

 4 Fig. 4: Total porosity in sintered materials vs. slurry concentration. The fit is plotted for the lamellar microstructures only.

Fig. 5 :

 5 Fig. 5: Effect of initial slurry concentration. Microstructure of porous HA samples, parallel to the ice front with (a) 64% porosity and (b) 47% porosity.

Fig. 6 :

 6 Fig. 6: Maximum compressive strength achieved by freeze casting (cooling rate of 5°C/min) vs. porosity. Comparison with compact bone[START_REF] Currey | Handbook of composites[END_REF][START_REF] Fung | Biomechanics. Mechanical properties of living tissues[END_REF] and literature values[START_REF] Almirall | Fabrication of low temperature macroporous hydroxyapatite scaffolds by foaming and hydrolysis of an alpha-TCP paste[END_REF][START_REF] Charriere | Hydroxyapatite cement scaffolds with controlled macroporosity: fabrication protocol and mechanical properties[END_REF][START_REF] Del Real | A new method to produce macropores in calcium phosphate cements[END_REF][START_REF] Kawata | Development of porous ceramics with well-controlled porosities and pore sizes from apatite fibers and their evaluations[END_REF][START_REF] Milosevski | Preparation and properties of dense and porous calcium phosphate[END_REF][START_REF] Ramay | Preparation of porous hydroxyapatite scaffolds by combination of the gel-casting and polymer sponge methods[END_REF][START_REF] Tian | Preparation of porous hydroxyapatite[END_REF][START_REF] Ioku | Preparation and characterization of porous apatite ceramics coated with beta-tricalcium phosphate[END_REF][START_REF] Landi | Carbonated hydroxyapatite as bone substitute[END_REF][START_REF] Liu | Influence of porosity and pore size on the compressive strength of porous hydroxyapatite ceramic[END_REF][START_REF] Ota | Preparation and compressive strength behavior of porous ceramics with beta-Ca(PO3)(2) fiber skeletons[END_REF][START_REF] Ramay | Biphasic calcium phosphate nanocomposite porous scaffolds for load-bearing bone tissue engineering[END_REF][START_REF] Sous | Cellular biocompatibility and resistance to compression of macroporous beta-tricalcium phosphate ceramics[END_REF][START_REF] Bignon | Optimization of the porous structure of calcium phosphate implants for bone substitutes and in situ release of active principles[END_REF][START_REF] Pilliar | Porous calcium polyphosphate scaffolds for bone substitute applications -in vitro characterization[END_REF][START_REF] Barralet | Preparation of macroporous calcium phosphate cement tissue engineering scaffold[END_REF][START_REF] Tancred | A synthetic bone implant macroscopically identical to cancellous bone[END_REF][START_REF] Dong | In vivo evaluation of a novel porous hydroxyapatite to sustain osteogenesis of transplanted bone marrow-derived osteoblastic cells[END_REF][START_REF] Chu | Mechanical and in vivo performance of hydroxyapatite implants with controlled architectures[END_REF] of porous HAP for tissue engineering, conventional and computer-assisted processing techniques. Each style of point corresponds to a different literature source.

Fig 8 :

 8 Fig 8: Compressive strength vs. cooling rate. Samples with 56% porosity. Sintering 4 hrs at 1300°C.

Fig. 9 :

 9 Fig. 9: Influence of sintering temperature (3 hrs at dwell temperature) on compressive strength, total porosity and grain size. Samples with 24 vol% of initial powder, frozen at 5°C/min.

Fig. 10 :

 10 Fig. 10: Microstructure of porous HA samples, close to the cold finger (cross section perpendicular to the ice front). The initial porosity gradient is related to the evolution of the freezing front morphology with time. The cellular to lamellar transition zone is clearly observed in the middle of the micrograph.

Fig 11 :

 11 Fig 11: Cross section parallel to ice front for initial slurry concentration of 42 vol%. The interconnected and open porous structure is completely lost. The particle-particle interaction becomes so important that particles cannot be repelled from the progressing ice front.

Fig 12 :

 12 Fig 12: Surface dendrites, oriented along the ice growth direction, cross section parallel (a) and perpendicular (b) to ice growth direction. These features might act as guidance pattern for cell growth, in addition to the already oriented porosity.

Acknowledgements

This work was supported by the National Institute of Health (NIH/NIDCR) under grant No. 5R01 DE015633 (Novel Scaffolds for Tissue Engineering and Bone-Like Composites).