Mohammad Ali Nazari

Pascal Perrier

Yohan Payan

The Distributed Lambda (λ) λ) λ) λ) Model (DLM): A 3-D, Finite-Element Muscle Model Based on Feldman's λ λ λ λ Model; Assessment of Orofacial Gestures

Keywords: Equilibrium Point Hypothesis, Feldman's λ Model, Muscle active force, Muscle passive force, finite element method, Speech motor control, Biomechanical orofacial model

Purpose:

The authors aimed to design a distributed Lambda model (DLM), which is well-adapted to implement three-dimensional (3-D) Finite Element descriptions of muscles.

Method: A muscle element model was designed. Its stress-strain relationships included the active force-length characteristics of the λ model along the muscle fibers, together with the passive properties of muscle tissues in the 3-D space. The muscle element was first assessed using simple geometrical representations of muscles in form of rectangular bars. Then, it was included in a 3-D face model, and its impact on lip protrusion was compared with the impact of a Hill-type muscle model..

Results:

The force-length characteristic associated with the muscle elements matched well with the invariant characteristics of the λ model. The impact of the passive properties was assessed. Isometric force variation and isotonic displacements were modeled. The comparison with a Hill-type model revealed strong similarities in terms of global stress and strain.

Conclusion:

The DLM accounted for the characteristics of the λ model. Biomechanically no clear differences were found between the DLM and a Hill-type model. Accurate evaluations of the λ model, based on the comparison between data and simulations, are now possible with 3-D biomechanical descriptions of the speech articulators because to the DLM.

INTRODUCTION

A number of studies using biomechanical orofacial models showed that the physical properties of the main speech articulators and their interactions with external structures are determining factors in a number of important characteristics of speech production: the prototypical articulatory configurations associated with each sound, the stability of these configurations, and the shape of the articulatory paths and the formant trajectories in transitions between these configurations. [START_REF] Perkell | Properties of tongue help to define vowel categories: hypotheses based on physiologically-oriented modeling[END_REF] found that the stability of the control of the two most frequent vowels in world languages ([i] and [a]) could be largely due to the strong nonlinearities of the relationships between muscle activation and the degree of constriction in the vocal tract. [START_REF] Perrier | Degrees of freedom of tongue movements in speech may be constrained by biomechanics[END_REF] showed that the main directions of tongue deformations observed in speech movements in the mid-sagittal plane -characterized by the well-known front-raising and back-raising factors found by [START_REF] Harshman | Factor analysis of tongue shapes[END_REF] -originated from intrinsic anatomical and biomechanical properties of the tongue muscles. [START_REF] Perrier | Influences of tongue biomechanics on speech movements during the production of velar stop consonants: A modeling study[END_REF] found evidence that the presence and shape of the articulatory loops observed in vowel-velar consonant-vowel sequences in a number of languages are strongly influenced by tongue biomechanics, including its muscular anatomy and contact with the palate. [START_REF] Stavness | Biomechanical modeling of English /r/ variants[END_REF] found that the most frequent /r/ variants in English are those that correspond to the minimum amount of volume displacement, relative strain, and relative muscle stress.

These results emphasize the need to include realistic and reliable biomechanical descriptions of the orofacial motor system in models of speech production. The speech articulators that are the most influential in achieving the fine tuning of the vocal tract's shape, which determines acoustic properties of speech (see [START_REF] Fant | Acoustic theory of speech production[END_REF][START_REF] Stevens | Acoustic phonetics[END_REF], are the tongue, lips, velum, and pharyngeal constrictors.

