Fabio Gabrieli

Riccardo Artoni

Andrea Santomaso

Simonetta Cola

Discrete particle simulations and experiments on the collapse of wet granular columns

Keywords: numbers: 45, 70, Ht 47, 57, Gc 83, 80, Fg

Small quantities of liquid in a granular material control the flow dynamics as well as the triggering and jamming phases. In order to study this problem some experimental collapse tests conducted in a rectangular box were reproduced with a 1:1 scale numerical model using the Discrete Element Method (DEM). In simulations the effect of the capillary bridges has been investigated implementing a mid-range attractive force between particles based on the minimum energy approach. Also a bonding-debonding mechanism was incorporated in the algorithm and the volume of each sessile drop on the particle surface was considered during its motion. The influence of some variables was investigated with respect to the final slope profiles and the runout lengths: the initial liquid content, the particle size, the solid density, the liquid surface tension and the liquid-solid contact angle. Also the crucial effect of the confinement walls on the collapse phenomenon was investigated: wet particles adhere to the lateral walls providing a higher flow resistance in comparison to the same material in dry conditions. It was observed that particles with largest path-lengths are localized near the movable wall at a middle-height of the initial column sample. Other particles at the surface moves in a rigid way especially if they were wet and with a low solid density. The "fidelity" of each particle with respect to the nearest neighbours was evaluated allowing to recognize the emergence of clusters of particles and rigid parts, to extract the failure surface and to localize where debonding mechanisms concentrate in the wet case.

I. INTRODUCTION

At least once in our lives we have played with sandcastles: experience has told us that in order to construct nearly vertical walls we needed to add a small amount of water to the material. Now we know that the stability of the sandcastles is due to the attractive forces between particles that arise only in partially saturated conditions.

In the last years this topic was studied by some authors [1][2][3] that realized the many existing links between the simple sandcastle and more complex civil engineering and geological problems like the assessment of the slope stability. A recent study 4 showed that the stability of cylindric sandcastles depends on the volumetric water content (in experiments with sand the authors found an optimum value of 1%) and on the aspect ratio (the maximum height is function of the base radius to the power of 2/3). They also reported the possibility to build underwater sandcastles using hydrophobic sand: in this case the "cohesive" properties of the partially saturated material were conserved by the presence of air menisci in place of water menisci.

Conventionally the saturation condition can be classified into five regimes 5 corresponding to five different physical conditions of the granular material. The dry state and the a) Electronic mail: fabio.gabrieli@unipd.it completely saturated conditions are the most known. The others are the pendular regime (when liquid is shared by pairs of grains), the funicular state (when liquid menisci coalesce grouping three or more particles) and the capillary state (when bubbles and menisci are present at the same time).

The completely saturated and the perfectly dry conditions (two-phase systems) share many similarities with respect to the mechanics of the granular material in static conditions. On the other hand, at intermediate conditions, the contemporary presence of liquid and air in the granular mass (three-phase systems) completely redefines the mechanical relationships between individual particles influencing the material behaviour at meso-and macro-scale.

At the contact scale, in static conditions, liquid-air interfaces are responsible of attractive (rarely repulsive, for convex menisci) forces between particles while lubrication, viscous and damping forces are reported in dynamic conditions 6,7 . At meso-and macro-scale some important secondary effects can be generally observed: a spontaneous agglomeration (granulation), an increase of the bulk density, a higher maximum packing density and an increase of the shear strength of the material. Many of these phenomena were studied at an "engineering scale", and usually the micromechanical interpretation is missing or unclear.

Also, the dynamic regime and the transitional phases (triggering and jamming) represent other important open issues that need to be addressed especially in the case of partially saturated conditions. The role of the liquid phase, the effect of the particle size and shape are only two of the several aspects that have important consequences in the landslide triggering in soil mechanics as well as in jamming/triggering problems in industrial processes such as handling, storage and transport. The results of the studies in this field can be the keystone to connect the classical solid mechanics (for quasi-static little deformations) with granular flow theories (for fully developed flows).

Collapse tests were considered as good reference tests in evaluating dynamic behaviour of liquids and granular materials in transient conditions 8 since they incorporate the triggering phase, the "fully-developed" flow and the final jammed state. Several authors studied the problem of the column collapse using dry granular materials [9][START_REF] Lajeunesse | Vilotte[END_REF][11][12][13][14] , investigating the effect of the particle size and shape, the material type, the initial column geometry, and observing the shape of the final slope, the segregation, the kinematic of the front toe and other variables. Recently, other works focused on the collapse in fluids, which is the first step for the comprehension of the fully saturated granular flows, highlighting the effect of other variables which emerge in this case like the dilatancy behaviour linked with the initial volume fraction 15 , and the opposing effect of viscous and lubrication forces 16 . Besides the laboratory tests many authors used particle numerical methods like the Discrete Element Method to model these experiments and to derive other useful contact-scale information [17][18][19][20][START_REF] Girolami | Granular Matter[END_REF] .

The aim of these experiments and of the simulations was to collect information for the setting up of a continuum framework, in order to solve granular flow problems at an "engineering-scale". The original strategy was to use the depth-averaged shallow-water equations and to treat the material as a sort of viscous-frictional fluid 18,[START_REF] Savage | [END_REF][23][24] . However, as remarked by some authors 18,24,25 , some difficulties arise using these approximations and other continuum approaches should be preferred to recognize the specificity of the spreading of the granular materials 26 .

Our work wants to validate a particle-based model for the collapse tests of wet beads in pendular regime and to recognize the main micro and macromechanical variables that affect this problem.

We started from the study of the capillary force between two spheres, to provide a particlescale model based on the Discrete Element Method taking into account liquid bridge failure, formation and redistribution. This model was used to reproduce and validate some collapse experiments at 1:1 scale with wet monosized glass spheres. After that some variables were modified in experiments and simulations: we focused on the use of different particle diameters, liquid amounts, surface tensions and boundary conditions. Other "virtual conditions", difficult to address in laboratory experiments, like the solid density of the particles and the liquid-air contact angle, were only simulated using the numerical model. The profiles provided by simulations were compared with the experiments from a qualitative and quan-

= Fcap/(ρgV sph ), V * = V /V sph )
titative point of view, highlighting the critical aspects of the model assumptions; then, the particle dynamics was examined extracting useful information about the mechanisms of collapse for dry and wet initial conditions at different solid density.

II. THE CAPILLARY BRIDGE BETWEEN TWO EQUAL SPHERES

Let's consider a small liquid volume between two bodies. When the liquid surface tension and the solid-fluid adhesion prevents the fluid to flow away under the action of external forces (e.g. gravity), a pendular bridge is formed. At a microscale level, the meniscus is constrained to a finite set of possible shapes since its stability is controlled by an equilibrium between pressure in liquid phase and air pressure depending on the solid surface geometry and the physical interface properties.

There are at least two theories concerning the calculation of these shapes: the first is the classic Young-Laplace problem 27 that consists in solving a boundary value problem containing the expression of surface geometry of the menisci and the boundary conditions. The second concerns the solution of a minimum total energy problem that considers the energy of the system as the energy of the liquid volume added to the energy at each interface.

