F Bartl

H Klaus

R Dallner
email: rudolf.dallner@fh-ingolstadt.de

O Huber

Material behaviour of a cellular composite undergoing large deformations

Keywords: material behaviour, cellular composite, finite element method, large deformation, yield surface, multiaxial material tests, crashbox

The intention of this paper is to give an overview of the recent investigations concerning the yield behaviour of a new cellular composite material. This novel material is a syntactic foam with a porous mineral granulate embedded in polyamide. The aim of the investigations is to develop a yield surface suitable for numerical simulations.

To acquire the necessary data for the initial yield surface, a number of different experiments is shown. These experiments include quasi-static uniaxial tension, shear and compression tests as well as hydrostatic tests and experiments with blocked lateral deformations. The evolution of the yield surface into the range of large deformations is investigated and it can be shown that the specific plastic deformation energy is a suitable parameter for the hardening of the material.

To exemplify the capabilities of the material and the material model, a prototype of a crashbox is tested under quasi-static conditions and compared to numerical results.

T P I R C S U N A M D E T P E C C A 1 Introduction
A modern engineer is often confronted with the task of designing components with a good ability to absorb energy at almost constant stress together with

T P I R C S U N A M D E T P E C C A
sufficient strength. Cellular materials are one type of material which is suitable for such a component. This fact in combination with a low specific weight makes them a suitable material for designing components relevant for crash applications. The material investigated in this paper is a new cellular composite material. It is a syntactic foam with porous mineral granulates embedded in a cast polyamide matrix (see Fig. 1). The material properties of the investigated cellular composite can be adjusted to the specific use by using suitable matrix-and/or granulate materials [START_REF] Dallner | Final report of the research project "IMVAL[END_REF]. For the investigated material Dennert Poraver GmbH Schlüsselfeld, Germany, delivered the cellular glass granulates with a diameter between two and four millimeters. The composite is produced by Quadrant engineering plastic products Almelo, Netherlands. Fig. 2 shows the casting process schematically. First the granulates are filled in the form and the form is then shaken to obtain a dense sphere packing of the granulates. After that, the granulate fill is infiltrated with the matrix material.

T P I R C S U N A M D E T P E C C A
The mechanical properties of the cast polyamide PA6 (ρ = 1.15g/cm 3 ) and the glass foam granulates (ρ average = 0.37g/cm 3 ) are described in [START_REF] Klaus | Experimentelle und numerische Untersuchung der Verbundwerkstoffes "Mineralschaumgranulate in Polyamidmatrix" im elastischen Bereich[END_REF]. The mechanical properties of the polyamide (both stiffness and strength) are dependent on the water content and the loading velocity. Experimental investigations of the granulates showed a Young's modulus of E = 1700M P a and a Poisson's ratio of ν = 0.3. The structure of the composite and the distribution of the glass foam granulates within the matrix were investigated by micro computer tomography analysis (µ-CT analysis) [START_REF] Godehardt | Analyse tomographischer Aufnahmen von Polyamidmatrizen mit keramischen Granulaten[END_REF]. Fig. 3 shows the results of the µ-CT analysis of a block with the dimensions 66mm × 45mm × 45mm. The average cell diameter is consistent with the average granulate diameter. The µ-CT analysis detected a matrix volume fraction of 50.8% and a granulate volume fraction of 49.2%.

The elastic behaviour of the material was investigated in [START_REF] Klaus | Lightweight Potential of Novel Cellular Spherical Composites[END_REF][START_REF] Klaus | Materialverhalten eines Kugelverbundwerkstoffes im linear-viskoelastischen Bereich[END_REF][START_REF] Klaus | Charakterisierung, Modellierung und Anwendung zellularer Verbundwerkstoffe im elasto-plastischen Bereich[END_REF][START_REF] Klaus | FEM-Modellierung für den Verbundwerkstoff "Mineralschaumkugeln in Polyamidmatrix" im elastischen Bereich[END_REF]. It was shown that the specific stiffness and the specific compression strength under uniaxial compression are improved significantly by the cellular granulates up to levels of strains, where the material starts to yield.

In the following the experimental investigations of the elasto-plastic material behaviour of the cellular composite are described. Special emphasis is laid on the description of the yield surface and its evolution into the range of large deformations. The yield behaviour is incorporated into an existing material model [START_REF] Bartl | Simulation eines neuartigen Mineralschaum-Verbundwerkstoffes für Crash-Applikationen[END_REF], *MAT BILKHU DUBOIS FOAM based on [START_REF] Nusholtz | A simple elastoplastic model for foam materials[END_REF], in the FEM software LS-DYNA [START_REF] Hallquist | LS-DYNA Theoretical Manual[END_REF]. This is a prerequisite for numerical simulations, which are an essential part of a modern product development process. The focus of this paper is on quasi-static investigations. However, the strain rate sensitivity of the material may be accounted for, e.g. by appropriate scaling of the quasi-static yield surface [START_REF] Bartl | Ein numerisches Werkstoffmodell für Mineralschaumverbundwerkstoffe für crashrelevante Anwendungen[END_REF].

T P I R C S U N A M D E T P E C C A 2 Experimental investigations
For all quasi-static tests the cast polyamide was conditioned to 2.5% water content (percent by weight). The tests were carried out with an electromechanical universal material testing machine (Zwick; Z 150). For the axial strain measurement at the surface of the specimen an incremental automatic extensometer (Zwick, MultiXtens) was used. The motorized arms of the extensometer were controlled by the test software (Zwick, testXpert) and allowed a variable gauge length. The test software controlled the test itself and the data acquisition during the test.

The composite has a unique infiltration direction due to the production process (Fig. 2). Uniaxial compression tests in the elastic range were made to investigate the isotropic behaviour of the mechanical properties on the basis of the Young's modulus. Two cubes with an edge length of 70mm were used as specimens, which were loaded between two parallel plates with a constant velocity of 10mm/min. The extensometer had a gauge length of 50mm. One specimen was loaded consistent with the infiltration direction, the other one was tested perpendicular to the direction of the infiltration. The Young's modulus was equal in good approximation in both cases. Therefore the composite is considered as quasi-isotropic.

Fig. 4a shows typical stress strain diagrammes of the composite under uniaxial compression. The dimensions of the specimen were 50mm × 50mm × 100mm. The specimens were loaded between two parallel plates with a constant velocity ( εstart = 1 • 10 -3 s -1 ). The strain was measured with the automatic extensometer (gauge length = 50mm). The principle mechanical behaviour of the cellular composite under compression load corresponds to conventional cellular solids. A range of transition

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
exists between the quasi linear elastic range and the range with low hardening (Fig. 4a).

For the uniaxial tension test (Fig. 4b) the width of the specimen (70mm × 70mm × 250mm) was reduced in its central part over a length of 150mm to a cross section of 45mm × 45mm. The load was transfered to the specimen via form locking. The strain was measured with the automatic extensometer (gauge length = 60mm). The specimens were loaded with a constant velocity ( εstart = 1 • 10 -3 s -1 ). Fracture occurred at 6% strain. Under uniaxial tension the range of transition was more pronounced than in the uniaxial compression tests.

Fig. 5 shows the experimental setup for the shear test according to DIN 53 294.

The vertical angle of the shear specimen was about 4 • higher compared to the DIN. Thus the distribution of the shear stress along the specimen is supposed to be more homogeneous [START_REF] Klaus | Charakterisierung, Modellierung und Anwendung zellularer Verbundwerkstoffe im elasto-plastischen Bereich[END_REF][START_REF] Huber | Semester: Konstruktion von Versuchsvorrichtungen für die Werkstoffprüfung eines neuartigen Leichtbauverbundwerkstoffes[END_REF] (Fig. 5b). The compression forces were transfered through cutting plates at the ends of the specimen (Fig. 5a). The dimensions of the specimen were 20mm × 50mm × 200mm. The specimens were loaded with a constant velocity ( γstart = 1 • 10 -3 s -1 ). The displacement between the two mounting plates was measured with an incremental measuring sensor (Heidenhain, ST3088). Fig. 6 shows the corresponding shear stress strain curves. Fig. 7 illustrates the failure behaviour of the composite. Under uniaxial compression the polyamide cell walls buckle in an elastic-plastic way and the composite shows a ductile behaviour, whereas the polyamide cell walls break under uniaxial tension at relatively low strains (Fig. 4b,7b). Pure shear load consists of compression and tension load in the principle axis (Fig. 7c). Therefore the Cellular materials as well as cellular composites compact under compression in a different way than materials without porosity. The yield criterion depends on hydrostatic pressure. Therefore a hydrostatic pressure test is essential to investigate this effect on the yield criterion of cellular materials. Loading the specimen through a surrounding fluid is an easy way to create hydrostatic compression. Fig. 8a shows the experimental setup schematically.