These articulators are all made of soft tissue, which is mostly muscle tissue. To model the biomechanical behavior of soft tissues, the finite element method [START_REF] Bathe | Finite Element Procedures[END_REF] has been proven to be extremely efficient, accurate, and reliable (e.g., [START_REF] Payan | Soft Tissue Biomechanical Modeling for Computer Assisted Surgery[END_REF]. This method uses a numerical technique to compute an approximate solution to a set of partial differential equations. It relies on a discretization of the continuum domain (Ω) to be simulated. This discretization, called a mesh, is the partition of Ω into simpler geometrical bodies known as elements defined by a set of vertices or nodes. Being able to formulate muscle models in the context of the finite element method is an important challenge for speech production modeling.

In biomechanics, the reference muscle model is the Hill-type model (Mc Mahon, 1984;[START_REF] Zajac | Muscle and tendon: Properties, models, scaling, and application to biomechanics and motor control[END_REF], and a number of finite element formulations of this model have been proposed [START_REF] Weiss | Finite element implementation of incompressible, transversely isotropic hyperelasticity[END_REF][START_REF] Cheng | Virtual muscle: a computational approach to understanding the effects of muscle properties on motor control[END_REF][START_REF] Blemker | A 3d model of muscle reveals the causes of nonuniform strains in the biceps brachii[END_REF] including for models of orofacial articulators [START_REF] Wilhelms-Tricarico | Physiological modeling of speech production: Methods for modeling soft-tissue articulators[END_REF][START_REF] Koolstra | A method to predict muscle control in the kinematically and mechanically indeterminate human masticatory system[END_REF][START_REF] Stavness | Coupled hard-soft tissue simulation with contact and constraints applied to jaw-tongue-hyoid dynamics[END_REF]. However, in motor control research in general, and in speech motor control research in particular, another muscle model, the lambda (λ) model, is often used, mainly because it is embedded in a theory of human motor control, the Equilibrium Point Hypothesis proposed by [START_REF] Feldman | Once more on the equilibrium-point hypothesis (λ model) for motor control[END_REF]. For a number of reasons that will be discussed in the next section, we believe that for speech motor control modeling, the λ model is more appropriate than Hill-type models.

We previously proposed using the λ model in finite element models of speech articulators [START_REF] Payan | Synthesis of v-v sequences with a 2d biomechanical tongue model controlled by the Equilibrium Point Hypothesis[END_REF][START_REF] Sanguineti | A control model of human tongue movements in speech[END_REF][START_REF] Buchaillard | A biomechanical model of cardinal vowel production: Muscle activations and the impact of gravity on tongue positioning[END_REF]; see also [START_REF] Sanguineti | A dynamic biomechanical model for neural control of speech production[END_REF]. In these models, the active part of the muscle was functionally modeled as a set of force generators acting as external forces applied onto the nodes of the finite element structure. However, this functional approach did not account for the fact that the muscles are part of the continuum. To increase the realism of the description, and, in particular, to provide a better account of the active-force-generation mechanisms of a muscle and their consequences on the mechanical properties of muscle tissue, we developed what is, to the best of our knowledge, the first finite-element formulation of the λ model, in which the muscle model is part of the factors that determine the stress-strain characteristic of the tissues.

In this paper, we first briefly describe the λ model and reasons why we consider it to be well-suited to speech production modeling, and then we describe the finite-element formulation of the model. Some classical evaluations of this formulation are then proposed, followed by a preliminary comparison of this formulation with a formulation of a Hill-type model.

THE λ λ λ λ MODEL: A WELL-SUITED MODEL FOR SPEECH MOTOR

CONTROL MODELING

In Hill-type models, the mechanical properties of the muscle (i.e., the relation between force and strain) are based on measures recorded from ex-vivo muscles that are artificially tetanized (i.e., maximally activated) with external electrical stimulations. The force-strain relations for smaller levels of activation are not measured. They are estimated and functionally modeled with various intermediate multiplicative or additive accounts [START_REF] Zajac | Muscle and tendon: Properties, models, scaling, and application to biomechanics and motor control[END_REF][START_REF] Winters | Hill-based muscle models: A systems engineering perspective[END_REF][START_REF] Shapiro | Control variables in mechanical muscle models: A minireview and a new model[END_REF]. In speech production, orofacial muscles generate levels of force that are far from their maximal force reached in tetanized conditions. The control variable of Hill-type models is the level of the force itself. It is known, however, that the actual muscle force is a consequence of a combination of influences due to descending commands from the central nervous system (CNS) and afferent signals associated with muscle length (via muscle spindles) and the rate of change in muscle length (via Golgi tendons) [START_REF] Mcmahon | Muscles, Reflexes, and Locomotion[END_REF]. Speech motor control has been shown to be very resistant to perturbations such as changes of the head position with respect to the gravity field [START_REF] Shiller | Effects of gravitational load on jaw movements in speech[END_REF], or unexpected perturbations of the mandible [START_REF] Folkins | Additional observations on responses to resistive loading of the jaw[END_REF] or the lower lip [START_REF] Abbs | Control of complex motor gestures: orofacial muscle responses to load perturbations of lip during speech[END_REF][START_REF] Gomi | Compensatory articulation during bilabial fricative production by regulating muscle stiffness[END_REF]. We believe that part of the stability of speech motor control is due to low-level feedback from the muscles to the motoneuron pool. Hence, we think that Hill-type models are not the most appropriate muscle models for speech motor control modeling. The λ model seems to be more appropriate, for several reasons.

The λ model includes hypotheses about the nature of the control variables and a description of the muscle-force-generation mechanisms. It is based on a study of physiological measurements of forcelength relationships in in-vivo muscles in deafferented cats (Feldman & Orlosvky, 1972) and on experimental data from human subjects in unloading tasks [START_REF] Feldman | Once more on the equilibrium-point hypothesis (λ model) for motor control[END_REF]. Both types of measurements have been done for different levels of muscle activation. The control variables do not specify force levels directly. The control variables are threshold muscle lengths (λ) above which active force generation begins. For a muscle, the specification of a λ value corresponds to the selection of a specific force-length relationship. Thus, the actual muscle force results from the combined influences of the λ value and of afferent inputs associated with muscle length, rate of change in muscle length, and cutaneous reflexes [START_REF] Feldman | Positional frames of reference in motor control: Origin and use[END_REF][START_REF] Pilon | Threshold position control of arm movement with anticipatory increase in grip force[END_REF]. In a given external force field, the λ value determines the muscle length at which the mechanical equilibrium point of the motor apparatus is reached. According to the Equilibrium Point Hypothesis, movements are controlled by shifting this mechanical equilibrium point. From this perspective, movements are the result of the attraction of the motor apparatus toward the specified equilibrium point, and continuous movements are the consequences of successive displacements towards a discrete sequence of equilibrium points. Such a movement-generation principle suggests the existence of a discrete representation of the motor task in the CNS.

A fundamental consequence of this hypothesis is that, once the CNS has specified the time variation of the λ values of all the muscles included in the motor apparatus, the movement's trajectory and its timing are fully determined by muscle mechanics, including feedback loops from mechanoreceptors, interacting with external dynamical constraints, such as external loads, frictions, or accelerations [START_REF] Ostry | A critical evaluation of the force control hypothesis in motor control[END_REF]. Trajectories are therefore assumed to be the consequences of the specification of the λ values; they are not the motor goals.

By controlling for movements and positions of the motor apparatus in a given external force field via the specification of a sequence of equilibrium points, the Equilibrium Point Hypothesis intrinsically predicts that a large amount of different motor command patterns and muscle forces can be associated with the same sequence of intended spatial positions. The Equilibrium Point Hypothesis accounts for the principle of co-activation of agonist and antagonist muscles, which enables the same mechanical equilibrium point of the motor apparatus to be reached for various levels of force in each muscle.

Different patterns of co-activation associated with the same goal will generate different movements toward the goal, with various velocity profiles and/or various trajectories. To account for this variability, two macro-parameters determine movement in the Equilibrium Point Hypothesis: an R (reciprocal) parameter that specifies the intended equilibrium point and a C (co-activation) parameter that determines the dynamic properties of the motor apparatus and influences the timing and velocity of the movement (i.e., the way to move toward the intended equilibrium point) [START_REF] Feldman | Positional frames of reference in motor control: Origin and use[END_REF].

We believe that the following basic principles of the Equilibrium Point Hypothesis make this theory very suitable for explaining speech motor control:

(1) The discrete representations underlying the generation of movements enable a link between the discrete phonological characterization of the speech sequence and the continuous articulatory and acoustic signals that carry the phonological information from the speaker to the listeners.

(2) The fact that movements are supposed to be intended towards dynamic attractors, defined in terms of mechanical equilibrium points, tends to provide stability to the motor system and to facilitate the preservation of equifinality under various movement conditions (see [START_REF] Feldman | Testing hypotheses and the advancement of science: recent attempts to falsify the Equilibrium Point Hypothesis[END_REF], for a discussion of equifinality in the context of the Equilibrium Point Hypothesis).

(3) The disassociation between the R parameter that specifies the intended equilibrium points and the C parameter that influences the dynamic conditions of the movement is very useful to explain the variability observed for the same sequence of phonemes, pronounced under various speaking-rate or clarity conditions (see for example [START_REF] Matthies | Variation in anticipatory coarticulation with changes in clarity and rate[END_REF].

Using a two-dimensional (2D) biomechanical model of the tongue [START_REF] Payan | Synthesis of v-v sequences with a 2d biomechanical tongue model controlled by the Equilibrium Point Hypothesis[END_REF] that is driven and controlled by the principles of the Equilibrium Point Hypothesis, some interesting properties of speech movements that had been observed in experiments with several subjects were simulated: looping articulatory patterns observed in vowel-velar consonant-vowel sequences [START_REF] Perrier | Influences of tongue biomechanics on speech movements during the production of velar stop consonants: A modeling study[END_REF], variability in articulatory trajectories for the same goals depending on the speaking rate or the stress [START_REF] Perrier | Biomechanical models to study speech[END_REF], and relations between trajectory curvature and velocity [START_REF] Perrier | Speed-curvature relations in speech production challenge the onethird power law[END_REF]. The results all suggested that the Equilibrium Point Hypothesis is well-adapted to model speech motor control.

The Equilibrium Point Hypothesis is also one of the most controversial theories for motor control.

The main criticisms are twofold. A first set of criticisms is about the fact that the Equilibrium Point Hypothesis rejects the hypothesis that the determination of motor commands that are adapted to a given movement would involve, in the brain, an inverse computation using a complex internal dynamic model of the motor apparatus [START_REF] Gomi | Equilibrium-point control hypothesis examined by measured arm stiffness during multijoint movement[END_REF][START_REF] Wolpert | Internal models in the cerebellum[END_REF].

This debate, which is beyond the scope of this paper, has been discussed elsewhere [START_REF] Perrier | About speech motor control complexity[END_REF][START_REF] Perrier | Gesture planning integrating knowledge of the motor apparatus's dynamics: A literature review from motor control and speech motor control[END_REF]. The second set of criticisms is about using the λ model for muscle-force generation, which is intrinsically associated with the Equilibrium Point Hypothesis. In simple terms, those who are opposed to using the λ model claim that simulations of experimental data with this model would require an unrealistically high stiffness in the model [START_REF] Gomi | Equilibrium-point control hypothesis examined by measured arm stiffness during multijoint movement[END_REF][START_REF] Hinder | The case for an internal dynamics model versus equilibrium point control in human movement[END_REF] and high gains in the feedback loops from muscle spindles, Golgi tendons, and cutaneous receptors [START_REF] Wolpert | Internal models in the cerebellum[END_REF]. Feldman and colleagues have argued that these authors may have used biomechanical models of the motor apparatus that were too simple and unrealistic (Gribble, Ostry, Sanguineti & Laboissière, 1998;[START_REF] Feldman | Testing hypotheses and the advancement of science: recent attempts to falsify the Equilibrium Point Hypothesis[END_REF]. Designing a finite element method formulation of the λ model is an important step to resolve this debate. Indeed it will enable extensive and quantitative tests of the λ model under many conditions, by comparing simulations and experimental data.

METHOD

We aimed to design an active muscle element that is suited to the Finite Element Method used to model soft tissues such as the tongue, the lips, the velum or the pharyngeal constrictors. This model has to account for the passive elastic properties and the active force generation mechanisms of a muscle, as described in the λ model. To be able to compare this model with Hill-type models, we also designed a muscle element accounting for the properties of Hill-type models.

Rationales Functional Muscle Models

In functional models, muscle behaviors are represented by macroscopic descriptions of the input/output relationships, without considering the details of the mechanisms underlying force generation [START_REF] Mcmahon | Muscles, Reflexes, and Locomotion[END_REF].

Adjustable-stiffness models

Basic Hill-type models belong to this group. In this model, the total muscle force (F m ) is expressed as the sum of a force in a parallel elastic element (F PE ) and a force in a contractile element (F CE ):

F m =F PE +F CE (1)
The force in the contractile element is a function of muscle length (l), muscle velocity (v), time (t), and a control variable called henceforth A c :

F CE =f(l,v,A c ,t) (2)
Classically, in Hill-type models [START_REF] Zajac | Muscle and tendon: Properties, models, scaling, and application to biomechanics and motor control[END_REF][START_REF] Cheng | Virtual muscle: a computational approach to understanding the effects of muscle properties on motor control[END_REF], this force is expressed in a multiplicative way as a product of three distinct and independent functions: a force-length function (F L ), a force-velocity function (F v ), and a time-varying, centrally specified activation function (called henceforth activation dynamics) (f ac ):

F CE =F L (l)×F v (v)×f ac (A c ,t) (3)
where F L and F v account for intrinsic biomechanical properties of the muscles, and f ac accounts for the time-varying central contribution to the level of the control variable.

In Hill-type models, the force-length characteristic of a muscle, F L (l), is based on a series of measurements obtained in isometric conditions for different lengths, when the muscle produces its maximum voluntary force (MVF) [START_REF] Zajac | Muscle and tendon: Properties, models, scaling, and application to biomechanics and motor control[END_REF]. The shape of this curve is parabolic and is similar to the force-length characteristic of the sarcomeres (Figure 1, thick solid line). In order to describe the contractile force at the same length but for smaller levels of activation, the usual method consists of multiplying the MVF-length characteristic by a factor that is less than 1 (Figure 1, dotted curves).

-

--------------------(Figure 1 around here) ------------------------
This multiplicative account of the contractile force can be categorized as an adjustable-stiffness model [START_REF] Shadmehr | A mathematical analysis of the force-stiffness characteristics of muscles in control of a single joint system[END_REF] in which the contractile element is assumed to behave like a nonlinear spring whose stiffness varies depending on muscle activation, rate of muscle length change, and time (see Appendix 1).

Adjustable-starting-length models: Feldman's model

In other muscle models, the contractile force is assumed to be a general nonlinear function of A c , muscle length, and velocity. [START_REF] Shapiro | Control variables in mechanical muscle models: A minireview and a new model[END_REF] proposed a multiplicative representation in which the control variable and the muscle length are not separable variables):

F CE =h(A c ,l,t)×g(v) (4)
In adjustable-starting-length models, the control variable A c is specified as a length quantity [START_REF] Shadmehr | A mathematical analysis of the force-stiffness characteristics of muscles in control of a single joint system[END_REF]. It determines the function h by specifying the muscle length value l for which the active force is equal to zero. This is the case for the λ model [START_REF] Feldman | Functional tuning of the nervous system with control of movement or maintenance of a steady posture -ii. controllable parameters of the muscles[END_REF][START_REF] Feldman | Once more on the equilibrium-point hypothesis (λ model) for motor control[END_REF]. In such models, in contrast to the Hill-type model (see Appendix 1), muscle stiffness is not directly controlled: it also depends on muscle length (l) and the rate of muscle-length change (v), due to reflex loops associated with mechanical receptors such as muscle spindles and Golgi tendons.

The λ model states that the muscle starting length (i.e., the zero-force-point length) λ is the activation command centrally specified by the CNS. For a given λ value, muscle behaves under the influence of the stretch-reflex mechanism. Based on experimental measurements from a deafferented cat gastrocnemius muscle [START_REF] Feldman | The influence of different descending systems on the tonic stretch reflex in the cat[END_REF], Feldman proposed that the stretch reflex mechanism is associated with a force-length characteristic in the form of an exponential curve, called the invariant characteristic [START_REF] Feldman | Once more on the equilibrium-point hypothesis (λ model) for motor control[END_REF]:

F active_Feldman =F max (exp([l(t-t d )-λ * +µv(t-t d )] + /l c )-1) (5) 1
where F max is the maximum force generated by the muscle, l(t) is the muscle length at time t, λ * is the starting (also called "threshold") length (see Appendix 2), l c is a characteristic length, v(t) is the rate of ---------------------(Figure 2 around here) ------------------------

To account for the sliding filament theory [START_REF] Huxley | Muscle structure and theories of contraction[END_REF] the force described by the IC curves is multiplied by a term that varies between 0 and 1 as a function of the rate of muscle-length change (see Appendix 3).

Comparison of the theoretical characteristics of adjustable stiffness (Hill-type) and adjustable starting-length muscle models

To study the difference between adjustable stiffness Hill-type models and the λ model, we plotted for each type of muscle model a set of muscle force-length curves for different activation levels (Figure 3). Let us now consider an arbitrary example of possible force and length variations of a muscle that would be moving voluntarily against an external load. The hypothetical length variation is represented in Figure 3 by the thick, bold black line that connects the resting length value of the muscle (length = 1) to the length value 0.6, when the force varies from 0 to 0.35 N. This thick line crosses the forcelength characteristic of both types of muscle model, and the intersection points are very close across models. Hence, the displacement of the external load corresponding to this path in the force-length plane is achievable by both types of muscle models. This can be done with realistic changes in the centrally specified control variable, namely the multiplying factor A c for Hill-type models and the threshold length λ for the λ model, in a way that generates similar movement patterns both in terms of displacement and muscle force.

---------------------(Figure 3 around here) ------------------------

Summary

In Hill-type models, the force is directly controlled with a control variable that scales the reference curve (Figure 1), and the stiffness is linearly related to the control variable. In the λ model, the control variable λ specifies the force-length characteristic (Figure 2). As a result, muscle force and muscle stiffness are under the combined influence of the control variable, the muscle length, and the rate of change in muscle length. In both kinds of models, the active force-length characteristic is combined in an additive way with the passive force length characteristic of muscle tissues, and the sliding-filament theory is accounted for with a non-linear function describing the relationship between force and muscle--change rate.

Hence, the principles underlying the choice of the control variables and the way muscle force and muscle stiffness are controlled are clearly different in these two models. The mechanical characteristics are based on two different kinds of experimental data. The comparison of their respective static force-length characteristics, as shown by their superimposition in Figure 3, does not show large differences in the range of force levels (moderate) and muscle length (somewhat smaller than the rest length) commonly used in speech production. This observation suggests that their respective mechanical impacts on tongue, lip and velum shapes and positions should not be very different.

Finite Element Modeling of the λ λ λ λ Model: Distributed Lambda Model

(DLM)

The finite element method requires a discretization of the soft body in a number of small elements.

Hence, formulating the λ model in such a modeling framework necessitates moving from the lumped original description of the model (equation [5] and Appendix 3, Equation [9]) to a "distributed description" in the direction of the muscle fibers. In a distributed model, all lumped quantities of the physical object are replaced with distributed representations:

-The force term is replaced with the Cauchy stress σ, which is the limiting value of the ratio of the force ∆F along the muscle-fiber direction to the cross-sectional area ∆A of the muscle (σ=

lim(∆F⁄∆A) ∆A→ 0 ).
-Length is replaced with the stretch ratio value, which is the ratio of the current muscle length to its initial length (SR=l/l 0 ).

Starting from equation ( 5), the active Cauchy stress becomes:

σ active_Feldman =σ max (A pcsa /A)(exp([SR(t-t d )-SR threshold +µ (t-t d )] + (l 0 /l c ))-1) (6)
In this relationship σ max =F max /A pcsa , where A pcsa is the physical cross-sectional area of the muscle, and σ max is the maximum stress generated by the muscle. The threshold length λ * is replaced with the threshold stretch ratio (SR threshold = λ * /l 0 ) and the velocity term becomes the strain rate ( =v/l 0 ).

With this new expression some parameters of the original λ model take on new meanings. For example, the ratio l 0 /l c shows how long the resting length l 0 is in comparison to the characteristic length l c . According to equation ( 5), the l c value influences the derivative of the invariant characteristics: the smaller the l c value, the stronger the derivative of the invariant characteristics. The derivative of the invariant characteristics corresponds to the stiffness of the active part of the muscle in the muscle-fiber direction due to the stretch reflex mechanism, when the muscle is stretched above the threshold length λ. Therefore, the ratio l 0 /l c can be physically interpreted as follows: the longer the resting length in comparison to the characteristic length, the stronger the stiffness of the stretch reflex, and the faster the attraction movement toward the equilibrium position determined by the resting length.

In the proposed three-dimensional (3D) distributed lambda (λ) model (DLM), a muscle is represented with a set of volumetric elements in which muscle fibers are embedded in a matrix of surrounding passive tissues. These surrounding passive tissues are modeled with hyperelastic materials in agreement with previous approaches [START_REF] Wilhelms-Tricarico | Physiological modeling of speech production: Methods for modeling soft-tissue articulators[END_REF][START_REF] Gerard | Non linear elastic properties of the lingual and facial tissues assessed by indentation technique. Application to the biomechanics of speech production[END_REF]. In line with our previous work [START_REF] Nazari | Simulation of dynamic orofacial movements using a constitutive law varying with muscle activation[END_REF][START_REF] Buchaillard | A biomechanical model of cardinal vowel production: Muscle activations and the impact of gravity on tongue positioning[END_REF]) a simplified, five-parameter Mooney-Rivlin model is used for such passive tissues. The muscle fiber/soft tissue matrix interaction is modeled with a hyperelastic transversely isotropic material [START_REF] Weiss | Finite element implementation of incompressible, transversely isotropic hyperelasticity[END_REF]. In this matrix, a specific passive property is modeled along the direction of the muscle fibers, to take into account the passive properties of the fibers themselves [START_REF] Mcmahon | Muscles, Reflexes, and Locomotion[END_REF], while the properties of the surrounding passive tissues are modeled in the directions orthogonal to the fibers. The equation for the specific passive property along the fibers is taken from [START_REF] Blemker | A 3d model of muscle reveals the causes of nonuniform strains in the biceps brachii[END_REF]. Interaction between the muscle fibers and the soft tissue matrix is taken into account by considering the shear terms to strain energy as proposed by [START_REF] Criscione | Physically based strain invariant set for materials exhibiting transversely isotropic behavior[END_REF].

In the literature, the characteristic length l c , in the form of the c parameter of the original formulation of the λ model (c=1/l c ) varies between 9 mm [START_REF] Laboissière | The control of multi-muscle systems: human jaw and hyoid movements[END_REF] and 25 mm [START_REF] Buchaillard | A biomechanical model of cardinal vowel production: Muscle activations and the impact of gravity on tongue positioning[END_REF]. In our model, we used [START_REF] Buchaillard | A biomechanical model of cardinal vowel production: Muscle activations and the impact of gravity on tongue positioning[END_REF] value.

This framework was implemented as a user-defined, 3D element in ANSYS ® mechanical software.

Further details about the mathematical formulation of the DLM and its formulation can be found in [START_REF] Nazari | Modélisation biomécanique du visage: Etude du contrôle des gestes oro-faciaux en production de la parole[END_REF].

RESULTS

The results are organized into two main parts. In the first part, the ability of the DLM to account for the main features of the original λ model is assessed. These main features are:

1) The shape of the invariant characteristics: for the assessment of the DLM, the force-length relationships were measured, with and without the passive contribution, and they were compared to the original invariant characteristics.