It can be proved that when two solid particles share a liquid phase the equilibrium of the system is not generally satisfied and attractive or repulsive forces emerge [START_REF] Israelachvili | Intermolecular and surface forces / Jacob N. Israelachvili[END_REF][START_REF] Willett | Granulation[END_REF] . The interparticle force in a certain direction can be computed in two ways: as derivative of the total energy in that direction or decomposing all the forces in a section normal to that direction. Each of the two ways generally leads to a numerical solution 27,[START_REF] Rabinovich | [END_REF]31 . Recently Lambert et al. 31 proved that the two approaches are equivalent for some simple geometrical shapes.

The most common solution of the menisci shape regards the liquid bridge between two equal spheres. Since generally no closed-form solution exists for this problem many other correlations were derived in the past. They can be classified in numerical and analytical/geometrical approximations 27,31,32 , solutions derived from fitting of the latter 33,34 and solution derived from fitting of laboratory measures 35 . In Figure 1 a comparison of some selected normalized solutions is reported for the sphere-sphere case and for different capillary volumes.

Regardless of the solution adopted, qualitatively all the models describe the capillary force as a function of the capillary volume V , the gap s, the contact angle θ and the filling angle β (see Fig. 2). Pressure information (of air and/or liquid) should also be included especially if gas domain is not connected (it is not the case of the pendular state). In the following we will adopt the Lambert solution 31,36 that is simple, close to the laboratory data, suitable for integration into DEM and theoretically based at the same time.

All the models agree that the capillary attraction force has a maximum value (2πRγ cos θ) for gap, s = 0 (particles in contact) which is constant and independent of the liquid volume. This evidence was experimentally confirmed for a wide range of capillary volumes excepted for very small liquid volume where thermodynamic considerations and details about the surface roughness play a major role 37,38 . The force decreases with the distance up to a certain gap s f at which the bridge breaks up. Lian et al. 27 observed that this critical distance is function of the liquid volume V and the contact angle θ:

s f = 1 + θ 2 V 1 3 (1)
On the other hand, the critical distance at which a given liquid volume can form a pendular bridge is [START_REF] Rabinovich | [END_REF] :

s c = 12V π 1 + 1 1-cos θ (2)

III. PENDULAR STATE DEM ALGORITHM

The algorithm used in reproducing the wet granular material behavior is the classical Discrete Element Method (DEM) 39 that integrates the motion of each particle starting from its initial position, imposing the external forces (like gravity) and computing contact forces derived from the constitutive model at the contacts. For dry conditions, the classical contact model consists in a linear spring with a viscous dashpot in normal direction and a linear spring plus a slider and a dashpot in tangential direction. DEM allows in a simple way the addition of variable custom forces thanks to the explicit nature of the algorithm which is based on the separation between the integration of the Newton's laws of motion and the computation of the forces at the contact. In our case, a liquid bridge model was implemented to study the material in the pendular regime.

The liquid bridge model is constituted of three parts: the first part manages the initial distribution of the liquid among contacts and particles and it is called at the initial steps.

The second part couples the classical DEM contact mechanics with the capillary force model. It treats the capillary attraction forces as mid-range forces in addition to the dry contact force mechanics that instead act only for overlapping particle. This part was treated in detail in Gabrieli et al. 36 considering specifically a two-particle system with different kinematic constraints.

The third part checks the formation of new capillary bridges from single or multiple drops and the failure of liquid bridges that reach the rupture distance (see Appendix B for the details). The total liquid volume was maintained constant during the simulations.

A. On the initial distribution of liquid

The liquid distribution in a granular material depends on many factors: the liquid amount, the initial solid packing density, the liquid properties and the "history" of the sample. As remarked by Scholtes et al. 40 , the initial capillary bridge distribution has an important effect on the mechanical response of the material, therefore an initial distribution that mimics the initial real sample is desirable.

For this purpose, some micro-tomographic tests were performed on small cylindrical samples of glass spheres mixed with water. The observation of the reconstructed slices allowed to highlight that the initial distribution of liquid volume in the mass was homogeneous and well distributed on about the 96% of contacts [START_REF] Artoni | Discrete Element Modelling of Particulate Media[END_REF] . No triple menisci could be recognized; this was verified for liquid contents typical of the pendular regime (w < 4%).

In order to obtain a liquid content in the numerical sample consistent with the experimental observations, the following strategy was adopted: a) the couples of particles sustaining a liquid volume (called "potential contacts") were identified as those with a distance lower than s f (from Eq. 1 this value in turn depends on the capillary volume); b) the total liquid volume was subdivided randomly among 96% of the potential contacts starting from a trial value of the capillary bridge volume and iteratively verifying the total amount of fluid, checking the critical distances and the homogeneity of the distribution. At the end of this procedure each capillary bond had the same liquid volume. Using this method, it was possible to find an initial stable particle configuration and the assembly was allowed to relax without any artificial localized tension.

B. The capillary force model

The capillary force model represents the heart of the pendular model. Let us consider two equal particles sharing a capillary bridge. Using the interface energy theory 28 it is possible to write down the analytical expression of the total energy of a bridge. Deriving this energy with respect to the gap between the particles, with some geometrical approximations, the attraction force exerted on two particles can be obtained. This finally depends on the capillary volume V and the gap s (see Appendix A for the details):

F cap (s, V ) = - 2πRγ cos θ 1 + [s/2b(s, V )] (3a) b(s, V ) = (s/2) -1 + 1 + 2V /(πRs 2 ) ( 3b
)
where γ is the liquid surface tension, θ the contact angle and b the wet spherical segment on sphere surface. Eq. 3a is valid for s < s f ; otherwise the bridge breaks up and no attractive force exists. With this assumption the capillary force model describes a midrange interaction, with a cutoff at s = s f . When particles are in contact, the attractive pendular force is added to the classical soft-sphere repulsive contact force; if not, only the capillary force (in addition to gravity) is considered for integrating the equation of motion. Because of their mid-range nature, the implementation of capillary forces involves an increase in the computational cost of the algorithm due to the enlargement of the Verlet cell of every particle. a Measured with the ring method.

b Measured taking telecentric photos of a sessile drop on a flat surface cleaned with ethanol.

IV. COLUMN COLLAPSE EXPERIMENTS

Collapse experiments are quick laboratory tests, quite simple to perform and that require only the use of a camera. From the time-lapse images of the collapse, using digital image analysis techniques, "in-plane" quantities like velocity field or profiles of the slope can be easily measured.

Our experimental program consisted of several collapse tests performed within the same initial cuboid geometry (size 7 x 8 x 5 cm) in a transparent glass box (see Figure 3). The initial granular column was created by pouring the wet materials with a funnel in the fixed box. To start the test one of the glass lateral walls was lifted by a weight/pulley system. Each test was captured using a semi-professional high-resolution camera (Casio EX-F1) using a frame rate of 30 fps. The analysis of the image sequence enabled us to obtain the time evolution of the granular mass and to measure the profiles at each time step.