T P I R C S U N A M D E T P E C C A
Glycerin was used as fluid. The end of the stamp, reaching into the fluid, was conical. Thus the included air spilled out of the pressure chamber with the displaced fluid until the o-ring entered the cylindrical part of the pressure chamber and started to seal. The feed opening and the opening for the stamp A wrapping protected the porous specimen against infiltration with the fluid. The wrapping was flexible to ensure the hydrostatic loading of the specimen. It consisted of two layers due to the high pressure. The first layer was a flexible heat shrink tube and the second a butyl rubber membrane (Fig. 8c and8d). The glass granulates on the boundary of the specimen were sliced. These granulates were loaded by the wrapping directly and after break down the flexible wrapping would have infiltrated in the emerging cavity (Fig. 8b). Therefore the sliced granulates were replaced by epoxy resin to avoid this effect.

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
The volumetric strain of the specimen was directly proportional to the penetration depth of the stamp, which was measured by an extensometer (Fig. 8a).

A reference test was required to consider the compliance of the pressure chamber, the sealing, the fluid and the wrapping. For this a steel specimen was used, the stiffness of which is much greater than the stiffness of the cellular composite. The reference specimen had the same dimensions and wrapping as the cellular composite. The volumetric strain, which was measured in the reference test, had to be subtracted from the volumetric strain of the cellular composite

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
at the corresponding pressure loads. Fig. 8e shows the experimental results. Under hydrostatic pressure the material strength is higher and the slope of the plateau is obviously steeper compared to the uniaxial compression test.

The experimental results are the basis for the description of a yield surface. The initial yield surface (Fig. 9) is determined for a plastic strain of 0.1%. For the description of the yield function the first invariant of the stress tensor I 1 = σ ii and the second invariant of the stress deviator J 2 = 1 2 s ij s ij are used. The yield surface is not symmetrical to the √ 3J 2 -axis (von Mises stress). Therefore a yield function [START_REF] Mahrenholtz | Zum elastisch-plastischen Übergangsverhalten metallischer Werkstoffe[END_REF] with the form

F = f (I 1 , J 2 ) -k 2 = a 0 • I 2 1 + a 1 • I 1 + J 2 -k 2 = 0 (1)
is necessary for a good approximation of the experiments (Fig. 9). If the yield function would contain only I 2 1 (without I 1 ), the elliptical function would be symmetrical to the √ 3J 2 -axis [START_REF] Deshpande | Isotropic constitutive models for metallic foams[END_REF], which is not sufficient for the material behaviour under consideration. The material law *MAT BILKHU DUBOIS FOAM in the commercial FEM software LS-DYNA, based on [START_REF] Nusholtz | A simple elastoplastic model for foam materials[END_REF], uses a yield function like equation (1). To investigate the strain rate sensitivity of the material, additional uniaxial compression tests were carried out. The experiments were conducted on an Instron VHS High Strain Rate System (VHS 50/20), a test rig which is capable of doing tests at velocities up to 20m/s (Fig. 10a). The test rig and data acquisition was controlled by the test software (High Rate, Instron). For the test at the strain rate ε = 0.01s -1 a sampling rate of 50Hz was used, for the strain rates ε = 1s -1 and ε = 10s -1 data was acquired with a sampling rate of 20kHz and the test at the strain rate ε = 100s -1 was sampled with a rate of 2MHz. It has to be noted, that due to the impact at high velocities of the piston onto the crosshead, vibrations are induced into the test rig. These vibrations will result in an oscillation of the recorded force signal. Therefore, to acquire reasonable data, suitable filtering techniques have to be used. In this case the data of the experiment with the strain rate ε = 100s -1 was filtered (digital lowpass filter, type Butterworth, threshold frequency 4kHz). Strain rate hardening can be observed (Fig. 10b). The stress level of the material is increasing with testing speed. This may be caused by the strain rate sensitivity of the polymeric matrix material 1 .

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

Experimental data at large deformations

In order to investigate the development of the yield surface under large deformations, the experimental data has to be expanded until the beginning of the compaction phase [START_REF] Bartl | Experimentelle und numerische Untersuchung ausgewaehlter zellularer Werkstoffe bei höheren Dehnraten[END_REF] (see Fig. 12).

To acquire additional points of the yield surface, a testing device was used, which blocked the transverse strain of the specimen and measured the resulting transverse force as well as the longitudinal force. The specimen (35mm × 35mm × 35mm) was fitted between four moveable rams (Fig. 11a) which constrained the lateral deformation. To measure the transverse force strain gauges were applied on the bolts connecting the rams to the outer walls. The testing device was placed in an electro-mechanical universal material testing machine (Zwick; Z 250), which provided the longitudinal loading of the specimen with a punch (Fig. 11b). In the following this experiment will be called 1 As it was impossible to install a extensometer, like the Zwick MultiXtens to the VHS, the displacements were measured by the piston movement. Although the piston movement is a suitable means of measuring the displacements in the elastoplastic range of the material, it lacks the needed accuracy for the measurement of the small displacements in the elastic range of the material. Therefore the strain rate influence on the quasi elastic properties of the material will not be discussed in this paper. As a result of the blocked lateral deformation the longitudinal stress of the quasi triaxial experiment (mean values of four specimens) is increased compared to the uniaxial experiment (Fig. 12a). In addition a significant stress in the transverse direction was measured. This is a result of the non-zero Poisson's ratio of the composite material within the plastic range. Fig. 12b shows a comparison between the quasi triaxial experiment (both lateral deformations were blocked) and the quasi biaxial experiment (only one lateral deformation was blocked, here mean values of two experiments). The stresses of the biaxial test are lower, because only one lateral deformation was constrained. The results in Fig. 12 are plotted over volumetric strain ε V (with: ε V = -ln V V 0 ), according to [START_REF] Nusholtz | A simple elastoplastic model for foam materials[END_REF].

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
These tests can be performed up to volumetric strains ε V higher than 0.6, that is into the compaction phase of the material. For the description of the

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
yield surface it would be beneficial, if the data for the hydrostatic test were also available for these large strains. As the hydrostatic test rig was not fit for the corresponding high pressures, a numerical FEM model of the material on the mesoscopic level was developed [START_REF] Huber | Herstellung und Eigenschaften syntaktischer Metallschäume mit unterschiedlichen Matrix-und Füllmaterialien (Teil I -III)[END_REF][START_REF] Bartl | Simulation eines neuartigen Mineralschaum-Verbundwerkstoffes für Crash-Applikationen[END_REF] (Fig. 13b) in order to simulate the test. Due to the low strength of the granulate under investigation, it has only a small influence on the material behaviour at large strains [START_REF] Bartl | Simulation eines neuartigen Mineralschaum-Verbundwerkstoffes für Crash-Applikationen[END_REF] and therefore was neglected in this model. The polyamide matrix was modeled with 3D tetrahedral elements using a piecewise linear elasto-plastic material behaviour and an element erasing criterion. The FEM programme LS-DYNA was used for the simulation. The mesoscopic model shows a similar behaviour as the experimental curve (Fig. 13a), available up to volumetric strains of 0.3. The offset of the two curves may be due to the simplifications in the numerical model. Nevertheless the general behaviour is the same and it may be found appropriate to extrapolate the experimental data into the range of larger strains by means of adding an offset to the numerical data (Fig. 13a). The results gathered in the uniaxial tension, the uniaxial compression, the quasi biaxial compression, the quasi triaxial compression and the hydrostatic pressure tests are plotted into the invariant space (von Mises equivalent stress versus mean stress, Fig. 14). Each of the symbols represents an experiment (or data obtained by numerical simulation)2 at a certain level of deformation. For this diagramme the individual experiments are compared at levels of the same volumetric strain. The initial yield surface shown in Fig. 9 is also drawn.