2) The ability of the λ model to account for the activation of isometric and isotonic agonist and antagonist muscles: to assess the DLM, a model including agonist-antagonist muscle pairs was designed, with which agonist-antagonist co-activation without movement (isometric), and movement with a constant agonist-antagonist global force level were simulated.

In the second part, the DLM and a Hill-type muscle model were integrated into a sophisticated, 3D biomechanical model of the face [START_REF] Nazari | Simulation of dynamic orofacial movements using a constitutive law varying with muscle activation[END_REF]Nazari, Perrier, Chabanas & Payan, 2011). A preliminary comparison of these muscle models was carried out by simulating lip protrusion gestures, similar to those used in the production of rounded vowels, such as /u/ or /y/. This gesture was selected because it corresponds to a quite complex shaping of the lips in which mechanical soft tissue properties play an important role [START_REF] Nazari | Shaping by stiffening: A modelling study for lips[END_REF]. The characteristics of the simulated movements obtained with these two models were compared in terms of stress, strain, and energy.

Assessment of the DLM

Force-length curves: Comparison with original invariant characteristics.

To measure the force-length relationships of the DLM, the following method was used. A simple fixed-end bar (with an arbitrary resting length equal to 0.1m) was designed as a series of DLM muscle elements. These elements all generated a force along the same direction, namely the main axis of the bar (Figure 4). An external force was applied to the free extremity of the bar along the direction of the fibers; when all the muscle elements were activated according to a given λ value, the corresponding displacement of the free extremity was measured. Figure 4 illustrates an example of this. The left panel shows the bar at rest, and the right panel shows the lengthening of the bar under the influence of the external force.