The sequence of images were analyzed to extract the surface slope of the material from the beginning of the collapse up to the final rest configuration of the pile. The mean final profile of a group of tests was obtained by averaging the z coordinates of the surfaces along the vertical axis. From the surface profiles, the slopes at different locations and the runout lengths L f were extracted (see Fig. 3).

Mono-disperse glass spheres with different size (2, 3 and 5 mm in diameter D) were used in various wetting conditions (liquid contents w = 0, 0.5, 1, 2 and 4% of the weight of the dry materials). Two wetting liquids were used: distilled water and a mixture of water and a fluorinated surfactant (1-perfluorotooctyl-3-propan-2-ol with a concentration of 0.5 g/l) that allowed to decrease the surface tension. The physical parameters for the liquid phase (liquid surface tension γ and contact angle θ), were measured with dedicated laboratory tests and were reported in Table I.

Figure 4 shows two typical final configurations of the pile of material after the collapse in the dry and the wet cases. It is evident that the main effect of the capillary forces is to give a shorter run out and as a consequence a larger slope of the final profile. However in order to better understand the phenomenon a more detailed analysis is required. For the sake of shortness only some selected results will be presented here while further details can be found in Artoni et al. [START_REF] Artoni | [END_REF] .

Figure 5 represents the mean normalized runout length obtained in all the experiments. Each group of tests with the same diameter was connected with a piece-wise line in order to facilitate the reading of the trends.

The first observation is that the collapse characteristics depended on both liquid content and particle size: the smallest particles (D = 2 mm) showed the strongest dependence on the moisture content while considering the coarser material (from D = 3 to 5 mm) the runout length decreased weakly with the moisture content (with exception of the case w = 0.5%, D = 5 mm).

Furthermore, it should be noted that for the finest material the variation of runout length with liquid content was greater passing from dry to w = 1% and it reached a quite constant value for higher moisture content values (w > 1%).

Lowering surface tension generally shifted the runout length values towards the dry results (with exception D = 2 and 3 mm, w = 0.5%). Even for the coarsest material the "wet" runout lengths are slightly higher than the dry values; this effect could be attributed to the lubrication of the contacts that could arise in these conditions.

As already noted in other works 43,44 , due to the rectangular geometry, an effect of lateral walls was found: the materials stopped with a slope profile and runout length that slightly depended on the ratio between the thickness of the box (W ) and particle size (D). This is evident considering the dry coarsest particles (W/D = 10) for which the normalized runout length was ∼ 13% lower than with the finest material (W/D = 25).

V. DEM SIMULATIONS

The experimental tests were reproduced with 1:1 scale simulations with the aim of validating the numerical model and investigating the effect of different physical variables (particle size, liquid content, surface tension) on the column collapse dynamics. Numerical simulations have been also extended to conditions not studied in the laboratory such as lower solid particle density, different contact angles, different initial geometries and different lateral wall removing mechanisms. The study was focused on two particle sizes, D=2 and D=5 mm.

The granular mass was simulated using the pendular state DEM algorithm described in section III. The number of particles was fixed in order to simulate the same weight of material adopted in the experiments. This corresponded to 2250 and 37500 particles for the 5 mm and 2 mm tests respectively.

In the simulations the wall was lifted up imposing an upward acceleration equal to the gravitational acceleration, as in the laboratory tests. However, since in literature many numerical collapse tests have been carried out applying a instantaneous horizontal displacement to remove the wall, or without mention about the technique, this aspect has been specifically addressed with dedicated simulations. In particular two different wall removal mechanisms have been simulated: the instantaneous removal of the wall and the removal by uplift in vertical direction with acceleration g. Only slight differences were observed for the two cases and therefore wall removal mechanisms will not be considered further.

The surface profiles were extracted from numerical simulations considering an average profile on a internal slice 2-cm thick located on the axis of symmetry of the test device.

The DEM model required the assignation of physical and contact parameters that were calibrated on the base of dedicated laboratory experiments. The solid density of the spheres was the glass density. For the contact parameters the elastic normal and tangential stiffness (k n and k t ) of the springs and the contact friction coefficient (µ = tan φ µ ) of the slider were calibrated by comparing results of triaxial tests performed with glass spheres at confinement stresses σ 3 equal to 25, 50, 100 and 200 kPa with true cubic triaxial test simulations 36 . The damping parameter was at first roughly estimated by simple restitution coefficient tests and then was finely tuned by comparing experimental and simulated final collapsed profiles of 2 mm spheres in dry conditions.

Table II lists the values of all the contact parameters used in the DEM simulations. The same parameters (stiffness, damping and friction) were applied in wet conditions and to the sphere-wall contact model.

A. Experimental vs numerical results

DEM simulations were first compared with laboratory tests focusing on the final slope profiles and run-out lengths in order to verify the behaviour of the model.

Figure 6 shows the qualitative behaviour of two simulations carried out in the dry and wet case, analogously to what was done in Figure 4 for the lab tests. A more direct comparison is presented in Figure 7 where some selected final profiles obtained for 2 mm spheres are shown. Finally, a summary of all the comparisons between the lab and simulated tests is reported in Figures 8a,b in terms of run-out length values for different liquid contents and liquid tensions.

Since the dry test with D = 2mm was used for calibration it is not surprising that a good matching was found in slope profiles (see Fig. 7a) and run-out values.

Also for the wet cases the shape of the slopes fairly agree with the experiments: the final repose profiles result similar both with distilled water and adding the surfactant (see Fig. 7b,c). For these two cases small discrepancies regard the top and the toe angles which appeared slightly lower than the real ones. This reflected, more generally, on the runout lengths of all the finer particles which were systematically overestimated for all the moisture contents, both with distilled water and with water+surfactant (see Fig. 8a,b). On the contrary, for the coarsest spheres, the run-out values were better estimated for the samples mixed with water (they were inside the error bars), while they resulted different in values and trend for the low surface tension case.

A closer examination of the mechanical behaviour of the particles at the toe can provide a partial answer to the overestimation of the runout values. The simulations revealed that at the toe there were particle rotations despite of the presence of many capillary bridges. This was supposed to be unrealistic since the liquid bridge between two moving surface should oppose a rotational resistance which is closely related to the hysteretic behaviour of the contact angle and the consequent energy dissipation. With respect to this aspect, a model for the rolling resistance of capillary bridges should be included in future analyses.

With regards to the discrepancies observed with the coarsest spheres in the wet case the lack of a lubrication model at the contacts may be the main reason. Lubrication force is a repulsive force that reach its maximum during collapse, lowering the frequency of frictional contacts. This is consistent with the fact that the differences regarded the coarsest spheres with the highest liquid contents and that the magnitude of the lubrication force increased with diameter and liquid volume. However, a clear and straightforward understanding of the underlying physics would require dedicated experiments and a more complex contact model which, at the time, cannot be included in this work.