The material shows a pronounced hardening behaviour, that is with increasing volumetric strain, the corresponding absolute values of the mean stress and the von Mises equivalent stress increase significantly. Although the modeling on the mesoscopic level of the material is possible [START_REF] Bartl | Simulation eines neuartigen Mineralschaum-Verbundwerkstoffes für Crash-Applikationen[END_REF], it is only suitable for basic research of the material, because it would lead to very large FEM-models for larger components, which are almost impossible to handle even on modern computers. Therefore a model at the macroscopic level is needed, which represents the material as a continuum for a specific representative element volume and is suitable for application in standard FEM codes. In this paper the focus is on the yield surface of such a material model. For further information concerning material modeling the reader is referred to literature (e.g. [START_REF] Nusholtz | A simple elastoplastic model for foam materials[END_REF][START_REF] Deshpande | Isotropic constitutive models for metallic foams[END_REF][START_REF] Chen | Plasticity for Structural Engineers[END_REF][START_REF] Simo | Computational Inelasticity[END_REF]). Due to the different diameters and the purely statistical distribution of the granulates during the manufacturing of the material, an isotropic behaviour of the material is assumed in the following, which could be proved by experiment [START_REF] Dallner | Final report of the research project "IMVAL[END_REF].

T P I R C S U N A M D E T P E C C A
In literature different shapes of the yield surface have been proposed [START_REF] Deshpande | Isotropic constitutive models for metallic foams[END_REF][START_REF] Hanssen | Validation of constitutive models applicable to aluminum foams[END_REF][START_REF] Chen | A phenomenological framework of constitutive modelling for incompressible and compressible elasto-plastic solids[END_REF][START_REF] Wicklein | Numerical investigations of the elastic and plastic behaviour of an open-cell aluminium foam[END_REF]. The most common formulation is that of an ellipsoid in the principal stress state [START_REF] Öchsner | Experimentelle und numerische Untersuchung des elastoplastischen Verhaltens zellularer Modellwerkstoffe[END_REF]. An elliptical yield surface can be described with an equation like [START_REF] Dallner | Final report of the research project "IMVAL[END_REF]. It can be seen in Fig. 14 that for this material a fitting of the yield curve at higher volumetric strains is not possible with an elliptic curve.

Different approaches may be used to compare states of multiaxial stress on the yield curve. In case of this specific material, a better fitting to elliptical curves is achieved by using the specific plastic deformation energy as a criterion. The increment of the specific plastic deformation energy can be calculated as

dw p = dε p ij σ ij (2)
T P I R C S U N A M D E T P E C C A
where dε p ij denotes the increment of the plastic strain tensor and σ ij represents the current stress tensor. The increment of the plastic strain tensor can be calculated by

dε p ij = dε ij -dε e ij . (3)
In this equation dε ij denotes the increment of the total strain tensor and dε e ij represents the increment of the elastic strain tensor, which can be determined by

dε e ij = 1 + ν e E • dσ ij - ν e E • dσ kk • δ ij . (4)
E is Young's modulus, ν e is the elastic Poisson's ratio and δ ij denotes the Kronecker symbol. Due to the damage to the granulates and the matrix material in the ongoing deformation Young's modulus is not a constant value, but changes. To investigate the quantity of the change, uniaxial experiments where conducted, in which the specimen was loaded and unloaded repeatedly. Due to the inhomogeneous structure of the material, a superposition of different effects results in the behaviour of the material shown in Fig. 15a. The cellular build-up of the material leads to bending and buckling of the cell walls during the loading process. The stresses induced into the matrix material by these deformations result in an elastic relaxation during the unloading process. These effects are overlapped by the viscoelastic behaviour of the matrix and the damage of the granulates. To separate these effects and evaluate the influence on the elastic behaviour, further investigations have to be conducted. In this paper, the change of the Young's modulus of the compressed material is of relevance. Therefore, the following approach was chosen: A modulus was approximated to the elastic unloading path of the specimen (Fig. 15a). The approximated modulus will be called current Young's modulus in the following. Fig. 15b shows the results of this investigation in values of relative Young's modulus. The relative Young's modulus in this case is the relation of the current Young's modulus E to the initial Young's modulus E 0 (E rel = E E 0 ). It can be observed, that at first, E rel falls, which corresponds to the loss of stiffness due to the break of the granulates. After the current Young's modulus has reached approximately 62% of the original Young's modulus value, E rel starts to increase to values which are even higher than the initial value (E rel > 1). This seems to be plausible, since the voids of the materials are mostly closed and the material behaves more and more like a solid matrix.

In Fig. 16 the experimental data gathered in section 3 is plotted at constant values of specific plastic deformation energy obtained by equation [START_REF] Huber | Herstellung und Eigenschaften syntaktischer Metallschäume mit unterschiedlichen Matrix-und Füllmaterialien (Teil I -III)[END_REF]. The initial yield surface, described in section 2, is plotted in the diagramme as well. The ellipses are calculated using equation [START_REF] Dallner | Final report of the research project "IMVAL[END_REF] and are fitted to the uniaxial, quasi biaxial, quasi triaxial and hydrostatic compression test. For the initial yield surface and the yield surface at w p = 0.25M P a the uniaxial tension test was used as well. In the uniaxial tension test, the material will fail at w p ≈ 0.28M P a and thus the yield surfaces for w p > 0.25M P a are not extended into the range of positive mean stress.

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
Hence, for this specific material the yield condition is given by equation (1).

To account for the work hardening of the material, the material parameters a 0 , a 1 and k are defined by the evolution of the yield surface, and are dependent on w p , as described in Fig. 16. Therefore, equation (1) can be written as

F = a 0 (w p ) • I 2 1 + a 1 (w p ) • I 1 + J 2 -k 2 (w p ) = 0 . (5)
T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
The increment of the plastic strain dε p ij can be calculated using the flow rule

dε p ij = dγ • ∂g ∂σ ij . (6)
In this equation dγ is the proportionality factor and g is the plastic potential function. For a cellular material, such as the presented composite material, the plastic potential function

g = 3J 2 + 1 2 I 2 1 (7)
was selected according to [START_REF] Bilkhu | Material modeling of Structural Foams in Finite Element Analysis Using Compressive Uniaxial and Triaxial Data[END_REF], which results in a nonassociated flow rule.

To represent the failure of the material a failure criterion for the tensile domain as shown in [START_REF] Bartl | Ein numerisches Werkstoffmodell für Mineralschaumverbundwerkstoffe für crashrelevante Anwendungen[END_REF] may be added to the material model3 . The elements are deleted from the FEM model, when a certain level of the first principal strain is reached. Due to the fact, that failure happens at different stages of deformation in the tension and the compression domain, the strain state has to be determined by checking whether the volumetric strain is positive or negative.

This continuum material model was implemented into the FEM code LS-DYNA (version 970, revision 5434a) as a user defined material routine on the basis of the *MAT BILKHU DUBOIS FOAM material model. The effect of the variation of Young's modulus was neglected in the continuum model, because it is supposed that its influence is small at large plastic deformations, for which this model was derived. In order to show the improvement of the result quality in simulations of the described material using the specific plastic deformation energy, the quasi triaxial test was modeled. In the first simulation, the described material model was used for the specimen. In the second simulation *MAT BILKHU DUBOIS FOAM in its original version was used, which uses the volumetric strain as hardening criterion (see Fig. 14).

The described material model is capable of modeling the complex multiaxial state of stress in the quasi triaxial test (Fig. 17a). This yields in a reasonable fitting of experiment and simulation both in longitudinal and transversal direction. The simulation using the volumetric strain as hardening criterion (Fig. 17b) shows a different result. As would be expected from Fig. 14, experiment and simulation are matching well as long as the deformation of the specimen is moderate. The results are starting to deviate, when a certain level of volumetric strain (in this case ε V ≈ 0.4) is reached (Fig. 17b). Repeated unloading allowed the documentation of deformation during the uniaxial compression test. The unloading curves are not included completely in Fig. 19b to improve the clearness. The force-displacement characteristic is relatively homogeneous up to a compression of 100mm. The force variation along the displacement axis corresponds clearly with the buckling characteristics of the crashbox (Fig. 19b). The cellular composite core buckled periodically. The