---------------------(Figure 4 around here) ------------------------By varying the amplitude of the external force within a sufficiently large range it was possible to plot the force-length curve corresponding to the selected λ value. Different values of λ were considered in order to cover a sufficiently large range of muscle activations, and the force-length curve was determined for each of them.

For a given external force and a given λ value, the displacement was influenced by the passive properties of the DLM along the direction of the external force, which in this specific example was also the direction of the fibers. These passive properties were characterized by the global passive force-length relationship of the muscle along the fibers' direction, which is shown in Figure 5a (dotted curve). They resulted from the combined influences of the passive property of the tissues along the muscle fiber direction (Figure 5b, dotted curve) and of the passive force-length relationship of the surrounding tissues (Figure 5c, dotted curve). To be able to evaluate the respective importance of these passive properties and their combination, these passive force-length curves were depicted together with an example of theoretical invariant characteristic in the λ model. It can be seen that the passive influences are not negligible.

---------------------(Figure 5 around here) ------------------------

In Figure 6a an active force-length curve generated with the DLM without the effect of passive properties (dotted curve) was plotted for a given λ value, together with the corresponding theoretical invariant characteristic (equation [5]) in the λ model (solid line). This shows a good match between the numerical formulation and the original λ model. In Figure 6b the effect of total muscle force including the passive properties (Figure 5a) is shown for the DLM. It can be observed that with the inclusion of the passive properties the muscle starting length was slightly shifted to the right, which corresponds to an increase of the actual threshold length, λ * , as compared to the centrally specified control variable λ. This passive influence comes in combination with the effect of the proprioceptive feedback, the inter-muscular interaction, and the cutaneous feedback (Appendix 2, equation [8]).