Despite of these differences, we retain this model sufficient to provide important informations about the effect of the single variables of the model to the overall behaviour of the material. In order to investigate the effect of each single variable on the final result we chose as reference case the test with the following parameters: D = 2 mm, w = 1%, γ = 72.75 mN/m, θ = 15 o , ρ s = 2532 kg/m 3 . The first group of simulations has been performed varying independently the particle size (D = 5 mm), the liquid properties (γ = 17.2 mN/m; θ = 30 o ) and the solid density (ρ s = 1266 kg/m 3 ). The analysis of results was focused on the final profiles, that were reported in Figure 9 and compared with the reference simulation.

With the same water content but with a larger particle size (D = 5mm), the runout length resulted to be longer and the angles lower (see Fig. 9a) than in the case with finer material. The same trends were observed with a low liquid surface tension (γ=17.1 mN/m as in the laboratory tests) (see Fig. 9b).

Instead, increasing the contact angle θ in a reasonable range (from 15 to 30 o ) slightly reduced the angle at the top while the runout length and the shape at the toe remained quite unvaried (see Fig. 9c).

The scarce influence of the contact angle on the final profile may be linked to the fact that its contribution on the capillary force reduced to cos θ and that for this hydrophilic surface its value is low, implying that cos θ ≈ 1 for the different values. Furthermore, it should be noted that the contact angle value reflects on the debonding-bonding distance which decreases for larger values leading a less frequent breakage of the capillary bonds. The counterbalanced effect of these two aspects probably yields a low dependence of the profiles on the contact angle values.

On the other hand, decreasing the solid density of one half strongly modified the final profile which was then characterized by a rectangular step on the rear and a curved toe (see Fig. 9d). With regards to this test a further inspection of the dynamics of the phenomena will be reported in section V E.

C. Effect of the liquid content

In the simulations, the effect of the liquid content on the final profiles resulted similar to those observed in laboratory experiments. It can be observed considering the case with D= 2mm and distilled water by superimposing the simulated slope profiles at different moisture contents as in Figure 10: the increase of liquid content produced steeper final slope profiles. Dry and wet conditions appeared well distinct while for w > 0 the differences seemed slighter and concentrated at the crest of the slope. At the toe the differences can be appreciated considering the runout length plot (see Figures 8a,b): it was clear that it decreases for increasingly higher water content reaching a sort of asymptotic value. Such behaviour was observed, although in a more gentle way, also for the case with the low surface tension value (see Fig. 8b).

On the other hand, the simulations with coarse spheres (D = 5 mm) exhibited a smoother decreasing trend of the runout length with liquid content and quite constant. In this case the runout values differed from the experimental ones especially for the low liquid tension for which the experimental data indicated an opposite trend with larger runout for higher liquid content.

In literature, the presence of an asymptotic behaviour for increasingly higher water contents was also reported in the shear strength curves of unconfined radial and axial compression tests 33 and direct shear tests 45 conducted both with spheres and other granular materials. In all these tests the effect of water content seemed to saturates for w values

0 0.1 0.2 0.3 0.4 0.5 1.2 1.4 1.6 1.8 2 2.2 2.4 (Bo -1 w 2/3 ) 1/3 Normalized run-out length, (Lf-L0)/L0 2D-Exp 2W-Exp 2T-Exp 5D-Exp 5W-Exp 5T-Exp Fit-Exp 2D-Sim 2W-Sim 2T-Sim 5D-Sim 5W-Sim 5T-Sim Fit-Sim FIG.
11. Runout length vs the dimensionless number Bo -1 w 2/3 for the experimental (Exp) and numerical (Sim) tests. The number in the legend denotes the particle diameter, while the letter D refers to the dry case, W the case with distilled water, T the case with water and surfactant.

between 1 and 3% depending on the material used.

For a micromechanical interpretation of this behaviour we can consider a simple twoparticle system. In Figure 1 the curves for capillary force vs. interparticle distance for different capillary volumes are presented. For a fixed gap the volume reduction cannot be larger than a certain value corresponding to failure of the capillary bridge. On the other hand, the volume increments slightly increase the capillary force (analytically we have

F cap ∝ -e - √ 1 V
following Soulie et al. 33 or F cap ∝ -V -1 2 for our capillary force model 31 ) and the number of potential capillary contacts (with also low capillary forces due to the larger gaps) because s f ∝ 3 √ V . Furthermore, as already mentioned in section II, it can be noted that for bonds with s = 0 the capillary force is constant and independent of the liquid volume and failure is not theoretically possible.

Considering the macro-scale, especially if the material is dense, an increment of liquid content does not increase the number of capillary bonds in a sensible manner (and the new bridges have low capillary forces due to the large gaps) and it doesn't increase the capillary forces of those already present (which have s = 0 for the majority). As result, the average of the capillary forces of the packing, in comparison to the other forces, tends to a quite constant value for increasingly higher liquid content.

Again, according to the dimensional analysis on the experimental results [START_REF] Artoni | [END_REF] , we can try to evaluate the contribution of the moisture content on the runout using a function including the Bond number, defined here as the ratio between the body and the capillary force (Bo = ρgR 2 /γ) and the liquid content w.

In Figure 11 the normalized runout values of laboratory tests and simulations were summarized in terms of a combination of Bo and w. All the data are quite well collapsed and the least square linear fitting provide a good data interpolation that confirms the adopted dimensional model. However the comparison with the fitting of the experimental data reveals the tendency for simulations to provide larger values of runout for the small diameters, and lower values for the coarse diameter and low liquid surface tension. Moreover it can be noted that the simulations of wet coarsest spheres and low liquid surface tension correlated with the liquid content while the experimental data did not, as previously commented. (g) as in (a) but with liquid content w = 0% (dry);

D. Effect of the boundary walls

In order to appreciate the effect of side walls on the collapse it is useful to observe the shape of the slope surface from a virtual top aerial view. Figure 12 collects several simulations results for this purpose.

The first observation is that the shape of the profiles depended on the distance from the lateral walls: the angles of the slopes were generally gentler and with a larger runout length in the middle section than near the walls. Moreover, observing the contour lines, two different zones related with different collapse mechanisms could be recognized: in the upper part of the slope the surface of the final deposit had preferentially a positive (sometimes low) concavity in transversal direction; on the other hand, at the toe of the profile the surface presented a more pronounced negative concavity.

Passing from the reference case to the case with D = 5mm and then to the case with a low liquid surface tension (Fig. 12a,b,c) we can observe how the retaining effect of the boundary walls tends to disappear becoming analogue to the dry case (Fig. 12g). The differences are magnified by the absence of particles glued to the lateral wall in the upper portion of the column which reveals the reduction in magnitude of the capillary force. As previously discussed, assigning θ = 30 o did not modify the final slope significantly and the FIG. 13. Example of path-length (li) and fidelity (fi) calculation for a generic particle pi.

boundary wall effect is similar to the reference case (see Fig. 12d vs. a). On the contrary, relevant differences were observed with a low solid density (Fig. 12e vs. a): a large portion of the initial column remained stable, glued to the rear walls while only a small curved wedge detached subjected to gravity. A similar result was obtained reducing the thickness of the model of one half (from W/D = 50 to 25): the final slope showed a vertical narrow column constrained on the rear and a smaller portion collapsed at the toe (Fig. 12f).