T P I R C S U N A M D E T P E C C A a) b)
T P I R C S U N A M D E T P E C C A
fracture of the cellular composite started at the tension side of the buckle. The coherence of the core material did not fail completely during the test. At the beginning of the deformation, a force level between approximately 50kN and 75kN was measured. At approximately 90mm deformation, the inner buckles of the core had contact. From this point on, the force displacement curve rose at a progressive rate. The crashbox with the composite core had a weight of 939g. A pure steel crashbox with the same weight would have a sheet thickness of 1.8mm. Because a sheet with the desired thickness was not available, two crashboxes with a sheet thickness of 1.5 mm and 2.0mm respectively were investigated. For a comparison of the different crashboxes the specific energy absorption e = s s0 =0 F ds m crashbox [START_REF] Klaus | FEM-Modellierung für den Verbundwerkstoff "Mineralschaumkugeln in Polyamidmatrix" im elastischen Bereich[END_REF] and the crush force efficiency

η = s s0 =0 F ds F max (s -s0 ) (9)
according to [START_REF] Hanssen | Static crushing of square aluminium extrusions with aluminium foam filler[END_REF] were used. Fig. 20 shows the comparison of the crashboxes with the crush force efficiencies for s = 100mm. For the pure steel crashboxes local force maxima were selected for this special comparison, because the first peak could be reduced by various techniques, e. g. corrugations. The

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
crashbox with a core of cellular composite showed the highest specific energy absorption, 6.64kJ/kg in relation to 5.36kJ/kg for the 2mm steel shell. Due to the pronounced force peaks of the pure steel shell solutions, the crush force efficiency of the crashbox with a core was better, 0.81 in relation to 0.69 and 0.59 respectively. An efficiency of 1 would be the ideal solution resulting in a constant force displacement diagramme. A FEM model of the described crashbox was developed for simulation using HyperMesh (version 8.0) as a pre-processor and solved with LS-DYNA (version 970, revision 5434a with the described user defined material model). The undeformed and deformed shapes are shown in Fig. 21. Because of the quasistatic conditions of the experiment, only the quasi-static parameters of the materials were used. Because of the explicit formulation of the used solver, it has to be ensured that the kinetic energy is negligible compared to the total energy of the model for a suitable representation of the quasi-static test. Fig. 21a shows the used FEM model of the crashbox. For the steel shell, the material model *MAT PIECEWISE LINEAR PLASTICITY and Belytschko-Tsay shell elements with a size of approximately 3mm were used. The core was modeled using the described material model and solid hexaedral elements with one integration point. Due to the under integrated solid formulation of the elements, a hourglass control using the Flanagan-Belytschko stiffness approach was used. A comparison of the force-displacement curve (Fig. 22) shows similar load levels. To improve the repeatability of the experiments and correlation between experiment and simulation, corrugations in the shell to trigger the deformation are required. Note: The onset of cracks in the core on the macroscopic level begins at deformations of about 50 -60mm (see also Fig. 19). Therefore the tensile failure of the material can be neglected up to this point. For simulations beyond that point, a failure model, as described in [START_REF] Bartl | Ein numerisches Werkstoffmodell für Mineralschaumverbundwerkstoffe für crashrelevante Anwendungen[END_REF] has to be implemented. The investigations of the crashbox were conducted in a quasi-static environment, as the focus of this paper is on the development of a yield surface of the material and not on other effects like strain rate sensitivity. In a test under higher strain rates like the AZT Test ( [START_REF]AZT Crashreparaturtest Front (Neuer RCAR Strukturtest -10[END_REF], impact velocity: 16km/h), strain rate effects (see Fig. 10) should be included, as done e.g. in [START_REF] Bartl | Ein numerisches Werkstoffmodell für Mineralschaumverbundwerkstoffe für crashrelevante Anwendungen[END_REF], where the yield surface is scaled by appropriate factors, corresponding to the current strain rate.

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

ARTICLE IN PRESS 6 Conclusion

A novel cellular material made of porous glass granulates embedded in a polyamide matrix has been investigated in this paper.

The elastic-plastic behaviour of the material up to large deformations was investigated by various experiments. The multi-axial state of stress plays a major role in the onset of plastic yielding. E.g. the blocking of lateral deformations during uniaxial compression results in significant increase of the compressive stress and in stresses in the lateral direction as well. This means that lateral contraction must not be neglected for this type of materials. The yield surface can be represented with an elliptical yield function. This will hold for larger deformations as well, if the specific plastic deformation energy is used as the evolution parameter for the yield surface.

The representation of the yield surface has been implemented in a material model for FEM simulations. It has been shown that the behaviour under a multi-axial state of stress can be simulated.

Using a crashbox under quasi-static loading as an example it could be shown that the cellular material may be beneficial for such applications and that FEM simulation of component behaviour is possible with the proposed material representation.

Besides this, it would be beneficial to incorporate into the material model additional effects like strain rate sensitivity, material failure, variation of Young's Modulus and Poisson's ratio over deformation history, which were not the scope of this paper. Part of this was shown in previous papers [START_REF] Bartl | Ein numerisches Werkstoffmodell für Mineralschaumverbundwerkstoffe für crashrelevante Anwendungen[END_REF][START_REF] Dallner | Final report of the research project "IMVAL[END_REF] or is the scope of current research.

Fig. 1 .

 1 Fig. 1. Cellular composite material: a) Granulates in Matrix; b) Cross Section of a single granulate [2]

Fig. 2 .

 2 Fig. 2. Composite production -schematically : a) filling of the form with granulates; b) compaction by vibration; c) infiltration of the form; d) cellular composite

Fig. 3 .

 3 Fig. 3. Micro computer tomography analysis of the composite: a) Average cell diameter; b) Detail of the 3D structure -upper picture: matrix -lower picture: granulates

Fig. 4 .

 4 Fig. 4. Quasi-static uniaxial tests: a) compression test; b) tension test

Fig. 5 .

 5 Fig. 5. Shear test: a) Experimental setup; b) Numerical Comparison between the test according to DIN and the modified test[START_REF] Klaus | Charakterisierung, Modellierung und Anwendung zellularer Verbundwerkstoffe im elasto-plastischen Bereich[END_REF][START_REF] Huber | Semester: Konstruktion von Versuchsvorrichtungen für die Werkstoffprüfung eines neuartigen Leichtbauverbundwerkstoffes[END_REF]

Fig. 6 .Fig. 7 .

 67 Fig. 6. Shear test: Stress strain curves material behaviour under uniaxial tension determines the shear behaviour.The shear strength corresponds approximately to the tensile strength (Fig.4b, 6).

Fig. 8 .

 8 Fig. 8. Hydrostatic pressure test: a) Experimental setup schematically; b) Specimen surface beneath the wrapping; c) First layer of the wrapping; d) Second layer of the wrapping; d) Hydrostatic pressure volumetric strain curves

Fig. 9 .

 9 Fig. 9. Initial yield surface of the celluar composite with experimental results

Fig. 10 .

 10 Fig. 10. Strain rate sensitivity a) Test rig (Instron VHS 50/20) b) Stress strain diagramme

Fig. 11 .Fig. 12

 1112 Fig. 11. Testing device for multiaxial experiments: a) Schematic view; b) Device in universal testing machinequasi triaxial. Two of the four rams could be dismantled. With this configuration the lateral deformation in one direction was left free, while the other one was blocked. This test setup will be called quasi biaxial test in the following.

Fig. 13 .

 13 Fig. 13. a) Hydrostatic data extrapolated with a numerical model; b) Numerical model

Fig. 14 .

 14 Fig. 14. Results in invariant space 4 Derivation of a continuum model for large deformations

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Investigation of compression damage: a) Uniaxial cyclic loading and unloading; b) Relative Young's modulus

Fig. 17 .

 17 Fig. 17. Simulation of the quasi triaxial experiment: a) Described material model b) *MAT BILKHU DUBOIS FOAM (original version) 5 Application: A crashbox under quasi-static compression load

Fig. 18 .

 18 Fig. 18. Stress strain curve of the steel DC01

Fig. 19 .

 19 Fig. 19. Crashbox -uniaxial compression test: a) crashbox dimensions; b) force-displacement curve of the quasi-static test

Fig. 20 .

 20 Fig. 20. Force versus displacement diagramme of the crashbox with and without the core out of cellular composite

Fig. 21 .

 21 Fig. 21. Crashbox: a) Before and b) After compression (displacement ≈ 60mm)

Fig. 22 .

 22 Fig. 22. Force versus displacement diagramme of the crashbox

All data points are evaluated from experimental data, only the two values for hydrostatic compression at volumetric strains of 0.4 and 0.5 are extrapolated data

The failure criterion in Fig.16is a basic illustration and not an exact graphical representation of the described procedures.

Acknowledgements