The effect of the passive properties is shown for various values of λ in Figure 7. When the λ value increases above the bar length at rest (r=0.1m), the effect of the passive properties on active forcelength curve becomes more important.

---------------------(Figure 6 around here) ---------------------------------------------(Figure 7 around here) ------------------------

Isometric force variations and isotonic displacements

An important and well-known characteristic of the λ model is its capacity to generate, with an agonist-antagonist muscle pair, isometric force changes (i.e., a change in force in the absence of movement) or isotonic displacements (i.e., displacements without any change in the global amount of force produced in the agonist-antagonist pair). Isometric force variation is obtained with the coactivation of the two muscles, (i.e., the coordinated changes of the λ commands in the same direction). Isotonic displacement is obtained by a reciprocal change in the λ commands to the two muscles. In this section, the DLM is evaluated along these lines.

We aimed to build a model in which the interaction between muscles within an agonist-antagonist pair could be investigated. Two muscles with the same length at rest (0.1m) were attached together at one extremity, while their other ends were fixed (Figure 8). The λ value was set to 0.08m for one muscle (agonist) and to 0.09m for the other muscle (antagonist). Since the two muscles had identical properties, the theoretical invariant characteristics predicted a displacement of 0.005m in the direction of the muscle controlled by the smaller λ value. In our simulations, with the force-length relationship represented in Figure 6a (i.e., without passive influences), a displacement was generated, as expected, toward the agonist muscle, but the amplitude was equal to 0.0047m. This amplitude difference arose from the finite element discretization. Indeed, with an increased mesh density (five times larger), the displacement amplitude became 0.0049 m, which is closer to the expected value.

These results were different when the global passive properties (including fibers and surrounding tissues, Figure 5a) were taken into account (see Figure 6b). Indeed, the reached equilibrium position corresponded to a displacement of the attachment point equal to 0.0044 m for the low density mesh. This is in line with the observations presented in Figure 6b, since the increase of λ* due to the inclusion of passive properties induces a decrease of the level of force.

---------------------(Figure 8 around here) - -----------------------Once the model reached the equilibrium position, the λ values of both muscles were modified in the same way, to simulate agonist-antagonist co-activation. Table 1 shows the position of the attachment point associated with variations in λ values. The more the λ values decreased, the more the global force in the muscle pairs increased. The position of the attachment point stayed essentially constant, with a maximum variation of 0.0001m, which is around 2% of the global displacement from the position at rest. This corresponds well with an isometric agonist-antagonist co-activation.

Starting from the same initial equilibrium position (0.0044m from the position at rest) the λ value of the agonist muscle decreased, while the λ value of the antagonist muscle increased by an equal amount. This was done in 10 successive steps with a maximum λ shift of 0.005m. Table 2 summarizes the results of the corresponding simulations. Each change in the λ values generated a change in the equilibrium position, and the final distance of the attachment point with respect to the rest position was nearly twice as large as in the initial equilibrium position (2 × 0.0044=0.0088m).

During this displacement, the global force level remained essentially constant, with a maximum variation of 0.16N, which is less than 1%. This result agreed quite well with the characteristics of an isotonic displacement controlled by an agonist-antagonist muscle pair.

---------------------(Table 1 around here) ---------------------------------------------(Table 2 around here) ------------------------

In sum, the behavior of the DLM is in very good agreement with the theoretical behavior specified by the original λ model.

Application of the DLM in a biomechanical face model to produce lip protrusion