All these results are consistent with the expected effect of the boundary walls: we can assume that with capillary bonds the coulombian maximum tangential forces at the contact are higher than in dry case. It means that maximum shear stress is higher. Given that the yield stress is generally higher in wet conditions, the incipient movement of the mass (for low shear rate, at the top of the column) will be favoured in dry than in wet conditions, while for larger shear rate (in the fully-developed phase of the collapse, at the toe) the two rheologies tend to be similar. This effect is larger in the upper part and lower at the toe where inertial forces prevail and a flow regime quite close to the collisional one produces frequent capillary debonding events. The solid density affects the normal and tangential stresses experienced by particle-particle and particle-wall contacts: the lower the density, the lower the mobilized stresses and the chance of local failures.

The stronger effect of the wall-particles meniscus than of the particle-particle meniscus can be explained by computing the capillary force magnitude for such cases. Following the same assumptions as in the in Appendix (relative the sphere-sphere contact) and assuming a gap s = 0, the capillary force between a sphere and a flat surface becomes twice the capillary force between two equal particles. This explains the higher probability to have particles "glued" to the lateral walls and the fact that the weakest failure plane tends to move inwards further modifying the velocity profile in transversal direction. From a qualitative point of view, for 2 mm spheres, numerical simulations predict that some triplets of particles may remain stick to the walls; similar structures were not present for coarser particles. The same phenomenon was observed in the experiments, where 2 mm spheres were found to form triplets in contact with the wall, while 3 and 5 mm ones did not.

E. Analysis of the paths and of the bond durability

The dynamics of the collapse is usually deduced from the slope profiles in terms of velocity of the front and evolution of the shape, or from lateral views that allow the measurement of the velocity field.

With the aim of reconstructing the history of the system in more detail, two new simple indexes able to track the motion of the particles and the "quality" of the particle-particle relationships will be introduced. Firstly, let's consider the final particle positions with respect to the initial ones. As sketched in Figure 13, we defined the path-length l i along x-axis as the x-distance between the final and the initial position of each particle; this measure helps us to recognize the static zones, the distribution of mobility of the mass in the runout direction and zones with affine motion.

Figures 14a-c depict the path-length of particles in a middle section 1-cm thick, plotted with respect to the initial particle position. Generally, particles near the movable wall at a middle-low height (in the darkest region) experience the larger path-length with respect to the others, since they are squeezed (compression for overburdened and not confined stress) and projected along x-direction 17 . During their motion the behaviour is almost collisional. Differently, on the upper part of the column the motion is more ordered and similar to a "pure shear". This can be argued noting the path-length value distribution: larger at the surface of the column with a regular decrease with depth and a little clockwise rotation of the gradient due to the progressive evolution of the slope surface.

In particular Figure 14a-c displays 3 cases with D = 2mm: dry, wet and half density. The globally mobilized wedge (with l i > 0) is very different in the three cases while the "squeezed zone" is not so much. In the dry case the wedge touches the rear wall while in the other two cases it does not, meaning that the mobilized zone was reduced. Moreover the sliding surface became steep and curved especially in the half density wet case. It is also evident that the half density case presented a quite rigid wedge at the surface, represented by a quite constant value of the path-length.

In order to combine particle proximity information with the durability of neighbouring relationships, let's now define the particle fidelity coefficient f i : for each particle we identify the near "fiducial particles" as the spheres that are inside a bounding spherical region of radius R b , centred on the considered particle, in the initial state and in the final stage (after the collapse). The ratio between the number of neighbouring particles (inside the bounding sphere) in the initial state over the number of the identical particles in the final condition was defined "fidelity" An example of fidelity calculation was reported in Figure 13. By definition this number spans from zero to one: it is equal to 1 if all the neighbour particles at the beginning are the same in the final state and zero otherwise.

It should be stated that particle fidelity depends on strain level with the advantage that is not sensitive to rotations. If convection phenomena are present the particle fidelity can have localized maxima (clustered particles).

As a consequence, higher fidelity values may be expected in wet than in dry case: during shear, couples of particles connected by capillary bonds with larger attractive forces remain in touch if the disturbance is not very high.

Again, Figure 15a-c displays 3 cases: dry, wet and half density. As expected, comparing Figure 15a and b it is evident that higher fidelity zones are more frequent in the wet than in the dry case. Especially in the wet conditions, minimum values of fidelity localize in a narrow surface that coincides with the failure surface predicted by the path-length distribution and the triangle over this surface seems to move rigidly. Instead, in the dry case the failure surface was not so clearly defined since every particles moved independently. The half solid density case (Fig. 15c) provided a further reduction of the mobilized zone and increased the frequency of high fidelity values. This supported the observation that in the wet case, especially for low solid density, the motion of the unstable wedge was quite rigid while frequent bond breakages were experienced by the particles in the shear plane, between the stable portion of the column and the rigid wedge.

VI. CONCLUSIONS

A micromechanical model for granular materials in the pendular state was defined. It incorporated a capillary force model that derived from the minimum energy approach, a physical law for the bonding-debonding mechanism and a model that stored the little droplet on the surface of the particles and that followed their motion.

Discrete Element simulations of column collapse tests were successfully performed to mimic 1:1 scale experiments on wet and dry spheres. The comparison resulted qualitatively good for different diameters and liquid tension. The general trends was well captured and the important effect of the presence of liquid was observed: wet conditions presented steeper final profiles while the reduction of the liquid tension or the increase in particle diameter tended to reduce this effect. All these aspects can be explained in terms of the ratio between body and capillary forces (i.e. Bond number). The estimation of the profiles, particularly at the toe of the slopes, was not completely successful: many particles at the base were free to rolls despite of the existence of a ball-wall liquid bridge. This fact overestimated the runout distance of wet spheres and it was ascribed to the lack of a rolling resistance for the capillary bonds. On the contrary, for the big spheres, runout values resulted underestimated and the behaviour is compatible with the lack of a lubrication force at the contact.

Collapse simulations in wet conditions highlighted that, from a micromechanical standpoint, the final slope can be divided into three different zones: the static zone (under the failure plane), the upper part, near the rear wall, and the toe of the slope. The top of the slope experienced "quasi-static" phenomena that highlighted the capillary force effects between particles and between particles and lateral walls. Conversely, at the toe, the dynamics was ruled by unbalanced forces that destroyed the capillary bonds minimizing their effect.

As consequence, greater differences between wet and dry states were observed in the rear of the slope where low shear rates prevailed in almost all the test duration. On the other hand, at the toe of the slope, the repose profiles were not so different for different liquid contents. Furthermore, in agreement with experiments, considering larger particle size, the capillary forces increased lower than the body forces. This explained why, in simulations and laboratory experiments, the larger the particle radius the lower the distance between wet and dry profiles.

Many variables influenced the final profile of the wet material: a part from the moisture content and particle size, some of them can be linked to the finite size effect (the thickness of the channel with respect to the grain size), others to the ratio between inertial and capillary forces (liquid tension, contact angle, solid density).