The presented investigations were made within the project IMVAL (Innovative Mineralschaum Verbund Applikationen für den Leichtbau). The authors thank the German "Bundesministerium für Wirtschaft und Technologie" for the encouragement of the project in the framework of the promotion programme InnoNet and all the partners for their specific support.


Frans Vainio

Timo Knuutila
email: timo.knuutila@it.utu.fi

Michael Maier

Esa Alhoniemi

Mika Johnsson

Olli S Nevalainen

Estimating PCB Assembly Times Using Neural Networks

Keywords: PCB ASSEMBLY, PRODUCTION PLANNING, NEURAL NETWORK APPLICATIONS, placement machines, assembly time PCB assembly, placement machines, assembly time, production planning, neural networks, multilayer perceptron, Bayesian framework

come   L'archive ouverte pluridisciplinaire

Introduction

Automated placement of electronics components on printed circuit boards (PCBs) is performed by various placement machines of different types. Each machine type is able to handle a certain set of component types. The components have different characteristics (e.g. dimensions) and this is reflected in the design of placement machines. Different mechanical solutions lead to differed nominal operation speeds, and the same machine can be run at different speeds according to the handled components types.

Production of PCBs is performed mainly in assembly lines consisting of several placement machines of different types connected to each other by the means of PCB conveyors. In addition to placement machines, a production line includes several other production machines, like glue dispenser, reflow-oven and special type component feeding systems. Assembly lines of such a high technology are very expensive and need to be utilized as efficiently as possible.

The aim of line balancing is to optimize the operations of the machines so that the number of PCBs processed in a time unit (throughput) is maximized, see [START_REF] Crama | Production planning problems in printed circuit board assembly[END_REF] and [START_REF] Smed | Techniques and applications of production planning in electronics manufacturing systems[END_REF]. This task is usually accomplished by maximizing the throughput of the bottleneck machine, i.e. the machine consuming the longest production time per PCB. Optimization of whole assembly lines is too complex task to be performed using one-level approaches as presented by [START_REF] Smed | Techniques and applications of production planning in electronics manufacturing systems[END_REF]. However, the complex optimization problems are connected to each other in such a way that they can be divided in a hierarchy of simpler ones. Since the single machine optimization is at the lowest level of this hierarchy, it has to be solved each time (for each different PCB type or for each group of PCBs of similar type) when higher level optimizations are made. Consequently rapid estimations, instead of time consuming optimization calculations of PCB assembly times, are needed to form the basis for higher level optimizations.

As first approach, the production time can be estimated from the nominal (empirical) component placement time. The total production time is the sum of a constant term for feeding to and fixing the next PCB on the placement position and from a term which is a product of the number of components and the average placement time per component. Benefits of this method are simplicity and immediate results. On the other hand, its accuracy is modest and it suits only for finding an initial approximate solution for line balance optimization. When a control software (i.e., result of the optimization) for an assembly robot has been implemented and tested in production, it becomes soon clear that there are situations where the above "rule of thumb" time estimate is not accurate enough. The optimization based on nominal component times may generate a low value for the objective (time) function but the observed real assembly times can be far from that and from optimal.

To get more detailed and accurate results, assembly engineering optimizer software is in broad use. Commercial software packages for assembly engineering optimization (direct from machine producers or from third party providers) are creating assembly lines of specific characteristics and choosing machines for selected component groups from a large collection of machine types. Optimizers obtain the particular assembly times from a built-in simulator for each particular machine type. The simulators use, among other things, information about the placement sequence, operation sequence, machine timings, component coordinates, and machine geometry.

Current simulators for assembly machines are relatively accurate; the usual error marginal of the placement machine simulators is less than 2 %. On the downside, a simulator does complex computations which take a substantial running time if excessive (feeder allocation and component placement sequence) optimization is needed. This limits the use of optimization in (multi-product) line balancing and production scheduling situations where results are needed in a short time. On the other hand, the accuracy of the estimates of the placement times is vital for the success of line balancing and scheduling; even the best possible optimizer is useless if its objective function is wrong or inaccurate.

Above considerations lead to the following question: How can we get rapid and accurate estimates on PCB assembly times? As a compromise between the two above methods for time estimates, multivariate linear regression model was derived for a rotary turret machine in [START_REF] Laakso | Estimating the Production Times in PCB Assembly[END_REF]. This model is based on the number of component types, the total number of placements and the dimensions of the PCBs. The model achieved good results in spite of its underlying hypothesis of linear dependency of time on the parameters. It is therefore tempting to consider the power of non-linear time estimates to these data.

In this paper, we consider multilayer feedforward networks of type MLP for estimating the assembly times. We train MLPs to approximate the objective (assembly time) function using Bayesian regulation with two very different machine types. For a turret type machine we use the same simulated data as in [START_REF] Laakso | Estimating the Production Times in PCB Assembly[END_REF]. The training data for gantry type machine is based on real production obtained from two different PCB assembly factories.

We begin in Section 2. with literature review. Section 3 contains a short description of rotary turret type machines and gantry type machines. Two commercial optimizer-simulators Valor [START_REF] Trilogy | Commercial assembly engineering software package[END_REF] and [START_REF] Siemens | Commercial assembly engineering software package[END_REF] which were used to evaluate the "true" assembly times for MLP training and testing are reviewed in Section 4. Non-linear estimation of assembly times by using MLP networks is presented in Section 5. Experiments with the MLP network are reported in Section 6. Concluding remarks appear in Section 7.

Literature review

Over the recent years, a rich body of papers have been published dealing with problems of optimizing manufacturing processes and scheduling operations in PCB assembly. For a comprehensive literature review on general manufacturing operations, see for instance [START_REF] Pinedo | Planning and scheduling in manufacturing and services[END_REF] and specially on PCB assembly operations, see [START_REF] Crama | Production planning problems in printed circuit board assembly[END_REF] and [START_REF] Smed | Techniques and applications of production planning in electronics manufacturing systems[END_REF]. Optimization of the feeder setup and component pick-and-place sequence is very important for the efficient utilization of placement machines and lines. The component pick-and-place sequencing problem is a travelling salesman problem (TSP), which is strongly NP-hard. As a result og this, most practical problem instances can not be solved to optimality in a reasonable time. Therefore, a common practice is to operate with suboptimal machine control programs as produced by suitable heuristic methods. These heuristic algorithms can generate good solutions at a reasonable computational time, but their drawback is that they are not able to guarantee optimality. For a comprehensive literature review on heuristic methods used in assembly optimization, see for instance Reeves, C.R., (1995), [START_REF] Ayob | Optimization of Surface Mount Placement Machine in Printed Circuit Board Assembly[END_REF], [START_REF] Ayob | Optimization for Surface Mount Placement Machines[END_REF][START_REF] Ayob | A survey of surface mount device placement machine optimization: Machine classification[END_REF].

In the everyday industrial practice planning of assembly processes is made using commercial assembly planning and optimization software, like Valor [START_REF] Trilogy | Commercial assembly engineering software package[END_REF] or [START_REF] Siemens | Commercial assembly engineering software package[END_REF]. This kind of proprietary software use heuristic algorithms to solve the general assembly problem. But, there is no information available about how these software operate and which heuristic algorithms are used in them. It is only known that the optimizer software use separate programs tuned for each particular placement machine type. The running times for the optimization procedures are typically from minutes until several hours. In the general case the optimization includes both pick-andplace sequencing and component-to-feeder assignment which make the PCB assembly optimization still more complex to solve, see [START_REF] Leipälä | Optimization of the movements of a component placement machine[END_REF].

When considering possible methods for estimating operation times of manufacturing machines one can trust on statistical techniques. Common non-linear regression methods are well known in classical statistics literature, see for instance [START_REF] Hald | Statistical theory with engineering applications[END_REF]. In recent years, neural networks have gained much popularity due to their flexibility, versality and simplicity in usage. There are many different types of neural networks supporting non-linear regression e.g. MLP (multi-layer perceptron), RBF (radial basis function) and SVM (support vector machine), for a comprehensive study of this topic, see for instance [START_REF] Bishop | Neural Networks for Pattern Recognition[END_REF]. Each one of these methods has its own strengths and weaknesses. It turns out that conventional learning methods are prone to overfitting when the number of hidden units of the network is large in relation to the number of training cases. To avoid this and complicated training and testing, one can use Bayesian regularisation in MLP to achieve efficient learning and good accuracy automatically, see for instance MacKay, D.J.C. (1992a and1992b) and [START_REF] Neal | Bayesian Methods for Machine Learning[END_REF]. The number of the neurons in the hidden layer determines the flexibility (the number of free parameters) of the model. [START_REF] Hornik | Multilayer Feedforward Networks are Universal Approximators[END_REF] have showed that "feedforward networks with as few as one hidden layer are capable of approximating any Borel measurable function from a finite dimensional space to another to any desired degree of accuracy, provided sufficiently many hidden units are available".