Hill-type models and the λ model differ in two basic features. First, the nature of the control variable is different. Hill-type models control the force directly, while the λ model controls the threshold muscle length that determines the relation between force and length. Second, the experimental foundation of the active force-length curve is different. In the Hill-type models it is based on measurements of isometric maximal voluntary force for various muscle lengths. In the λ model it is based on experimental measurements of force and length variations in unloading tasks for various levels of activation. Evaluating the consequences of the nature of the control variable is beyond the scope of the current study. In previous papers, we have shown why we consider the λ model to be useful in the context of speech motor control (Perrier et al., 1996a;[START_REF] Perrier | The Equilibrium Point Hypothesis and its application to speech motor control[END_REF][START_REF] Payan | Synthesis of v-v sequences with a 2d biomechanical tongue model controlled by the Equilibrium Point Hypothesis[END_REF][START_REF] Buchaillard | A biomechanical model of cardinal vowel production: Muscle activations and the impact of gravity on tongue positioning[END_REF]. Building realistic biomechanical muscle models will be part of our methodology to further address this issue in the future.

In the current study, we aimed to provide a first quantitative comparison of the mechanical behavior of the two muscle models. To do so, we implemented the DLM and a Hill-type model in a sophisticated biomechanical finite element model of the face [START_REF] Nazari | Simulation of dynamic orofacial movements using a constitutive law varying with muscle activation[END_REF] (Figure 9a). The implementation of the Hill-type model was based on Blemker et al.'s (2005) formulation. In line with our preceding work [START_REF] Nazari | Shaping by stiffening: A modelling study for lips[END_REF] we studied more specifically the consequences of the orbicularis oris peripheralis muscle (Figure 9b) on the lip protrusion gesture, which is crucial in the production of rounded vowels. This gesture is particularly interesting in the context of this comparison, since it involves a complex behavior of muscle tissue acting as a sphincter while generating a forward movement. It has been shown (e.g., [START_REF] Kim | Model-based investigation of control and dynamics in human articulatory motion[END_REF][START_REF] Nazari | Shaping by stiffening: A modelling study for lips[END_REF] that the correct achievement of this gesture is highly dependent on the stiffness properties of the muscle tissue.

---------------------(Figure 9 around here) ------------------------

The values of the control variables were carefully selected, in order for the face model to reach very similar final lip shapes for both muscle models. The obtained lip shapes were realistic and correspond well to the shapes that were experimentally observed in human subjects when lips were protruded.

The simulated lip shapes are plotted in Figure 10 (both muscle models provided the same shape).

---------------------(Figure 10 around here) ------------------------Our mechanical evaluation of the results was based on the computation of the 3D von Mises stresses and strains, which provide information about the global amount of stress and strain in the 3D soft body. The relation between these stresses and strains provides information about the required energy to shape the 3D soft body: for a given strain, the higher the von Mises stress, the higher the energy.

These mechanical variables were measured for all muscle activation levels in the final shape. The measurements were collected on three nodes of the face mesh, namely the lip corner, the central node of the upper lip, and the central node of the lower lip. These three points are classically considered to assess lip shaping [START_REF] Abry | Laws" for lips[END_REF].

The time patterns of stress and strain behavior and the stress-strain curves are shown in Figures 11 and 12 for the three nodes and for different levels of muscle activation. Some small differences can be observed between the muscle models. For the Hill-type model, the stress and strain tended to vary more quickly at the beginning of the movement and slow down at the end of the movement (Figures 11 and12, the two upper panels). This aspect was stronger for low activation (Figure 11) than for high activation (Figure 12). This is consistent with the fact that for the Hill-type model, the control variables act directly on the force level. Consequently these differences in timing do not seem to be primarily related to the specific biomechanical properties of the muscle models. Looking at the relation between stress and strain (lower panels), there was no clear difference between the two models. The curves were very similar, which suggests that the mechanical behaviors of the two models were very close, in terms of deformation of the lips and displacement of the selected nodes, as well as in terms of the energy required to generate the movement.

---------------------(Figure 11 around here) ---------------------------------------------(Figure 12 around here) ------------------------

DISCUSSION

Different muscle models have been proposed to study the complex control mechanisms of muscle mechanics by the CNS. The models all aimed to implement and/or explain the state of the muscle in response to CNS commands under the influence of external physical factors. In all the early models, for the sake of simplicity in calculations, a muscle was modeled as a one-dimensional (1D) actuator, generating active force, connected in parallel and in series with some springs and dashpots to simulate passive physical influences. With the development of fast computers and powerful numerical methods, it has become possible to design more sophisticated models with 3D details. All those general, 1D muscle models could be reformulated and fine-tuned to match the new demands of more realism in physical descriptions. In doing so, some parameters in 1D models would get new meaning.

At the same time, it would shed light on some details that may have been neglected due to oversimplification.

Currently, the 3D extension of Hill-type models is a well-developed subject in the field of biomechanics. This is not the case for the λ model, where the design of improved numerical descriptions of this model is required for its correct assessment in the context of debates about the Equilibrium Point Hypothesis. Our main goal in the present study was to fill this void and develop an extension of the λ model to meet these new needs. To do so, it was necessary to elaborate a new description of the model as a distributed version, which we called the Distributed Lambda Model or DLM. The DLM has given more physical meaning to some parameters of the λ model. It has helped to extract these parameters from a microscopic view and from experiments done at microscopic scales. It has also provided ideas for the design of potential new experiments to characterize these parameters. The main advantage of this modeling approach is the use of the λ model in powerful 3D numerical models such as the finite element method.

The DLM was shown to provide a very good account of the main properties of the λ model. In the absence of passive influences, the static force-length relations along the fiber direction, obtained for various values of the control variables, were in very good agreement with the original invariant characteristics proposed by [START_REF] Feldman | Once more on the equilibrium-point hypothesis (λ model) for motor control[END_REF]. It was possible to simulate, with a very good precision, isometric force changes and isotonic position changes. In addition, the DLM has allowed the quantitative study of the effect of muscle passive properties on the force-length characteristic along the fiber direction via their impact on the actual threshold length λ* (Figure 6b). When muscle contracts (i.e., when its length becomes smaller than its resting length), passive tissue properties cause a small increase in threshold length: the force generated for a given muscle length is somewhat smaller than the force predicted by the invariant characteristics. In muscle extension, when the length is larger than the resting length, this effect becomes rapidly prominent and it decreases the threshold length considerably. In the context of speech production, when muscles are active they generally contract. If they become elongated and active, their length stays quite close to their resting length.

Hence, in speech production the impact of passive tissue properties on force-length relationships appears to be limited. Therefore the approximation of the force-length relations with equations ( 5) and (8) (Appendix 2) appears to be accurate. For other motor tasks, the effect of passive tissue properties should be considered with care in the range of movements involving muscle length beyond resting length.

We also provided a first biomechanical quantitative comparison of Hill-type models and the λ model.

It suggested that in the range of voluntary movements used in speech production, in which muscles act mainly in contraction, there is not much difference between these two approaches in terms of stress, strain, and energy. Since we have shown that these differences are negligible, the use of the λ model in speech production appears to be perfectly justified from a biomechanical point of view as an alternative model to the Hill-type models. For motor tasks, in which muscles are likely to work in elongation, the Hill-type models give a better account of the trend for the muscle force to saturate above a certain level of elongation. Hence, combining the biomechanical properties of the λ model for contraction or small elongation with those of Hill-type models for significant elongation could be a powerful way to model muscle mechanical behaviors for a large range of movements. Proposals for Hill-type models and adjustable starting length muscle models along these lines have been made in the past by [START_REF] Winters | Hill-based muscle models: A systems engineering perspective[END_REF] and [START_REF] Shapiro | Control variables in mechanical muscle models: A minireview and a new model[END_REF].

CONCLUSION

We have elaborated a distributed version of Feldman's λ model of muscle mechanics and its control.

This model, the DLM, proved to behave in good agreement with the experimental data provided by [START_REF] Feldman | Once more on the equilibrium-point hypothesis (λ model) for motor control[END_REF], while realistically integrating the influence of the passive properties of muscles and their surrounding tissues. The integration of both a Hill-type model and the λ model in a 3D biomechanical face model provided a useful basis for studying the impact of muscle mechanics on speech facial gestures. Our preliminary results suggest that there are only few, small differences in normal speech articulatory movements, such as protrusion or rounding, between these two models.

The DLM opens the way to further speech production studies that associate simulation work and experimental studies, to provide quantitative evaluations of motor control hypotheses related to the nature of control variables (e.g., force, stiffness, or starting length), their specification from the definition of the motor task (e.g., with/without internal models), and their variations in time during the production of speech sequences.

APPENDIX 3

The sliding filament theory in the λ λ λ λ model

In order to account for the sliding filament theory [START_REF] Huxley | Muscle structure and theories of contraction[END_REF] the force described by the invariant characteristic curves was also scaled multiplicatively by [START_REF] Hill | The heat of shortening and the dynamic constants of muscle[END_REF] hyperbolic force-velocity term [START_REF] Laboissière | The control of multi-muscle systems: human jaw and hyoid movements[END_REF][START_REF] Payan | Synthesis of v-v sequences with a 2d biomechanical tongue model controlled by the Equilibrium Point Hypothesis[END_REF], St-Onge et al., 1997) whose values vary between 0 and 1 (see Equation 4). This term increases the damping characteristic of the model. The new equation below gives the final expression of the muscle force, including all the different contributions:

F=F passive +F active_Feldman *(f 1 +f 2 tan -1 (f 3 +f 4 v/l 0 )+f 5 v/l 0 ) (9)
where l 0 is the resting length i.e. the length at which the muscle can generate its maximum voluntary force, and f 1 to f 5 are constants used to fit the force-velocity characteristic of the muscle. Note that the velocity used in hyperbolic term is the current value of the velocity and not a delayed one as in Equation 5. Indeed, according to the sliding filament theory, this hyperbolic force-velocity term does not correspond to a feedback contribution, but to the direct impact of velocity on the mechanical properties of the actin-myosin bridges.

Tables

  muscle length change at time t, and µ is a damping coefficient. In this equation, both muscle length and muscle-length-change rate are delayed by a duration t d , to account for the physical delay in feedback propagation along afferent fibers.To take into account the passive mechanical properties of the muscle, a passive force should be added to the active force, similar to what is proposed for Hill-type models with the force F PE , in equation (1).Examples of schematic representations of the impact of the passive components are shown in Figure2. The invariant characteristic curves at zero velocity (static conditions) are plotted without (dotted curves) and with the inclusion (solid curves) of the passive component (represented by the dashed curve). Each force-length curve is fully specified by its activation parameter (i.e., by the threshold length λ, above which active muscle force is generated). If movement occurs due to external force, without any change in the control variable, the stretch-reflex mechanism will maintain the forcelength relation on one of the invariant characteristics.

Figure 3

 3 FIGURES CAPTIONS

Figure 4 .

 4 Figure 4. Experimental procedure for the measurement of the force-length relationship in the DLM: An external force is applied to a fixed-end muscle bar at rest and for a given λ λ λ λ value

Figure 5 .

 5 Figure 5. Force-length relations in the DLM: Influence of the passive tissues. For matter of

Figure

  Figure 6. Force-length relations in the DLM -Global assessment: (a) For same λ λ λ λ value,

Figure

  Figure 7. Global force-length curves in the DLM for various λ λ λ λ values (dotted) as compared to the corresponding ICs in the original λ λ λ λ model. The effect of passive properties tends to increase when λ λ λ λ increases, especially when λ λ λ λ becomes larger than the length at rest (0.1 m).

Figure 8 .Figure

 8 Figure 8. Modeling of an agonist-antagonist muscle pair.

Figure 10 .

 10 Figure 10. Example of the shape of the face reached under the activation of the OOP for both

Figure 11 .

 11 Figure 11. Comparison between a Hill-type model and the λ λ λ λ model for low levels of activation:

Figure 12 .FiguresFigure 1 .

 121 Figure 12. Comparison between a Hill-type model and the λ λ λ λ model for high levels of activation:

Figure 2 .

 2 Figure 1. Force-length relations in Hill-type muscle models: Changes in activation levels are

Figure 3 Figure 4 .Figure 5 .

 345 Figure 3 Comparison between the λ λ λ λ model and a Hill-type model: the bold solid path shows an

Figure 7 .

 7 Figure 7. Global force-length curves in the DLM for various λ λ λ λ values (dotted) as compared to the corresponding ICs in the original λ λ λ λ model. The effect of passive properties tends to increase when λ λ λ λ increases, especially when λ λ λ λ becomes larger than the length at rest (0.1 m).

Figure 8 .

 8 Figure 8. Modeling of an agonist-antagonist muscle pair.

Figure 10 .

 10 Figure 10. Example of the shape of the face reached under the activation of the OOP for both

Figure 11 .

 11 Figure 11. Comparison between a Hill-type model and the λ λ λ λ model for low levels of activation:

Figure 12 .

 12 Figure 12. Comparison between a Hill-type model and the λ λ λ λ model for high levels of activation:

Table 1 . Isometric agonist-antagonist co-activation in the model of Figure 8

 1

			Lambda Threshold	Lambda Threshold
	Displacement (m)	Force (N)	(Agonist) (m)	Antagonist (m)
	-0.0044	17.70	0.08	0.09
	-0.0044	18.47	0.0795	0.0895
	-0.0044	19.26	0.079	0.089
	-0.0044	20.06	0.0785	0.0885
	-0.0044	20.88	0.078	0.088
	-0.00444	21.71	0.0775	0.0875
	-0.0044	22.56	0.077	0.087
	-0.0045	23.43	0.0765	0.0865
	-0.0045	24.31	0.076	0.086
	-0.0045	25.20	0.0755	0.0855
	-0.0045	26.12	0.075	0.085

Table 2 . Isotonic agonist-antagonist displacement in the model of Figure 8

 2

			Lambda Threshold	Lambda Threshold
	Displacement (m)	Force (N)	(Agonist) (m)	Antagonist (m)
	-0.0044	17.70	0.08	0.09
	-0.0049	17.71	0.0795	0.0905
	-0.0053	17.73	0.079	0.091
	-0.0057	17.74	0.0785	0.0915
	-0.0062	17.76	0.078	0.092
	-0.0066	17.77	0.0775	0.0925
	-0.0070	17.79	0.077	0.093
	-0.0074	17.81	0.0765	0.0935
	-0.0079	17.82	0.076	0.094
	-0.0083	17.84	0.0755	0.0945
	-0.0088	17.86	0.075	0.095

The expression [G] + is equal to zero if the quantitity G is negative, and it is equal to G, otherwise.

Acknowledgement

This work has been partly funded by the French Agence Nationale de la Recherche (Project ANR-08-Blan-0272-01) and by the French German University (Project DFH/UFA G2R-FA-16-07).

We'd like to address special thanks to Matthieu Chabanas from Gipsa-lab, for fruitful discussions at early stages of the project.

APPENDIX 1 Adjustable stiffness muscle models

According to Equation 3, since the three functions, F v (v), F L (l) and f ac (A c ,t), are supposed to be independent, the stiffness (i.e., the length derivative of the force) is expressed as:

Thus, for given mechanical properties of the muscle (i.e., for given functions F v (v) and F L (l)), change in muscle activation directly results in change in stiffness. This multiplicative account of the contractile force is called adjustable-stiffness [START_REF] Shadmehr | A mathematical analysis of the force-stiffness characteristics of muscles in control of a single joint system[END_REF]: the contractile element is assumed to behave like a nonlinear spring whose stiffness varies with muscle activation, velocity and time. Under isometric condition, (i.e., when v is equal to zero; muscle length does not change during contraction), the stiffness varies linearly with f ac .

APPENDIX 2

Starting muscle length in the λ λ λ λ model