Other variables instead showed very small effects: the liquid-solid contact angle (in a reasonable range of values) and the mechanism of wall removal seemed not influential. Using virtual conditions like the half solid density case was useful to predict the collapse mechanism at different gravities and to appreciate different ratios between capillary and body forces.

In order to recognize the internal mechanism of collapse we used a simple quantity, called path-length, that tracked each particle from its initial position to the final "destination". Following this strategy we were able to observe that particles positioned at a middle height of the column and near the movable wall reached the largest runout distance: at the first stages of the collapse they were squeezed in horizontal direction, like in a plug flow velocity field. This important aspect as well as the observation of different regimes (quasi-static and collisional) during time in different zones of the mass are some key issues that should be taken into account when passing to continuum approaches in which classically the depthaverage assumption is employed.

To investigate the cohesive behaviour of the material and the possible aggregation of the particles due to capillary bonds we introduced a new particle quantity called "fidelity". The value of this variable describes the level of damage of the bonds near each particle: no clusters formed in dry conditions with exception of the static zone. In wet case, instead, fidelity revealed the presence of many pendular clusters that were concentrated on the upper part of the column where rigid-like mechanism were observed. Furthermore, low values of fidelity showed the location of the failure plane, whose angle was steeper and well defined in wet case.

The adopted model gave good qualitative results which can be further improved, also quantitatively, once the rolling resistance and contact lubrication will be incorporated. In order to complete the understanding of the collapse in partial saturation conditions, the effects of initial capillary distribution and of higher water content (funicular state) need to be addressed in the future.

FIG. 1 .

 1 FIG. 1. Comparison of some capillary force models for a liquid bridge between two equal spheres (D = 2 mm, γ = 72.75 mN/m, s * = s/D, F * = Fcap/(ρgV sph ), V * = V /V sph )

FIG. 2 .

 2 FIG. 2. Capillary bridge geometry.

FIG. 3 .

 3 FIG.3. Layout of a collapse test.

FIG. 4 .FIG. 5 .

 45 FIG. 4. Snapshots of the final profiles for D = 2 mm: (a) dry conditions; (b) wet conditions (γ = 72.75 mN/m; liquid content w = 0.5%).

FIG. 6 .FIG. 7 .

 67 FIG. 6. Simulated final profiles for D = 2 mm; (a) dry; (b) γ = 72.75 mN/m (water); liquid content w = 0.5%.

FIG. 8 .

 8 FIG. 8. Comparison of the normalized runout length from experiments (Exp) and DEM simulations (Sim) using (a) water (γ = 72.75 mN/m) and (b) water + surfactant (γ = 17.1 mN/m).

FIG. 9 .

 9 FIG. 9. Slope profile for different simulations: reference simulation (dashed line) [D = 2 mm; γ = 72.75 mN/m (water); liquid content w = 1%; contact angle θ = 15 o ; solid density ρs = 2532 kg/m 3 ; wall lifted with acceleration = g] vs a simulation (solid line) with: (a) coarse particles (D = 5 mm); (b) γ = 17.2 mN/m (water+surfactant); (c) contact angle θ = 30 o ; (d) ρs = 1266 kg/m 3 (half of the glass density); (e) initial thickness W = 2.5 cm;

FIG. 10 .

 10 FIG. 10. Simulated final profiles for different water content w (D = 2 mm, γ = 72.75 mN/m (water)).

FIG. 12 .

 12 FIG. 12. Contour lines of the surface for: (a) the reference simulation [D = 2 mm; γ = 72.75 mN/m (water); liquid content w = 1%; contact angle θ = 15 o ; solid density ρs = 2532 kg/m 3 (glass); wall lifted with acceleration = g] (b) as in (a) but with coarse particles (D = 5 mm); (c) as in (a) but with γ = 17.2 mN/m (water+surfactant); (d) as in (a) but with contact angle θ = 30 o ; (e) as in (a) but with ρs = 1266 kg/m 3 (half of the glass density); (f) as in (a) but with initial thickness W = 2.5 cm; (g) as in (a) but with liquid content w = 0% (dry);

3 FIG. 14 .FIG. 15 .

 31415 FIG. 14. Path-length values for: (a) D = 2 mm, dry glass sphere (ρs = 2532 kg/m 3 ); (b) D = 2 mm, wet glass sphere (γ = 72.75 mN/m, w = 1%, ρs = 2532 kg/m 3 ); (c) D = 2 mm, wet sphere half solid density of glass (γ = 72.75 mN/m; w = 1%; ρs = 1266 kg/m 3 ).

FIG. 16 .

 16 FIG. 16. (a) Failure and (b) formation of a capillary bridge from a single drop and (c) from two drops.

TABLE I .

 I Measured liquid parameters

	Liquid type	Surface tension a γ [mN/m]	Contact angle b θ [ o ]
	water	72.75	15
	water + surfactant	17.1	40

TABLE II .

 II Parameters used in DEM simulations.

	Solid	Normal	Tangential	Contact	Norm. and tang.
	density a	stiffness b	stiffness b	friction coeff. b	viscous damping coeff. c
	ρs [kg/m 3 ]	kn [kN/m]	ks [kN/m]	tan φµ	n
	2532	400	100	0.62	0.1

a Measured by liquid picnometry. b Calibrated with triaxial tests. c Estimated and then calibrated on the dry collapse tests.

Appendix A: The capillary force

The total energy to be spent to form a liquid bridge between two equal spheres is equal to the difference between the energy required to form the solid-liquid interface and that lost at the solid-air interface 46 :

where γ sl and γ sv are the solid-air and solid-liquid interfacial tension.

For small s and β, and considering the validity of the Young-Dupre equation (γ cos θ = γ sv -γ sl ) the Eq. A1 reduced to [START_REF] Israelachvili | Intermolecular and surface forces / Jacob N. Israelachvili[END_REF] :

The attractive force results as:

The meniscus volume may be represented in a cylindrical geometry by:

assuming small β and the wet profile of the sphere as a parabola (b ≈ Rβ 2 /2) the liquid volume becomes 30 :

Assuming the invariance of the liquid bridge volume with the inter-particle distance (dV /ds = 0), one obtains:

and finally:

where b(s, V ) represents the wet spherical segment on particle surface (Fig. 2), defined as:

Appendix B: The capillary rupture and formation

This part of the algorithm manages capillary bridge rupture, redistribution and formation which are expected for large deformation of the wet granular material. These mechanisms are sketched in Figure 16 for clarity.

When two particles reached the rupture distance s f the liquid was splitted into two equal volume drops and assigned to the two particles while the originary capillary volume of that specific potential contact was set to zero. The effect of gravity on the redistribution was neglected.

At each time-step a simple liquid drop kinematics on the particle surface was considered: the motion of each drop was associated to the motion (displacement and rotation) of the relative particle at which it was attached neglecting the variation of the particle centre of mass.