When working with neural networks, there are two different operation phases; constructing the optimal layout and parameters of the network on the basis of the training data, and using the constructed network to determine the output for a new input. The training phase may take a long time (many hours) but it is done off-line. More importantly the use of a given neural network is very fast which is a necessity in the context of line balancing where one has to determine the assembly time of numerous different subsets of the components on the machines of the production line. The neural network model has also the further advantage of being adaptive; one can use the measured true assembly times to quickly retrain the model parameters.

Operation principles of placement machines

While different placement machine types normally include parts with rather similar functionality, many details of the machine design cause that each type should be analyzed separately. Here we concentrate on two rather different, well known and widely used machine types and discuss their design on the level which is necessary for understanding their main time factors.

Turret type placement machine

Turret-type placement machines (Figure 1) have been used in particular for high volume production. In the following, the discussion refers to the high-speed Universal HSP 4795/96 placement machine. The main body of this older machine type includes a conveyor belt, a PCB holding table and a rotary turret placement head for holding the vacuum nozzles. The feeder unit moves horizontally in the x-direction on the plane of the PCB holding table and comprises the required component reels. The number of component reels in the feeder unit is limited and each reel occupies a certain number of storage places, called feeder slots. This number depends on the width of the particular feeder reel. The task of the feeder unit is to bring the proper reel to the pick-up position where the placement head picks up a component with the appropriate vacuum nozzle. After that, the turret rotates stepwise around the vertical axis until the nozzle is at the placing position (which is 180° from the pick-up position). The intermediate positions are for inspection, rotation etc. of the components. Because the rotary turret contains several placement heads, there are several components ready in the turret and the feeder movement is performed a given be moved by two independent step motors to the proper location for placing the component. Usually for heavy components the speed of the holding table must be limited so that the components will stay in their proper places during the table movements. These components are usually placed last so that the limited speed does not need to be used with the smaller components.

The turret rotates on steps, and each step lasts a fixed time during which the feeder unit and the holding table can move freely (i.e., their movements are simultaneous to the rotation step and do not increase the total time for the assembly task). Mathematically speaking, we have the maximum or Chebyshev metric (L 4 ), where the largest time of component rotation, x-transfer, y-transfer and feeder transfer gives the final time of a particular placement operation. The control program tries to minimize the assembly time and, consequently, arranges the feeder tapes and the pick-up and placement operations so that 1. the feeder and PCB movements are short enough to be carried out during the rotation steps, and 2. if this cannot be realized, the overall length of the movements is short and long jumps occur parallel in both the feeder unit and the holding table.

Obviously, we can get accurate time estimates of the placement operations either by performing and observing the operations in the device in question itself or by simulating the operations using an optimizer-simulator. High time costs often rule out both of these methods in everyday practise and that is why there is a need for rapid and accurate estimators.

Gantry type placement machines

Nowadays, turret machines have been frequently replaced with collect-and-place machines. These are equiped with (one or more) gantries which carry rotating heads that collect the components from feeders and place them to the PCB. Collect-andplace machines can usually manipulate also bigger components that are not accepted by turret machines. The two-gantry placement machines of Siplace S-Series (Fig. 2) are typical examples of this design. Each gantry carries a collect-and-place revolver head with either 12 or 6 nozzles according to the component heights (max. component height for 12-nozzle-head is 6 mm and for 6-nozzle-head 8,5 mm). Both head types can collect a number of components and then place them sequentially on the PCB which is motionless during the placement operation. The placement heads operate in parallel on the two gantries and have their own stationary feeder slots and

Optimizing simulators for assembly engineering

To get an idea about how much time is needed to simulate and optimize the assembly operations of a single placement machine, we take a closer look on two commercial assembly engineering software packages. In both software packages the optimization process is similar. The process is complex and time consuming in practice but simple in theory; initial information on the assembly equipment configuration (e.g., assembly lines, placement machines, heads, nozzle changers, etc.) functions as raw data. This information is given as an input to the optimizer software which optimizes the configuration of the production line (feeders, nozzles etc.) and generates the placement programs/rules for the machines.

Optimizing simulator for rotary turret machine

The "Line Engineering Package" is a part of the "Trilogy" assembly engineering software of Valor Corporation Ltd, see Valor [START_REF] Trilogy | Commercial assembly engineering software package[END_REF]. The software package contains machine specific optimizers, also for our case, the Universal HSP 4795/96 turret machine. The Line Engineering Package operates on three optimization levels. At level 1 the optimizer generates a quick-and-dirty control program for the placement machine. This option suits for testing prototypes when the amounts of PCBs are small or the result of this optimization level is used as an initial solution for higher level optimization. The running time is typically tens of seconds and it is limited to one minute. At level 2 the optimizer is suited for normal production situations when the PCB batches are of a reasonable size. The running times are typically few minutes and they are restricted to one hour. The optimizer determines four reasonable feeder setups and chooses among them the most promising one. Level 3 is for the optimization of mass production and it takes typically one working day. The optimizer performs a full-scale optimization of the feeder setup. The difference of the value in objective function in comparison to level 2 is typically one percent. The training dataset has been generated by [START_REF] Laakso | Estimating the Production Times in PCB Assembly[END_REF] using the level 2 optimization (includes reasonable component-to-feeder allocation) and the optimization time has been restricted to one hour per PCB job.

The line engineering package bases its operation on a machine simulator tuned for each particular machine type. When the solutions found by the package are carried out in a real placement machine, the observed time differs slightly from the simulator (usually less than two percent) but the difference is stochastic with a zero mean. Nevertheless, the simulator assumes that the machine operates in ideal conditions (e.g. all components are accepted by the optical quality check and component reels in the feeder are in the sequence they are needed) and there are no production interrupts caused by some real world phenomena.

Optimizing simulator for gantry type machines

To simulate and optimize the assembly times for Siemens gantry type Siplace placement machines we used real data from PCB factories and Siemens' own, inhouse optimization software, see [START_REF] Siemens | Commercial assembly engineering software package[END_REF]. Siplace Pro assembly line control and optimization software converts basic CAD data into placement programs on Siemens Siplace placement machines and lines. Siplace Pro includes an optimizer as a separate program having another data model to optimize the placement programs at the machine and line level. There are various optimization scenarios: optimize setups for a single product, multiple products in a single family setup, or multiple products in multiple family setups. There are, similar to Valor's Trilogy package, also 3 different optimization levels. The running time for the optimization procedure can be limited by giving the number of the optimization steps or by giving a maximum running time. The training datasets for gantry type machines has been generated by Siplace Pro using the level 3 optimization (which includes also feeder allocation optimization) and the optimization time has been limited to ten minutes per PCB job.

Non-linear regression estimation of assembly times

Even though the PCB assembly time can be roughly estimated by using linear models, it is actually a non-linear function of the input parameters. As can be seen in Figure 8 (left panel), the input-output mapping seems to be rather smooth and continuous. However, it is also obvious that a linear estimator of the processing time is not sufficient. Therefore we now turn our attention into non-linear regression models.

Artificial neuron (perceptron) networks with linear activation function are known to be closely related to linear regression models. The need to be able to handle more complex dependencies has resulted in development of the perceptron models that include more layers and apply a non-linear function to make the activation of the perceptrons non-linear, scaled and differentiable.

The perceptron computes a single output from multiple real valued inputs by forming a linear combination according to its input weights and then putting the output through some nonlinear activation function, see Figure 3. The operation of the perceptron can be mathematically described as y = Ν(∑ w i x i + b) = Ν(w T x + b), where w denotes the weight vector, x is the vector of inputs, b is the bias and Ν the activation function. In this study, we evaluate the feasibility of a standard neural network based approach to the estimation of assembly times by using an MLP network in the regression. Our multilayer perceptron network (Figure 4) consist of three layers, there is one input layer (k input variables = k neurons), one hidden layer having n neurons and one output layer (1 output variable = 1 neuron). The number of the neurons in the hidden layer determines the flexibility (the number of free parameters) in the model. The three layer model is selected because [START_REF] Hornik | Multilayer Feedforward Networks are Universal Approximators[END_REF] has showed that feedforward networks with as few as one hidden layer are capable of approximating any Borel measurable function from a finite dimensional space to another to any desired degree of accuracy, provided sufficiently many hidden units are available". We calculate later on in Section 6.2.2 how many hidden units are sufficiently.