In successive versions of the λ model [START_REF] Feldman | Positional frames of reference in motor control: Origin and use[END_REF]St-Onge et al., 1997), inter-muscular interaction resulting in activation or inhibition of interneurons was taken into account functionally via a ρ parameter that modifies the λ value into a λ* value (see Equation 8). Very recently [START_REF] Pilon | Threshold position control of arm movement with anticipatory increase in grip force[END_REF] the cutaneous feedback, which acts in addition to the feedback arising from muscle spindles and/or Golgi tendons, was also integrated into the ρ parameter. [START_REF] Pilon | Threshold position control of arm movement with anticipatory increase in grip force[END_REF] also propose another shift (noted ε(t) in Equation 7) due to history-dependent changes in intrinsic properties of the motor neurons. Under the combination of these effects, the actual starting length ( λ ' ), which determines the force level for a given value of the centrally specified control variable λ, is given by: λ ' =λ-ρ+ε(t) (8)


Loïc Favergeon
email: favergeon@emse.fr

Jacques Morandini

Michèle Pijolat

Michel Soustelle

Kinetic Modeling Of Solid-Gas Reactions At Reactor Scale: A General Approach

A rigorous simulation of industrial reactors in the case of solid-gas reacting systems is a complicated task due to several difficulties linked to the kinetic problem at the scale of the solid grains and to the problem of gas and heat transfers within the powder bed. Firstly it requires the knowledge of the kinetic model for the calculation of the speed of reaction in one part of the reactor and for given conditions of temperature and gas composition, and secondly it necessitates solving the material and heat balance equations for the thermohydraulic conditions settled in the reactor.

Introduction

Solid-gas reactions are of great interest in many industrial fields such as nuclear, chemistry, metallurgy, CO 2 capture , etc… Industrial reactors where these reactions take place are difficult to understand. Indeed the solid phase is a granular medium through which circulate gaseous reactants and products. The properties of such a medium are modified in space and time due to reactions occurring at a microscopic scale. The thermodynamic conditions are driven not only by the operating conditions but also by the heat and mass transfers in the reactor.

Several models have been developed to account for the complexity of these transformations such as the grain model [1] and the pore model [2] and all their improved derivatives. However most of these models are based on the law of additive reaction times of Sohn [3] which is only valid in some special cases: if the order changes from the gas is equal to one, and for a shrinking core model. In all other cases, it is necessary to consider another approach to understand the kinetics of these reactions.

Method

In this goal, we have developed a multi-physic approach combining the resolution of the thermohydraulic equations with kinetic laws describing the heterogeneous reactions, for aggregated or separate particles [4].

Indeed, rather than the restrictive equation

) ( ) exp( α - = α f RT E A dt d (1)
and with about a tenth of different laws only, depending on the expression of the function f(α), we propose a more general approach for kinetic modeling at microscopic scale, including instantaneous nucleation models, instantaneous growth model, surface-nucleation and growth models (anisotropic or isotropic growth), for the three usual symmetries (spheres, cylinders, planes), and other possible kinetic assumptions depending on where the rate-determining step occurs, and the sense of development of the new phase. This approach is based on the following equations ( ) ( )

,... , t S P T dt d m i φ = α (2)
where φ(T,P i ) is the areic reactivity of growth and S m (t,…) is a molar space function. It allows to calculate about 45 various kinetic models (obviously including the classical ones as R3, D3, …) [5]. At macroscopic scale, heat and mass transfer terms entering in the balance equations depend on the kinetics evaluated at the microscopic scale. These equations give the temperature and partial pressure in the reactor, which in turn influence the microscopic kinetic behavior. The various equations corresponding to the thermal model, hydrodynamic and mass transport models can be found in [4].

At microscopic scale, the reaction fractional conversion is followed for a representative population of grains. By finite difference it is possible to evaluate the reaction rate dα/dt. This rate allows to calculate heat and mass sources produced by the reaction. So, using these sources terms, microscopic reactions have an impact on the spatial and temporal evolution of the thermodynamic processes at the reactor scale. Inversely since it modifies the kinetic parameters (areic frequency of nucleation γ and areic reactivity of growth φ), thermodynamics influences fractional conversion of the microscopic reaction.

Results and Discussion

In order to validate our approach we have chosen to study the dehydroxylation of kaolinite by means of thermogravimetry. Isothermal and isobaric curves were first obtained with about 25 mg of powder (corresponding to a 1 mm height powder bed). During the entire transformation, the sample weight was low enough to prevent any effect of pressure gradients in the powder layer. As expected, the reaction rate increases with temperature and decreases with water vapor pressure.

A second series of experiments was done with various heights of the powder bed. In the same crucible as previously, two distinct TGA experiments were performed at 450°C under a water vapor pressure of 7hPa: the first one with 50 mg of powder corresponding to a powder bed height of 3 mm; the second one with 300 mg corresponding to 10 mm. The kinetic curves and rates curves obtained with both bed heights are presented on Figure 2. One can see that the powder bed height has an important effect on the kinetic curves since the maximum reaction rate obtained for 3 mm is about twice more that obtained with 10 mm (and is about twice less than that obtained with 1 mm).

In order to interpret the experimental curves obtained with thin powder beds (1 mm height), we used a model involving random nucleation followed by anisotropic growth of the nuclei first proposed by Nahdi et al. [6].

The optimization procedure of CIN4 provided the values of γ (in nb. nuclei m -2 s -1 ) and φ (in mol m -2 s -1 ) for each experiment so the laws of variation of these quantities with temperature and water vapor pressure could be expressed.

Concerning the results obtained with thick powder beds, the "macro" part of CIN4 was used in order take into account the heat and mass transfers between all the various places inside the reactor (sample bed, crucible, gas flow). The first step consisted in defining the domains of calculation. The system being axi-symmetric, we used weak forms written in cylindrical coordinates.

Inside the granular medium, heat and mass sources were calculated simulating the reaction phenomena at the microscopic level. Figure 1 shows the gas flow rate obtained as well as temperature and fractional conversion fields at a given time for bed height of 3 mm. Figure 2 shows that the CIN4 simulation is able to discriminate the height of the powder bed from 3 to 10 mm. It is important to notice that these results at reactor scale are obtained by means of numerical simulation without any adjustment of parameters. Nevertheless the calculated curves are not perfectly superimposed with the experimental ones. Some improvements should be done in order to get better simulations. For example several points could be improved: at the grains scale ("micro" optimization), the particles size distribution could be taken into account to obtain more precise values of g and f, at the reactor scale, the change in the physical properties of the granular medium (including the porosity) due to the transformation of kaolinite into metakaolinite should be considered using appropriate relationships between the physical properties of the granular medium and fractional conversion. Fig. 2: Reaction rate vs. fractional conversion from the experiment and from the simulation for both powder bed heights.

Concluding remarks

CIN4 is a new software tool able to couple kinetics at grains population scale with both heat and mass transfers by means of the finite element method. Kinetic models are based on a general method agreeing with physical reality, offering the possibility to simulate more than 40 models depending on nucleation and growth processes, the shape of the grains, the growth type (isotropic or anisotropic), the development direction of the new phase and the localization of the rate determining step of growth (internal interface, external interface or diffusion volume). Heat and mass transfers are governed by differential equations and solved by a finite element method. Powder beds consisting in separate individual grains (micro and macro scales) or aggregates (micro-meso-macro scales) may be considered. CIN4 allows simulating not only reactions at laboratory scales such as thermobalances, but also industrial reactors.

Fig. 1 :

 1 Fig. 1: Representation of the temperature field, the fractional conversion and the flowing rate (arrows) at a given time for 3 mm-height bed.