The capillary bridge formation was ruled as in the following: two single drops of volume V 1 and V 2 , were assumed to lay on the surface of two different particles of radius R with two different "direction" u 1 and u 2 (see Fig. 16b,c). β 1 and β 2 were the two positive "filling angles" of the drops (that in general are different from the "filling angles" of the capillary bridges). The normalized vector joining the center of particle i with particle j was dij . When the two particles approached, the formation of a new liquid bridge could happen only if d ij -2R ≤ s c . As second step, the "directions" of the two drops u 1 and u 2 and their filling angles β 1 and β 2 were checked with respect to the unit vector dij leading to three possibilities: no formation of capillary volumes; the formation of one bridge of volume V = V 1 (or V = V 2 ) from a single drop (see Fig. 16b); the formation of one bridge of volume V = V 1 + V 2 when the two drops coalesced (see Fig. 16c).

The angles β 1 and β 2 should be related to the volume of the sessile drops, the particle radius and the value of the contact angles. However, for the sake of simplicity, it was assumed a fixed drop-filling angle ψ = β 1 = β 2 = 45 o as a constant in the model.


R Yang
email: yang.rui@ensta-bretagne.fr

J G Wang

B Clement
email: benoit.clement@ensta-bretagne.fr

A Mansour
email: mansour@ieee.org

FPGA Implementation of a Parameterized Fourier Synthesizer

Field-Programmable Gate Array (FPGA) offers advantages for many applications, particularly where missions are complex and time performance is critical. For small-production digital acoustic synthesizers, FPGA can achieve the abovementioned tighter system requirements with low total system costs on single chip. In this manuscript, a real-time acoustic synthesizer is implemented using Fourier series algorithm on Altera's Cyclone II FPGA chip. This work emphasizes systematic designs and parallel computations. The proposed system includes a flexible processor and a parallel parameterized acoustic module. On one hand, the Nios II embedded processor, which is relatively low-speed component, is used to generate commands and configure high-speed acoustic module parameters. On the other hand, acoustic module which should require high-speed components contains 4 parallel architectures to gain high-speed simultaneous calculus of 4 independent digital timbres. Every timbre is equivalent to 16 parallel high-precision harmonic channels with 0.3 % frequency error. Experimental results corroborate the fact that a single FPGA chip can achieve complex missions and attain real-time performances.

I. INTRODUCTION

Since the last two decades, FPGA has become widely used in electronic engineering applications after its first product in 1984. Comparing with Application-Specific Integrated Circuit (ASICs), FPGA is considered as the best technical solution for small-production applications. According to [START_REF] Kuon | Measuring the Gap between FPGAs and ASICs[END_REF], first FPGA circuits are slower, less energy efficient and larger size than fixed ASICs. However, FPGA has many advantages: Flexibility, reconfigurability and low-cost. New technology and devices enhance these features and make FPGA more and more popular in many electronic applications.

It is well known that FPGA has the ability to implement any combinational function implemented on ASICs. Programming FPGA to form complex system logical functions can be achieved by interconnecting Look-Up Table (LUT) elements on a single chip [START_REF] Compton | Reconfigurable computing: A survery of system and software[END_REF].Comparing to digital computer architectures, such as Von Neumann architecture or Harvard architecture [START_REF] Smith | The Scientist and Engineer's Guide to Digital Signal Processing[END_REF], FPGA can provide more flexibility and achieve better calculation performance. On one hand, FPGA allows designers to make serial and parallel mixed designs that are more flexible for complex missions. On the other hand, parallel designs can produce more simultaneous calculus for highspeed applications. In this paper, we focus on the case of digital acoustic synthesizers.

The primary aim of a synthesizer is to generate different acoustic timbres. Each timbre can be distinguished from others by the quantity of different harmonics. Acoustic synthesizers can be divided into two kinds: Analog and digital synthesizers.

Actually, analog synthesizer applications have good sound quality and are based on many well-defined algorithms, such as Frequency Modulation (FM) which is popular and simple to be implemented on electronic devices. Recently, more convenient digital approach emerges with similar ideas. Main advantages of digital synthesizers comparing analog ones are robustness, cheap prices and simple architectures.

Two popular digital synthesizers focus on Direct Digital Synthesizer (DDS) [START_REF] Kroupa | Direct Digital Frequency Synthesizers[END_REF] and FM synthesizers [START_REF] Chowing | The synthesis of complex audio spectra by means of frequency modulation[END_REF]. The crucial technique in DDS is the phase-amplitude conversion architecture [START_REF] Vankka | Digital Synthesizers and Transmitters for Software Radio[END_REF]. The main drawback of DDS is still the larger size of needed ROM tables for higher performance. FM approaches do not suffer from these inconvenient and they intend to modulate the frequency of carrier wave for generating the timbres. In fact, FM methods can save more hardware resources than DDS. The major drawback of FM techniques is still the need for an expert to tune the multitude parameters. To evaluate the ability of FPGA to handle a large amount of data, we select the Fourier series, namely an extended DDS application, instead of FM which can be easily implemented on FPGA. Our work challenges FPGA on special aspect of solving complex missions and parallel calculations. It is worth mentioning that the Fourier series designs don't require the knowledge of experts and parameters can be tuned automatically.

Using Fourier series, we can decompose a target timbre F (2πf 0 t) as follows:

F (2πf 0 t) = A 0 + ∞ i=1 A i sin(2πif 0 t + ϕ i ) (1)
where f 0 represents the frequency of the target timbre,A i and ϕ i represent respectively the amplitude and phase of i th harmonic. A 0 is the average signal value. From the engineering point of view, equation (1) need to be truncated at a finite order of harmonic instead of using infinite frequency components.

In our case, all the high order harmonics that have less energy than human hearing threshold should be ignored in the synthesizer implementation. According to [START_REF] Qiwen | Principle of Keyboard Maintenance[END_REF], the sum in equation ( 1) can be truncated in worst case up to the 10 th harmonic. To get better precision, we assumed the sum truncated at the 16 th harmonic is rich enough to generate a single timbre.

In order to implement Fourier series on FPGA boards, we designed a complete system which includes three main parts: Human-Machine Interface, Coordinator Module (Embedded Nios II Processor [START_REF]Altera Inc. Nios II Software Developer's Handbook[END_REF]) and Acoustic Module. The Nios II Processor and the Acoustic Module are implemented on a single FPGA chip; further details are given in the following section.

II. SYSTEM DESIGN

As it was early mentioned, the implementation of a complex application, such as acoustic synthesizer, on hardware components is still very challenging. Serial as well as parallel structures are both required. In the following, the whole system architecture is described and the Timbre Module using parallel structure and requiring high computational efforts is emphasized.

Nios II

I2C Bus

Cyclone II FPGA EP2C20

A. System Architecture

In order to generate 16 harmonics at each of the implemented timbres, our system can be divided into three major parts, (the first one is realized using a PC, the two other using FPGA), see Figure 1: 1) Human-Machine Interface: This part is done using a classic PC architecture and it should handle the human commands, sound generation and recording,sound analysis tools, etc. The synthesizer user interface, Commands and Parameters Panel, is developed using a VC++ software to monitor and manage the synthesizer (For example, users can modify frequency or amplitude parameters on-line via this interface). The audio analysis tools have been done using available commercial software, i.e. Adobe Audition. 2) Embedded Nios II Processor: Flexible Nios II embedded processor is relatively low-speed comparing to parallel architecture as it runs codes in serial sequence. In fact, its ultimate processing speed is limited by the CPU clock rate. However, Nios II shows better flexibility and robustness in responding to random environmental changes with interrupt mechanism which are ideally for human-machine interface tasks. In our system, Nios II is dedicated to communicate with the PC, to manage audio codec, to respond to keyboards and to configure acoustic module. 3) Acoustic Module: High-speed Acoustic Module can reach real time performance by mean of parallel architectures which can solve large calculations in short time.