- --------> INSERT FIGURE 4 HERE <-------------Figure 4. MLP network having one hidden layer.

Input values of the input layer are fed as such into the layer upstream, into the hidden layer. Once the neurons for the hidden layer are computed, their activations are fed upstream to the output neuron. In a fully connected multilayer feedforward network, each neuron in one layer is connected by a weight to every neuron in the layer upstream. A bias is also associated with each of these weighted sums. The activation functions of the hidden layer are sigmoids whereas the activation function of the output neuron is linear.

When the network structure is fixed, there still exists a large variety of different algorithms and heuristics that can be used in the training of the neural network. The problems that are faced in the learning, i.e. estimation of the network parameters, include overfitting to the training data (the network learns the characteristics of the training data but is not able to generalize), slow learning rate and getting stuck in local minima of the cost function that the learning algorithm is minimizing (sum of squared errors).

One known way to avoid the overfitting problem is to use Bayesian framework having automated regularisation in practical MLP network calculations. The Bayesian framework presents uncertainties in the values of the parameters as probability density functions. Before presenting the data, the parameters are described by a prior probability density, which is typically quite broad to reflect the fact that we have little idea of what values the parameters should take. Once the data are presented, we can compute the corresponding posterior probability density. Since some values of the parameters are more consistent with the data to be presented than others, we find that the posterior distribution is narrower than the prior distribution. This phenomenon is known as Bayesian learning and is illustrated in Figure 5. --------> INSERT FIGURE 5 HERE <-------------Figure 5. A prior distribution of weights p(w) and the posterior distribution p(wD). The most probable weight vector w MP corresponds to the maximum of the posterior distribution.

In our experiments, we used the Bayesian regularisation of the network parameters, which is roughly equivalent to a technique called weight decay. The regularisation aims at keeping the weights small, and therefore prevents overfitting. As the training algorithm we used a scaled conjugate gradient method which is a faster substitute for the classical gradient descent algorithm. The classical gradient descent algorithm adjusts the weights in the steepest descent direction (negative of the gradient), the direction in which the error function is decreasing most rapidly. It turns out that, although the function decreases most rapidly along the negative of the gradient, this does not necessarily produce the fastest convergence. In the conjugate gradient algorithms a search is performed along conjugate directions, which produces generally faster convergence than steepest descent directions.

The number of neurons in the hidden layer should be selected so that the model has good generalization properties, that is, its prediction error for testing data (not used in training of the model) is as low as possible. We estimate the suitable number of neurons in Section 6.2.2. Because our data sets are quite small, we were forced to use the leave-one-out (LOO) validation of the generalization error in our experiments.

Experiments

Because the two machine types of Section 3 are suited for different kinds of PCB assembly tasks, we had tree different sets of training data for placement time estimations. For the turret type machine the data set was taken from an earlier study [START_REF] Laakso | Estimating the Production Times in PCB Assembly[END_REF] and for the gantry type machine two sets of real production data were obtained from two different PCB assembly factories.

Training data

The data we use in the context of rotary turret machine were originally generated by [START_REF] Laakso | Estimating the Production Times in PCB Assembly[END_REF] for studying multivariate linear regression estimators. The data were formed by creating artificial PCBs for placing n components having c different component types on a square-like PCB with a side length of l millimetres. The locations of the components on the PCB were chosen by sampling the x-and ycoordinates independently and uniformly from the interval [0, l]. The parameter values were n {50,100, 400, 700, 1000}; c {1, 30, 100}; l {0, 50, 200, 500}. The generated assembly tasks were then simulated using the level 2 optimizer of the "Line Engineering Package" for each instance of the 70 PCB jobs. The simulation results were considered as the "true" placement times, which we use to train MLPs to estimate the times.

The two training data sets for gantry type machine were real proprietary production data, obtained from two different PCB assembly factories. The first dataset (dataset 1) contained 63 different PCB jobs that were rather homogeneous since there were a lot of different versions of the same product. The second dataset (dataset 2) contained only 10 jobs which were from a rather different type of electronic product, in comparison to the jobs in the first dataset. Dataset 2 was not used alone since it is too small for the MLP to be used as training material. Each of these jobs was then optimized using the level 3 optimization of the Siplace Pro software to get the objective function ("true" placing times). The benefit of using this kind of real world assembly data is that the MLP can be trained using realistic inputs with corresponding simulated placing times.

Estimators for rotary turret type machine

Both linear regression models (univariate and multivariate) and the training data were taken directly from the earlier study made by [START_REF] Laakso | Estimating the Production Times in PCB Assembly[END_REF] to be compared with our non-linear MLP estimator.

Estimators based on linear regression

Two linear regression estimators for the HSP 4795/96 turret type assembly machine were presented in the study by [START_REF] Laakso | Estimating the Production Times in PCB Assembly[END_REF]. Both of these were based on Valor Line Engineering Package simulations. The univariate linear regression model for the assembly time (t) is of the form t = 0.117 ⋅ n + 6.991, (1) where n is number of components to be assembled. The model explains 85.1 percent of the time variation. A multivariate linear regression model is of the form t = 0.115 ⋅ n + 0.163 ⋅ c + 0.000536 ⋅ l -6.467, (2) where c is number of different component types and l the side length of the PCB. This model explains 92.2 percent of the time variation. Graphical results of this multivariate linear regression (green line) compared with simulated "true values" (blue line) and with absolute error values (red line) are shown in Figure 6.

- --------> INSERT FIGURE 6 HERE <-------------Figure 6. Multivariate linear regressions on assembly times for rotary turret machine (HSP 4795) by formula (2).

Estimator based on non-linear regression

In our experiments with non-linear estimators, we first tested how many hidden neurons our MLP model should use for non-linear regression of the PCB assembly times. The input features of the MLP were the same as in the multivariate linear regression in [START_REF] Laakso | Estimating the Production Times in PCB Assembly[END_REF]. In the training of the MLPs we used the Netlab toolbox for Matlab functions software created by [START_REF] Nabley | Netlab toolbox of Matlab functions[END_REF].

---------> INSERT FIGURE 7 HERE <------------- Figure 7.
The effect of the network complexity (number of hidden neurons) on the error magnitude of the estimation. Sum of the squared errors (SSE) is given for 16 network models. The solid line (with dots) corresponds to the mean of the predictive distribution. Dashed lines represent the standard deviation of the predictive distribution around the mean.

We estimated a suitable network complexity (number of neurons in the hidden layer) by using the leave-one-out (LOO) validation method. As the learning of the network is stochastic, we repeated the validation 10 times for each network structure. The results of the model complexity estimation measured in terms of the sum of the squared errors (SSE) are summarized in Figure 7. The model which gives the best results has complexity of 12 hidden neurons. The errors here are the differences between the "true" assembly times and estimated assembly times measured in seconds. To get more information about the behaviour of the 12 hidden neurons model we calculated also some other statistical characteristics of the models MSE (Mean squared error), MAE (Mean absolute error), RSE (Relative squared error) and RAE (Relative absolute error) for this model, see Table 1.

T ab le 1 . Co mp ar i so n o f st at is ti cal c har act er i s tic s fo r mu l t i var ia te li ne ar r eg r e ss io n ( 2 ) a nd t h e M LP mo d el fo r r o ta r y t ur r et ma c hi n e ( seco nd s) .

---------> INSERT TABLE 1 HERE <-------------Assembly times of the input data, the MLP model with 12 hidden neurons and the multivariate linear model are shown in Figure 8 for different number of component types. As it can be observed, the dependency between the input variables and the output variable is non-linear. The amount of non-linearity increases as the number of different component types in the assembly grows; if there is only one component type, the assembly time is in practice linear with respect to the two other input variables.

- --------> INSERT FIGURE 8 HERE <------------- -

--------> INSERT FIGURE 9 HERE <------------- Figure 9.
Results of the MLP model having 12 hidden neurons for rotary turret machine compared to the true values.

In Figure 9 the results of MLP having 12 hidden neurons are presented on the same visual manner as earlier for the multivariate linear model (Figure 6). The model having 12 neurons was found in Figure 7 to have the smallest generalization error. The graphs of the true values (blue) and the MLP approximations (green) are so closely one on the other that they can not be seen separately.

These results show that we can make apparently good assembly time estimations using MLP networks. However, one should note here that our data set from [START_REF] Laakso | Estimating the Production Times in PCB Assembly[END_REF] was artificially created for a particular machine type as described in Section 6.1.

Estimators for gantry type placement machines

Preliminary test showed that the use of a combination of datasets 1 and 2 as a training data resulted in less accurate results than the use of dataset 1. That is why we decided to do two sets of tests, one for dataset 1 and the other for the combination of dataset 1 and 2. Dataset 2 was not used alone since it is all too small for the MLP to be used as training material. We tested the multivariate linear regression and MLP techniques with more parameters and multiple input parameter combinations to find the best combination.

T ab le 2 . I np ut p ar a met e r s u sed fo r t he ga nt r y t y p e mac h i ne cal c ul at io ns.

- --------> INSERT TABLE 2 HERE <-------------

The accuracy of the estimation was measured with three error values; MSE, MAE and MAX. MAX is the maximum value of the estimation error which is the absolute difference in seconds between the correct value and the estimated value.

Estimators based on linear regression

Some key results from the multivariate linear regression with different input parameter combinations from Table 2 are shown in Table 3. It is natural that the accuracy of the linear estimators decreases with the number of different job types since presumably the variance of the jobs will also increases due to the use of different component shapes. That is especially the case when dataset 2 is used in conjunction with dataset 1. However, the same effect can be later seen also with the non-linear MLP.

T ab le 3 . M SE, M AE, M AX a nd t h e sq u ar e o f P e ar so n p r o d u ct mo me n t c o r r ela tio n ( r 2 ) o f t he e s ti ma ted p la ci n g t i me fo r d i f f er e n t i np u t p ar a met er co mb i na tio n s in t he mu l ti v ar i ate li n ea r r e gr e s sio n. S y mb o l s a to f r e fer to t he p ar a me t er s o f T ab l e 2 .

---------> INSERT TABLE 3 HERE <-------------As it can be seen from Table 3, the parameter combination a, b, d, e, f seems to be the best (at least when looking at the MSE value) for both datasets. The best parameter combinations include the inputs e or f, so it is rather obvious that the variance in the component sizes affects the placement speed. That is probably because the usage of the nozzles may not be optimal when there are a lot of different sized components that require nozzles that can be used to place only a small subset of all components. The differences between some of the input parameter combinations are marginal and it is difficult to say which combination would give the smallest errors with a greater amount of data. The dependency between the number of components and placing time is shown in Figure 10. It is clear that the total number of components heavily dominates in the regression of the total placing time but improvements can be achieved by selecting additional input parameters.

- --------> INSERT FIGURE 10 HERE <-------------Figure 10. Dependency between the number of components and the placement time (asterisk = dataset 1 and square = dataset 2).

6.3.2 Estimator based on non-linear regression MLP (one hidden layer with 12 neurons) was trained and tested in the same way as in Section 6.2.2. Table 4 shows some of the results obtained with the MLP. T ab l e 4 . M SE, M AE a nd M AX o f t he es ti ma ted p la ci n g ti me fo r d i f fer e n t i np ut p ar a me ter co mb i n at io n s i n t he M LP . S y mb o ls a to f r e f er to t he f eat ur e s o f T ab le 2 .

---------> INSERT TABLE 4 HERE <-------------

The results for MLP are interesting since for most input combinations the accuracy decreases a lot when the dataset 2 is included in the training set. Also the best input combination is different when the dataset 2 is included. It must be noted that for dataset 1+2 the linear regression performs generally better than the MLP. The reason for the bad performance of the MLP may derive from the sparseness of the training set. Also the data may be so intense in certain areas that the MLP is distorted by those "clusters".

The Bayesian regularisation of the network parameters in the combination with the scaled conjugate gradient turned out to be somewhat problematic in our case due to the high number of inputs in comparison to the size of the training data. For high number of inputs the number of the regularisation cycles was reduced from three to two which decreases the accuracy of the estimation.

The two data sets included a number of cases where there was a basic PCB type and its minor variants. We therefore clustered the data in order to remove the possible effect of intense data spots in MLP training. When doing this, the number of clusters should be rather high so that only the most intense data spots would be identified as a single cluster. Clustering was done so that each identified cluster was replaced by a single data point which was the arithmetic mean of the points belonging to the clusters.

- --------> INSERT FIGURE 11 HERE <-------------Figure 11. The effect of clustering to the error values in MLP with datasets 1+2 for gantry type machine.

We tested the MLP with different numbers of clusters using the input combination a, b, e, f (this combination had the best average accuracy for both input datasets), see Figure 11 for a summary of results. It appears that clustering gives no actual benefits since there is no remarkable improvement on the estimation accuracy while the cluster count is high. When 73 clusters are used the case is the same as without clustering because there are total 73 objects in dataset 1+2. Significant decrease of the amount of clusters will naturally decrease the accuracy since objects that have a rather large distance from each other are placed into the same cluster.

Conclusions

Production control in fast-paced electronic assembly industry requires proper solutions for production scheduling and line balancing to ensure that the promised dead lines are met and expensive machines are used efficiently. To solve these problems in a short time a fast single machine placement time estimation scheme is needed. Not only that the estimation must be fast enough it must also be accurate since better estimation leads to better line balancing. Methods presented in this paper are directed for factory scheduling where heavy machine level optimizations are all too slow or need information that is not yet available.

We introduced the use of multi-layer perceptron neural network for the placement time estimation in PCB component assembly. The method was compared to multivariate non-linear regression for two rather different machine types. For the rotary turret the new method was clearly superior to the multivariate linear regression (Table 1). This mirrors the non-linearity of the processing times as demonstrated in Figure 8. For gantry type machines the results show the importance of a proper training set. When having a relatively homogeneous training set (dataset 1) the MLP neural network works again efficiently (Table 3 compared to Table 4). With very heterogeneous data (datasets 1+2) multivariate non-linear regression gives somewhat smaller MSE. On the other hand, this kind of data is an extreme case which is not very suited for modelling any how.

When considering the practical usability of MLP neural networks it is observed that the application of the technique includes four steps; collecting a training data set, selecting a proper network layout and optimizing the network weights, using the network in an application (like line balancing or cost accounting), and adapting the network to new observations. Out of these the first step, collecting a proper training set, has been most challenging to us. While there are plenty of machine data available from producers, creating a well-founded set of observations (PCB assembly definition along with observed manufacturing time) has to be done with care. Of course, this note holds also for the competitors of the neural network approach. The second step, selecting the layout and optimizing network weights, is standard one and when the layout has been chosen the same general structure can (hopefully) be used longer time. The step takes several hours but this is again true for other techniques, too. The third step, using the network, is fast as the application of the ready network can be written as a function. This even makes it possible to use the method in conjunction of exact mathematical optimization software; research on exact line balancing is currently going on in our group. A further advantage of the neural network approach is the possibility to retrain (step four) the network weights in an easy manner when more production data have been collected. This naturally presupposes the preservation of the general layout of the network.

As a topic of further studies one should consider the modelling of so-called fast and slow components in the same neural network. For turret machines one may have to slow down the machine speed in order to avoid unwanted displacements of some heavy components due to high accelerations of the PCB holding board. Modelling of two operation speeds with the same network seems to be possible but needs some care. An other open question is the consideration of the capacities of the placement heads (to hold different nozzle types) in relation to the number and kind of component types of PCBs. In the current version the head capacity is fixed but one might let it be a parameter in order to generalise the model.

Figure 1 .

 1 Figure 1. Main components of a turret type placement machine (top view).

  Figure 2. Main components of the Siplace S-Series two-gantry placement machine (1: gantries, 2: feeders, 3: revolver heads, 4: nozzle, 5: PCB).

Figure 3 .

 3 Figure 3. Artificial neuron (perceptron).

Figure 8 .

 8 Figure 8. The input-output mappings of the test data of the MLP network for rotary turret machine. The columns of panels show the dependency of the assembly times on different combinations of system parameters (number of components, number of component types and PCB side length) for input data, MLP network having 12 hidden neurons and the multivariate linear model (2).

  Figures and tables for "Estimating PCB assembly times using neural networks"

  Figure_3.tif

  Figure_6.tif

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