In the Acoustic Module we have 4 independent Timbre Modules which simultaneously generate 4 different timbres. The real time timbre data is processed in the Polyphonic Module to get a 32 bits polyphonic audio stream. The I2S Driver converts digital audio stream into serial I2S [START_REF]I2S bus specification[END_REF] format and finally analog audio signal is generated using an audio codec (such as WM8978).

B. Timbre Module

Previously, we mentioned that a timbre can be synthesized as a truncated sum of equation [START_REF] Kuon | Measuring the Gap between FPGAs and ASICs[END_REF]. In this section, we focus on the realization of the j th timbre by electronic circuits. By adjusting A 0 to be zero, adding the shape of ADSR (µ j (t), Attack-Decay-Sustain-Release envelope techniques are used to reshape signal in time domain to mimic the acoustic sounds of real instruments.) and multiplying by the strength parameters (σ j ), equation ( 1) can be rewritten as:

V j = µ j (t)σ j 16 i=1 α ij sin(2πf ij t + ϕ ij ) = µ j (t)σ j T j (2)
In this case, each timbre is characterized by three main features: a) Harmonics: Up to 16 harmonics. b) Strength (σ j ) is equivalent to volume parameters. c) Attack-Decay-Sustain-Release envelope µ j (t). Each of the above features is implemented in different stage of the pipelined Timbre Module, see Figure 2. Concerning the first feature, we need to generate T j , as described in equation (2). To generate the main sinusoidal wave, we may use paralleled conventional DDS concept as shown in Figure 3. In conventional DDS method, Frequency Divider is a digital counter that divides system clock by a chosen factor into desired frequency. Phase Accumulator is another digital counter that generates addresses to activate the RAM table which used to generate Sin Wave Table (SWT). The output T j of the DDS component is obtained as the sum of all partial outputs. It is worth mentioning that the above structure becomes very expensive as it used 16 SWT and 16 multipliers per Timbre Module. This problem is a major one because the FPGA resources are limited.

Phase Accumulator

To solve this problem we proposed another time-division multiplexing (TDM) structure to share one SWT for all the 16 parallel channels. As shown in Figure 2 (a), Frequency Dividers and Phase Accumulators still work in parallel mode, because they require high speed processing and their operating clock rates are usually comparable to the Nios II processor clock. SWTs of the conventional DDS on the contrary run at comparatively lower speed. They are reduced into one sequential access table in our TDM architecture. In order to prevent access conflict among channels, a Table Router was added to control the access priority. The added router can solve half of the problem, because using the router we can't obtain the total sum of the sinusoidal waves as described in equation ( 2). The total sum can be however obtained by using the Multiply Accumulator (MAC), seeing in previous Figure 2 (a). In our architecture, each Timbre Module will only consume one SWT , one extra small size Amplitude Tuning Table (ATT) and one MAC.

To provide the details of this parallel and sequential mixed structure, we separate the TDM architecture into Process Control Channel (PCC) and Data Processing Channel (DPC), seeing in Figure 4.

Finite State Machine

Frequency Dividers & Phase Accumulators

Sin Wave PCC consists of a finite state machine, its function is detecting and coordinating the workflow of DPC. It gathers states information of the DPC, especially the Flags posted by the parallel counting channels. This Flag is triggered to apply for access of Sin Wave Table . When PCC recieve a Flag, PCC send coordinating orders to MUX, SWT, ATT and MAC in the piplined DPC to make a calculation cycle.

DPC is used to receive orders from PCC and process data as fast as possible. It separates into a parallel region and sequential region, they deal with high speed missions and low speed missions respectively. 16 parameteried channels run simultaneously and independently. This design gains more speed in calculating the sum in equation [START_REF] Compton | Reconfigurable computing: A survery of system and software[END_REF]. Whenever a channel finishes counting with a result of table address, it will require PCC to make a sequential calculation cycle. PCC is the coordinator and core mechanism connecting parallel with serial.

The PCC follows the mechanism that low frequency channels get higher priority to Sin Wave Table . In each calculation cycle, MUX connect one channel to SWT, and corresponding Phase Accumulator provide an address pointing to the SWT memory unit, which is current normalized sinusoidal amplitude. Because every harmonic channels share the same sine wave that discreted into 64 points in length and 16 bits in amplitude. Therefore, ATT is essential to provide different amplitude for different harmonics. This coefficient from ATT is multipied with SWT result and storaged in the corresponding buffer in MAC. Then PCC ask MAC to update the sum to the following circuits. After the PCC clear the responded Flag, it will be ready for the next cycle. This kind of pipeline structure is also used to add Strength and ADSR parameters.

.

III. EXPERIMENTAL RESULTS

Many experiments have been conducted in order to evaluate operational properties and performances of the proposed Fourier synthesizer (such as high-speed, high-precision, realtime and flexibility). We also implemented a test board according to Figure 1 and Tab I . FPGA design software Quartus II 9.0 [START_REF]Altera Inc. Quartus II Handbook[END_REF] is used to make a logic synthesis. Step-by-step experiments were carried out to test frequencyamplitude accuracy and the functionality of our ADSR. In the test of frequency-amplitude accuracy,we compared the output signal frequency spectrum with an off-line pre-prepared target signal. The match between the two signals is accepted if they share the same first five harmonics, as shown in Figure 5 that indicates a good accuracy in both frequency and amplitude of each harmonic. In this paper a systematic design of Fourier Synthesizer is proposed and implemented on a single FPGA chip. The proposed circuit, especially the parallel design of Timbre Modules, shows efficient flexibility and large calculation capabilities in achieving real-time, high-speed and high-precision applications, as the implementation of digital acoustic synthesizer in FPGA field. We should highlight the fact that the small-size and low-cost FPGA boards are efficient to perform small-producion applications where tasks are complex and time performance is critical.

Fig. 1 .

 1 Fig. 1. System Diagram of Fourier Synthesizer

Fig. 3 .

 3 Fig. 3. Simple Fourier Synthesis by Sinusoidal DDS

Fig. 5 .

 5 Fig. 5. Frequency accuracy analysis

Fig. 6 .

 6 Fig. 6. Time domain signal and its frequency spectrum

Table Depth

 Depth

			MAC
		: 64		Tj
			Buffer 1
	Frequency Dividers & Phase Accumulators	Amplitude Tuning Table	Buffer 16	DPC
	States/Flags		Commands	PCC

MUX

Fig. 4. SWT and MAC Shared Fourier Synthesis

R. Yang is a joint PhD student between ENSTA Bretagne and Ocean University of China.

This work was supported by the China Scholarship Council.