Behavior of a net of fibers linked by viscous interactions: theory and mechanical properties

This paper presents an investigation of the macroscopic mechanical behavior of highly concentrated &ber suspensions for which the mechanical behavior is governed by local &ber–&ber interactions.Theproblem is approached by considering the case of a net of rigid &bers of uniform length, linked by viscous point interactions of power-law type. Those interactions may result in local forces and moments located at the contacting point between two &bers, and respectively power-law functions of the local linear and angular velocity at this point. Assuming the existence of an elementary representative volume which size is small compared to the size of the whole structure, the &ber net is regarded as a periodic assembly of identical cells. Macroscopic equilibrium and constitutive equations of the equivalent continuum are then obtained by the discrete and periodic media homogenization method, based on the use of asymptotic expansions. Depending on the order of magnitude of local translational viscosities and rotational viscosities, three types of the equivalent continua are proved to be possible. One of them leads to an e;ective Cosserat medium, the other ones being usual Cauchy media. Lastly, formulations that enable an e;ective computation of constitutive equations are detailed. They show that the equivalent continuum behaves like an anisotropic power-law <uid.


Introduction
Understanding the behavior of short ÿber-reinforced uids is of great interest in many industrial applications such as the processing of polymeric materials and other composites. In short ÿber systems processing, ÿbers noticeably in uence the ow of the uid matrix and conversely, the ow of the matrix determines the spatial distribution and orientation of ÿbers, which makes the modeling of such a problem complex.
Many relevant theoretical studies on this topic, forming an important part of the theory of suspensions, have tried to establish the relation between the microscopic properties of such materials, given by the behavior, geometry, orientation and distribution of ÿbers, and macroscopic mechanical behavior, under some restrictive assumptions. They generally apply to the case of rigid straight ÿbers immersed in a Newtonian uid when short range interactions between ÿbers may be neglected, that is to say for dilute and semi-dilute suspensions (Batchelor, 1971;Advani, 1994).
However, those theories cannot account properly for the e ects of local interactions between ÿbers, that may be due to dry friction e ects or localized viscous forces (Toll and Manson, 1994), and their validity is therefore limited to quite dilute suspensions, which is not the case of many industrial processes such as compression of Sheet Molding Compounds (SMC) or Glass Mat Transfer (GMT), or injection of Bulk Molding Compounds (BMC) (Dumont et al., 2003;Le Corre et al., 2002). Moreover, the complete solving of theoretical problems often requires further statistical assumptions that make the e ective calculation of the behavior possible only in the case of perfectly unidirectional or perfectly isotropic orientations of ÿbers (Fredrickson and Shaqfeh, 1989;Shaqfeh and Fredrickson, 1990).
Furthermore, almost no analytical solutions can be found in the literature in the case of non-Newtonian uids reinforced with ÿbers. This is due to the di culty of calculating velocity ÿelds around ÿbers as in the Newtonian case (Je ery, 1922), and to the non-applicability of the superposition principle often used in that case. However, according to the work of Batchelor (1971), the behavior at high concentration regimes is mainly dominated by short range interactions between ÿbers, the contribution to the total stress of the matrix and of the ÿber-matrix interaction becoming then negligible. The critical ÿber concentration for such an assumption to be valid cannot be established in a general way. It largely depends on the ÿber's geometry but it is also conditioned by the nature of the matrix and of the interactions.
Assuming a highly concentrated regime, the macroscopic behavior of a suspension can be drawn from some simple micro-mechanical considerations as it is done in the works of Toll and Manson (1994), Gibson and Toll (1999) and Servais et al. (1999) in the case of planar ÿbers linked by a combination of dry and non-linear viscous interactions. However, the interesting results obtained by these authors lack generality for they only apply to simple viscometric ows such as biaxial extension or simple shear ow.
As a ÿrst step towards a more general approach to the modelling of the behavior of short ÿber systems, this paper exposes an homogenization method suitable for the modelling of highly concentrated ÿber suspensions linked by non-linear viscous interactions of power law type. This method is an application of the homogenization method of discrete and periodic media, initially developed by Moreau and Caillerie (1995), Tollenaere and Caillerie (1998) and Pradel (1998) for the modelling of periodic trusses or foams in the scope of elasticity. It is based on the use of asymptotic expansion methods for periodic homogenization proposed by Bakhvalov and Panasenko (1989), Bensoussan et al. (1978) and Sanchez-Palencia (1980), adapted to discrete problems. Starting from a discrete problem at the scale of the ÿbers (microscopic scale), the proposed method enables ÿnding the essential properties of the equivalent continuum, that is to say the general form of its balance and constitutive equations at the macroscopic scale.
In Section 2 are detailed the basic assumptions of the upscaling method, the notations adapted to the discrete geometrical description of the net and the assumptions relative to the modelling of interactions between ÿbers. In Section 3, the upscaling process, based on asymptotic expansions, is discussed and preliminary results are exposed. This process then enables the discrete balance equations of the ÿber net to be transformed into continuous ones: this leads to the deÿnition of the macroscopic stress tensors of the equivalent continuum (see Section 4). Depending on the ÿber-ÿber interaction laws at the microscopic level, three types of macroscopic constitutive equations are obtained, the e ective computation of which are then discussed (see Section 5). Finally, in Section 6, simple considerations on those constitutive equations show that in some cases, the equivalent continuum may be modeled by an anisotropic power law uid.
For the sake of simplicity, the problem is exposed in the case of a planar ÿber net, but, as it will be clear in the equations, the extension to a three dimensional problem would be straightforward.

Notations
Boldface symbols denote tensors, the order of which is indicated by the corresponding number of underlinings. Dots and colons are used to indicate tensor products contracted over one and two indices respectively. In the usual Cartesian frame, this leads to using Einstein's summation convention over repeated indices.
Tensorial product is denoted by the symbol ⊗, e.g.: The gradient of a vector a with respect to space variables x i will be denoted ∇ a and deÿned as The same convention will be used for higher order spatial derivatives, e.g.: (∇ a) ijk = 9 2 a i 9x j 9x k :

Basic assumptions
The structure under consideration is a planar net made of rigid cylindrical ÿbers of uniform length l, that will be called bars. The plane is characterized by the Cartesian reference frame R = (O; e 1 ; e 2 ), where O is an arbitrary reference point (see Fig. 1).
The discrete homogenization technique we propose requires two further assumptions. The ÿrst assumption is that the structure may be considered as periodic. Thus the ÿber net is regarded as an assembly of identical cells as shown in Fig. 1, which are characterized by the two vectors of R , d 1 and d 2 , called periodicity vectors. The ÿber net dimensions are therefore L 1 and L 2 which are such that where N 1 is the number of cells in the d 1 direction and N 2 the number of cells in the d 2 direction. The second assumption is that the ÿber net is made of a huge number of cells N c , so that the scale separation parameter , deÿned as may be considered as a very small parameter. Condition (2) is equivalent to the following, and so can have another meaning: where s and S are respectively the surface of one cell and the surface of the whole ÿber net, d = √ s represents the characteristic length of the microstructure and L = √ S is the characteristic length of the whole net. Let us assume that d has the same order of magnitude as the length of ÿbers. Condition (3) is therefore a condition of scale separation; it implies that the size of the ÿber net should be large in comparison with the size of ÿbers.

Numbering system
The net of ÿbers is made of a discrete set of rigid bars. Each bar is supposed to be linked to the rest of the net by one or several interactions, located at the contacting point between two ÿbers. The geometry and topology of the net is therefore entirely deÿned by the spatial and angular position of each ÿber, and by the connectivity of each ÿber with the rest of the assembly.
Centers of cells are ÿrst located by a vector of integers a whose components are (a 1 ; a 2 ). Their positions are located by vectors p 0 (a) (see Fig. 1): p 0 (a) = a 1 d 1 + a 2 d 2 : (4) The assumption of periodicity then suggests a system of numbering of bars and links re ecting the regularity of the microstructure. Each bar of the net is numbered bỹ b = (b; ab), which means that the barb is the b th bar of the cell ab. The set all bars of the net is denoted B.
In the same way, the action of the barc on the barb is denoted either byk = (c=b) or byk = (k; ak ), the set of connections of the net being denoted C. Barsb andc can respectively be considered as the interior and the exterior of the interactionk, so the following notation will be used: b = I(k);c = E(k) and sok = (E(k)=I(k)): We will consider that the cell to whichk belongs is the belonging cell of the bar on which the action is exerted, i.e. the barb = I(k), so that ak = ab. The reciprocal interaction (b=c) will be denoted tk = ( t k; a tk ) and will therefore belong to the cell ac.
In the notationsb = (b; ab) andk = (k; ak ), integers b and k are the numbers of bars and connections in a reference cell the sets of which are denoted B R and C R . In this work, we will assume that the size of cells d is about the same order as the length of ÿbers l, and that d ¿ l. This implies that interactions can take place only between two ÿbers of the same cell or between ÿbers of two neighboring cells. The exterior of the connectionk = (k; ak ) will therefore be located in the cell ak or in a neighboring one. Anyway, it will belong to the cell ak + k , where k is a vector of integers whose components take their values in {−1; 0; +1}. Notice that the periodicity assumption causes k to be independent of the position of the cell.

Geometrical description
As illustrated in Fig. 2, the position of a barb is deÿned by the position of its center Pb, located by the vector p(b) = OPb, and by its unit vector eb. The periodicity of the net implies that p(b) can be partitioned as where p 0 (a) is the macroscopic position of the belonging cell ofb in the ÿber net, and p b 1 is the local position of bar b in that cell. It is to be noted that p b 1 only depends on the considered bar b and not on the macroscopic position for the net is assumed to be perfectly periodic. Some extension to this restriction could be achieved by considering a quasi-periodic ÿber net, as done by Tollenaere and Caillerie (1998), but is not in the scope of this paper, which is only interested in obtaining the equivalent continuum's constitutive equations.
Vectors d i are small compared to the size L of the net, so they can be written as d i = i where the vectors i are macroscopic vectors, independent of . Let us now introduce the pseudo continuous variable b deÿned by b = ( 1 ; 2 ) = ab. According to this notation, Eq. (6) becomes p(b) = 1 1 + 2 2 + p b 1 = p 0 ( ) + p b 1 : As illustrated in Fig. 3, macroscopic positions of bars are now parameterized by the vector of reals which is a variable of ⊂ [0; 1]×[0; 1]. is the reference parametric space, where the net is made of square cells of size . As is assumed to tend to zero in the use of asymptotic expansion methods, tends to become, and will be used as a continuous variable, even if, strictly speaking, it only takes values such as = a for a ÿnite value of . Let us introduce G , the gradient ∇ p 0 of the geometrical transformation → p 0 ( ). According to Eq. (7), the components of G relative to e i ⊗ e j can be expressed as and the Jacobian g of the transformation is

Local kinematics and interactions
The problem being planar in R and ÿbers being considered as rigid bars, their motion is a rigid body planar motion. Kinematics of a barb can therefore be deÿned by: • the velocity of its center Pb, denoted C(b) • its angular velocity !(b) = !(b)e 3 , where e 3 is the unit vector normal to the plane of ÿbers. The velocity of a point Q of barb is then where is the curvilinear abscissa of Q on the barb with respect to the center Pb as shown in Fig. 4. is a small variable; its order of magnitude is about the size of the cell d so it can be written as = d = L , where is independent of the size of the problem L. In the following, it will be convenient to introduce the new variable = L!, which enables us to rewrite Eq. (10) as The actionk of bar E(k) on bar I(k) is supposed to take place at a point and to be of non-linear viscous type, due to the relative velocities of both bars at their contacting point. In this paper, we will consider the case where those interactions follow a power law of the relative velocities. The interactionk can therefore be partitioned into: • a force fk proportional to the di erence of velocity of the two bars Ck at their contacting point, denoted Qk : where the power law index m is a real scalar ranging from 0 to 1. By deÿnition, Ck , so making use of Eq. (11), its expression reads In the last equation, k (resp. t k ) is the normalized local abscissa of Q k on the bar b = I(k) (resp.c = E(k)).
• a moment Mk (Qk ) relative to Qk , proportional to the di erence of angular velocities of both bars: where !k denotes the di erence of angular velocities !(E(k)) − !(I(k)).
In order to normalize moments with respect to forces, vectors mk (Qk ), such as Mk (Qk )= Lmk (Qk ), will be used. In the following, they will also be called moments, even if they are homogeneous to forces. The moment interaction law can then be rewritten as where ÿ k = B k =L m+1 . The moment of actionk expressed at point Pk , the center of bar I(k), will be denoted mk and will have the following expression: In Eqs. (12) and (15), the coe cients k and ÿ k are respectively called translational and rotational viscosity. According to the assumption of periodicity, they only depend on the considered interaction k inside the reference cell. Those viscosities are chosen such that and where k 0 and ÿ k 0 are strictly positive real scalar constants of same order of magnitude and where the parameter q is a real number, characteristic of the relative order of magnitude of ÿ k with respect to k . In the following, only cases where q ¿ 0 will be examined, for situations where ÿ k would be much greater than k are not physically likely to happen in the considered suspensions. For simplicity, three cases will be discussed. They may represent three di erent situations one might expect when studying a speciÿc application: • Case 1: q=0, rotational viscosities have the same order of magnitude as translational ones. • Case 2: q = 1, rotational viscosities are "quite" small compared to translational ones. • Case 3: q = 1 + m, rotational viscosities are very small or negligible.
As it will be clear in the following (see Section 5.3), cases where q = 1 + m + j, j being a strictly positive integer are equivalent to Case 3 with ÿ k 0 = 0, so the results that will be drawn in Case 3 will also apply to such physical situations.
The analysis of intermediate cases, where q is a real positive number, was also carried out. They were shown to lead to three types of macroscopic descriptions (i.e. equivalent continua) which are identical to the descriptions that can be deduced from the analysis of the three proposed cases. The reader should therefore keep in mind that results exposed in the following sections are not restrictive.

Upscaling process
At this stage, it is not possible to give an explicit formulation of forces and moments in terms of the macroscopic velocity gradient of the suspension, as it is generally done in the classical theory of ÿber suspensions (Batchelor, 1971;Gibson and Toll, 1999). Such a process would require a further assumption on the velocity ÿelds relative to each bar which is not required here. As it will be clear in the following, the upscaling process used in our homogenization method will provide local forces and moments as implicit functions of a macroscopic velocity gradient which might be assimilated to the bulk velocity gradient of the suspension, as done in the previously cited works. But this property will be a result of the homogenization process and not an assumption.

Velocity ÿelds
Taking advantage of the smallness of the parameter and of the assumption of periodicity, velocity and angular velocity ÿelds relative to a barb are expanded in discrete asymptotic series expansion in powers of (Moreau and Caillerie, 1995;Tollenaere and Caillerie, 1998). Thus, they are written as where functions C b n and b n are continuous functions of that generally depend on the bar b. re ects the macroscopic variation of those functions whereas index b is a local variable, depending on the corresponding bar ofb in the reference cell.
For a pair of connected barsk = (c;b), let us recall thatb is supposed to belong to the cell located by , whereasc may belong to a neighboring cell located by + . Thus C(c) = C E(k) ( + k ) is expanded using a Taylor expansion around , using the fact that k is very small with respect to : Carrying asymptotic expansion (19) into Eq. (21) and identifying each order of powers of , the expansion of C(c) reads The expansion of the angular velocity (c) is obtained exactly in the same manner.

Velocity di erence ÿelds
In the same way, velocity di erences Ck (13) and k = E(k) − I(k) can be expanded: By substituting expansion (19) into Eq. (13), and making use of Taylor expansions (22), velocity di erences C k i can be identiÿed at each order of powers of . For the two ÿrst orders ( 0 and 1 ), we obtain For the two ÿrst orders of angular velocities di erence, one obtains

Equilibrium of the net
In the scope of this paper, inertia e ects will be neglected and no external forces or torques will be supposed to act on bars of the net. Therefore, the equilibrium equations of a given barb can be described by the two following sets of equations: where C(b) is the set of connections of barb. The expansion of equilibrium equations of the whole ÿber net can be simpliÿed by the use of a virtual power formulation. This technique, developed by Moreau and Caillerie (1995) and Tollenaere and Caillerie (1998) enables the macroscopic behavior to be obtained fast and avoids the need to consider expansions of the equilibrium equations at higher orders of powers of . Thus using two sets of virtual ÿelds denoted ( and '), problem (28) is equivalent to the following virtual formulation: Within such a summation, it can be noticed that if the connectionk = (c=b) belongs to C(b), connection tk = (b=c) belongs to C(c). Summation on the bars can therefore be transformed into a summation on the set of connected pairs (c=b) such thatc ¿b. Doing so, formulation (29) can rewritten as This last equation was obtained using the action and reaction theorem which in this case can be shown to imply where vector L k denotes the vector PbPc and is therefore deÿned as

From discrete to continuous
Suppose that ub is any vectorial ÿeld of the spatial variable with the following two properties: Then it can be shown that when tends to zero, we have using the deÿnition of integrals by Riemann sums as detailed in Moreau (1996). The latter property is important and will often be used in the following developments. It enables us to transform the discrete problem into a continuous one.

Self-equilibrium at the lower orders
In a ÿrst stage let us assume that both order zero velocities C b 0 and angular velocities b 0 actually depend on the considered bar, so that order zero velocity di erences C k 0 and k 0 are non null functions.
3.4.1. Self-equilibrium of forces According to asymptotic expansions (23), expressions such as Ck m−1 have an asymptotic expansion of the following form: Subsequently, at this stage, the asymptotic expansions of interaction forces, deÿned by Eq. (12), can be written as Their asymptotic expansion may therefore be written in the following way: with, at the lower order of powers, Taking into the virtual power formulation (30) virtual velocity ÿelds (b) and '(b) such that ∀b; then carrying into Eq. (30) asymptotic expansions (36) and making use of property (33), we obtain, when tends to 0: This relation being veriÿed for any continuous function 1 , forces of the lower order are solution of the lower order self-equilibrium equation of the reference cell: Then carrying into this last equation the constitutive relationship (37), we obtain a relation that must be satisÿed by the ÿrst order velocities of bars in the reference cell. Velocities C b 0 are actually solutions of the problem In the special case where b = C b 0 , one checks that order zero velocities have to satisfy the condition Translational viscosities being strictly positive quantities, Eq. (42) means that all the bars of the same cell have the same linear velocity at the macroscopic scale: According to such a property, we also deduce that Such a result is however inconsistent with the asymptotic expansion (36) proposed for interaction forces which was based on the assumption that Ck 0 = 0. Actually, according to property (43), expressions such as Ck m−1 have an asymptotic expansion of the form: Subsequently, asymptotic expansions of interaction forces can be written as with, at order 0 ,

Self-equilibrium of moments
In the assumption of non null order zero angular velocities, by substituting (46) and expansion (23) into the interaction law (16), one obtains the following expansion of interaction moments: Let us now distinguish the three cases mentioned in Section 2.4. Cases 1 and 2 can be treated in the way as was done for forces. Taking into the virtual power formulation of equilibrium (30) virtual velocity ÿelds (b) and '(b) such that ∀b; and then carrying into Eq. (30) asymptotic expansions (49), making use again of property (33), we obtain, when tends to 0 According to property (33), all the higher order terms vanish when tends to zero. The last relation being veriÿed for any function 2 , order zero angular velocities are then shown to solve the following self-equilibrium problem on the reference cell: This problem is strictly equivalent to the problem satisÿed by order zero velocities so, for cases 1 and 2, we also obtain ∀b ∈ B R ; b 0 ( ) = 0 : No such simple conclusion can be drawn in case 3, and order zero angular velocities should actually depend on the considered bar b in the reference cell. However, asymptotic expansions of moments are now expressed as described in Eq. (49), but the exponent 1 + q − m now equals 2. For the three cases, asymptotic expansions of moments thus can be written as with the following expressions of ÿrst order terms: • Case 1: • Case 2: • Case 3:

Equilibrium of the equivalent continuum
Before detailing the obtaining of the equivalent continuum of the ÿber net, let us recall some fundamental results from the mechanics of continuous media with internal rotation. The equivalent continuous description deduced from the discrete homogenization process will actually prove to be consistent with such general results.

Continuous media with internal rotation
A medium for which microscopic kinematics may imply rotational motion such as granular materials, foams or ÿber suspensions, may be considered as a Cosserat medium (Cosserat and Cosserat, 1909), or more generally as a continuous medium with internal rotations. Local kinematics of such a medium is characterized by a velocity ÿeld C(x) and an angular velocity ÿeld (x) and, as described by Eringen (1968) and Germain (1973), the virtual power formulation of its equilibrium reads ∀(C * ; * ); − ( : ∇ C * + · * + Ä : ∇ * ) d where C * and * are virtual velocity and angular velocity ÿelds. In Eq. (61), represents the stress tensor, Ä the couple-stress tensor (due to local rotations) and the micro-stresses vector. Vectors f and c respectively denote the external volume forces and moments. The principle of objectivity causes to be directly linked to by the relation: tensor e being the permutation tensor whose deÿnition is given in Appendix A and A is the antisymmetric part of . Then by splitting ∇ C * into its symmetric and antisymmetric parts D * and R * and accounting for relation (62), formulation (61) may also be rewritten as Such a formulation leads to the following local equilibrium equations: and implies, for a viscous material, constitutive equations of the generic type: In the special case where Ä and external moments c actually happen to be null or negligible, ∇ no longer in uences the behavior. It is then clear from Eqs. (64) that vector e · A is also null or negligible. In such a case, constitutive equations (65) simplify to The symmetry condition (67)

Equilibrium of the ÿber net
For the ÿber net, assuming that external forces and moments are null, the virtual power formulation of the equivalent continuum is obtained by taking in formulation (30) smooth "macroscopic" virtual velocity ÿelds and ' such that ∀b; (b) = 0 ( ) and '(b) = ' 0 ( ): Then when tends to zero, using Eq. (33), forces and moments can be shown to be solutions of the problem: which, using the property (a · b) · c = b ⊗ c : a, can also be written as In this last formulation, the following stress tensors, deÿned in the parametric space , were introduced: The formulation of the macroscopic equilibrium of the net in terms of the velocity gradients in the physical space are simply obtained by making the change of variables → p 0 ( ) in formulation (71). According to deÿnition (8), the relation between the gradient of a vectorial ÿeld u with respect to , denoted ∇ u, and its gradient with respect to x = p 0 , denoted ∇ x u, is given by The macroscopic equilibrium of the ÿber net in the physical space then reads This last problem is the virtual power formulation of the equilibrium of a Cosserat continuous medium without any external forces or moments, as detailed in Section 4.1. Its local equilibrium is governed by the following balance equations: Therefore, the state of stresses inside the ÿber net is deÿned by the three following tensors: 0 = g −1 G · S 0 : stress tensor; Ä 0 = g −1 G · M 0 : couple stress tensor; 0 = g −1 Z 0 : micro-stresses vector: In accordance with the general theory of Cosserat media exposed in Section 4.1, the macroscopic stress tensor 0 is non symmetric, and its antisymmetric part is directly linked to 0 by the relation 0 = e : 0 = e : A 0 : This general property of Cosserat media can easily be checked in the case of the ÿber net problem, as shown in Appendix A.
As visible in Eqs. (78) and (72)- (74), the state of stress of the equivalent continuum happens to be directly related to ÿrst order interaction forces and moments. The sti ness of the ÿber net will therefore by closely linked to its density of connections. This remark is consistent with the results obtained by Servais et al. (1999) in the case where dry friction between ÿbers may be neglected as well as local moments.

Constitutive equations of the equivalent continuum
As Eqs. (72)-(74) show, the determination of the macroscopic state of stress ( 0 ; Ä 0 ; 0 ) in the ÿber net requires the determination of forces and moments of order 0 , f k 0 and m k 0 . Nevertheless, with the chosen local interaction relations, those quantities directly depend on the value of kinematic variables C b 1 , b 0 and b 1 , as well as on the macroscopic velocity gradient ∇ C 0 . To fully determine the constitutive equations corresponding to the three cases deÿned in Section 2.4, further equilibrium formulations are required in order to enable the determination of local kinematic variables.
Such formulations will necessarily depend on the local interaction laws. Here again, one has to distinguish three cases, depending on the relative magnitude of rotational viscosity ÿ k with respect to translational one k , characterized by the parameter q (see Section 2.4).

Case 1: q = 0
According to results (43) and (53), order zero interaction laws become Carrying into the general equilibrium formulation (30) the virtual functions and making tend to zero, one gets This relation being satisÿed for all ÿelds 1 and 2 , forces f k 0 and moments m k 0 are solutions of the problems: Formulations (86) and (87) are respectively force and moment order zero self-equilibrium equations. They are strictly equivalent to the following non-linear systems: Constitutive equations of the equivalent continuum can then be calculated considering as given the macroscopic ÿelds ∇ C 0 , ∇ 0 and 0 . The computation gives velocities C b 1 and b 1 as functions of ∇ C 0 , ∇ 0 and 0 , then Eqs. (72)- (74) and (78) give 0 , Ä 0 and 0 in terms of those ÿelds, which provides constitutive laws. Such a computation obviously requires a numerical implementation of problems (88) and (89). This will explicitly provide the macroscopic stress tensors as functions of the macroscopic ÿelds and can be achieved by the use of any suitable numerical methods for non-linear systems. Subsequently, in case 1, according to local behavior equations (80) -(83), the equivalent continuum is a general Cosserat medium whose constitutive relationships are of the following type: which is consistent with the general formulation (65) obtained by continuous media theory (see Section 4.1).

Case 2: q = 1
In this case, as shown in Section 3.4 order zero moments are null. The macroscopic couple-stress tensor Ä 0 is therefore automatically null, and from Eqs. (77), the micro-stresses vector 0 is also null. Property (79) then immediately causes the antisymmetric part of the macroscopic stress tensor to be null. Thus, the equilibrium of the equivalent continuum does not imply any local moment and 0 is a symmetric tensor.
In this case, the only constitutive equation to determine is therefore Eq. (72), which only requires the determination of forces f k 0 . Thanks to property (53), their expression is the same as in case 1; they are deÿned by Eq. (80), with C k 1 given by Eq. (82). Then, adopting the same technique as in case 1, order zero forces f k 0 can be proved to solve the self-equilibrium (86), which is strictly equivalent to the following non-linear system: Furthermore, property (79) implies an additional relation on forces f k 0 which reads This last relation shows that 0 can theoretically be expressed as a function of ∇ C 0 . One then notices that Eqs. (91) and (92) enable C b 1 and 0 to be calculated in terms of the macroscopic ÿeld ∇ C 0 . Here again, an explicit determination of the macroscopic stress tensor 0 can be achieved by the numerical solution of the non-linear system formed by both equations.
Finally, in this case, the ÿber net's equivalent continuous medium happens to be analogous to the special case of the continuous medium discussed in Section 4.1, governed by the classical local equilibrium equation 20 with a constitutive equation of the type 0 = 0 (∇ C 0 ; 0 (∇ C 0 )); where 0 is a symmetric tensor. The relation between 0 and ∇ C 0 , is obtained in an implicit way by ensuring the symmetry condition. Furthermore, in accordance with the general theory, because of the symmetry of 0 , only the symmetric part of the macroscopic velocity gradient ∇C 0 contributes to the total dissipated mechanical power. This causes 0 to depend only on the macroscopic strain rate D 0 , deÿned as So in Eq. (94) ∇ C 0 may be replaced by D 0 and the equivalent continuum exhibits a general uid-like behavior.

Case 3: q = 1 + m
In this case, moments of order 0 are immediately null, which causes Ä 0 and 0 to be null. As in case 2, the equivalent continuum is a Cauchy medium, governed by local balance equation (93), and the stress state is deÿned by the single tensor 0 . Its determination, according to Eq. (72), requires the determination of order zero interaction forces f k 0 . Here again, the set of self-equilibrium formulations (86) and (87) can be obtained by the same process as in the two previous cases, but Eq. (87) no longer brings further information on order zero angular velocities b 0 . Those variables now depend on the considered bar in the reference cell.
Forces f k 0 are therefore deÿned by They solve Eq. (86), which is equivalent to the non-linear-system: the unknowns of which are now C b 1 and b 0 . It therefore brings 2 equations per bar whereas 3 unknowns per bar are to be determined, ∇ C 0 being considered as data.
It is to be noted that in this case, order one moments m k 1 , deÿned by imply the same kinematic variables as f k 0 . The missing equations thus can be obtained by taking in Eq. (30) the virtual functions (b)=0 and '(b)= 2 ( )' b leading to the moments order one self-equilibrium following formulation: Making use of the action-reaction theorem (31), this equation can be transformed and proved to be equivalent to a new non-linear system, which reads Vectors C b 1 , and b 0 can therefore be computed by the simultaneous solving of Eqs. (98) and (101) in terms of ∇ C 0 , what will provide f k 0 , and then 0 in terms of this macroscopic velocity gradient.
Finally, as in case 2, the ÿber net's equivalent continuous medium is also a classical continuous medium governed by the local equilibrium equation (93), and its constitutive equations are of the type: where 0 is a symmetric tensor and D 0 is the macroscopic strain rate tensor. Cases 2 and 3 ÿnally happen to lead to the same type of equivalent continuous medium, even if the calculation of their constitutive equations leads to somewhat di erent resolution schemes. As mentioned above, it is now clear that case 3 also includes cases where q = 1 + m + j, where j is any positive integer. In such cases, rotational viscosities vanish from the macroscopic constitutive equations, and only the determination of the translation one is required for a full solution of the problem.

Fundamental properties
In the case of the power law interaction relations (12) and (16) discussed in this paper, further properties of constitutive equations can be drawn from results exposed in the previous sections.
In a ÿrst stage, let us focus on case 3, where rotational viscosities are assumed to be very small compared to translational ones. Self-equilibrium equations (98) and (101) form a system that enables the calculation of local variables C b 1 and b 0 in terms of ∇ C 0 . Such a system can be summarized in the following way: where F is a vector that contains the force and moment equilibrium equations of each bar of the reference cell, and X is a vector containing the kinematic unknowns C b 1 and b 0 relative to each bar. F is a block vector where blocks of components [3b − 2 : 3b], relative to bar b, are denoted F b and read 22 whereas X may be assembled in blocks X b such as In the following, the uniqueness of the solution of problem (103) will be assumed. Nevertheless, when internal mechanisms exist (isolated bars or isolated groups of bars) this will not be the case anymore, but we will assume that the concentration regime is su ciently high for every bar to be connected with some of its neighbors.
Let us now consider that X is the solution of Eq. (103) and thatX is the solution of F(X; ∇ C 0 ) = 0; where is a non-null real scalar, and study the relation between X andX.
0 are the respective solutions of problems (103) and (106), forces and moments solutions of Eq. (106) are then They can therefore be rewritten aŝ (111) and (112) in Eq. (106) and according to the expression of vector F, one checks thatX is solution of the problem: The uniqueness of the solution then immediately causeŝ Subsequently, forces and moments resulting from problem (106) are such that where f k 0 and m k 1 result from the initial problem (103). Now forming the macroscopic stress tensorsˆ 0 and 0 deÿned by Eqs. (78) and (72) corresponding to both problems ÿnally leads tô 0 = m 0 : This result therefore enables one to deduce the following important property: The macroscopic stress tensor is a homogeneous function of degree m of the macroscopic strain rate tensor. Such a result shows that, in case 3, the equivalent continuum is a power law uid, with a strain rate sensitivity equal to the strain rate sensitivity postulated at the level of interactions between ÿbers. Actually, if one deÿnes a norm : eq in the space of second order tensors, property (117) enables us to write the following relation, characteristic of power law uids: The same property can be deduced for case 2, thanks to interaction relation (80) and admitting the uniqueness of the system formed by Eqs. (91) and (92).
In case 1, macroscopic ÿelds ∇ C 0 and 0 (directly linked to ∇ 0 ) can be imposed separately so no such simple property can be deduced. However, as evident from the interaction relations (80) and (81), and from the formulation of self-equilibrium problems (88) and (89), the following property can be written: The equivalent Cosserat medium exhibits a degree m homogeneity property in terms of the pairs (∇ C 0 ; 0 ).

Conclusions
This theoretical work on the behavior of a net of rigid ÿbers linked by punctual power law ÿber-ÿber interactions shows several interesting results.
If the scale separation assumption (3) is satisÿed in any practical application, an equivalent continuous description of the behavior of the net is possible and its general equilibrium equations are typical of a Cosserat continuous medium. The state of stress of this medium is entirely deÿned by Eqs. (72)-(74) that explicitly provide the link between the local forces and moments and the macroscopic stress tensors.
Furthermore, the analysis of three di erent ÿber-ÿber interaction laws leads to 2 main di erent types of equivalent continuous media depending on the relative order of magnitude of rotational viscosities with respect to the translational ones. If rotational viscosities are of the same order of magnitude as translational ones, the ÿber net is actually equivalent to a Cosserat medium, its state of stress is given by the usual, but non-symmetric, Cauchy stress tensor and by a couple stress tensor accounting for local moments generated at ÿber-ÿber interactions. Such a case would probably be relevant for almost rigid interactions (m ≈ 0). If rotational viscosities become smaller, that is to say in cases like cases 2 or 3, the equivalent continuum is a usual Cauchy medium, deÿned by a single symmetric stress tensor.
Constitutive equations of the ÿber net's equivalent continuous medium cannot be obtained in an explicit form. They require the numerical determination of each order zero forces and sometimes of order zero or order one moments. Nevertheless, such a computation can be achieved quite simply and does not imply huge numerical problems, thanks to the periodicity assumption.
In a last stage, another fundamental property of the equivalent continuum was drawn. Thanks to the power law nature of the ÿber-ÿber local interaction laws, in cases 2 and 3, the macroscopic stress tensor could be proved to be a degree m homogeneous function of the macroscopic strain rate tensor. This shows that the equivalent continuum is a power law and anisotropic uid with the same strain rate sensitivity m as the one postulated at the scale of ÿbers. Such a feature shows the way any appropriate phenomenological continuous constitutive model should be chosen. Analogous results could be deduced from the analysis of case 1, but with no such simple interpretation, because the behavior of the equivalent continuum depends on two independent macroscopic ÿelds.
In the present work, an application of the method of homogenization of periodic discrete structures was presented. The method was shown to provide fundamental theoretical results on the structure of macroscopic constitutive equations suitable for a continuous modelling of a speciÿc net of ÿbers. It requires almost no restrictive physical assumptions and enables an easy, computer time e cient, analysis of the behavior for ÿber nets eventually including a great number of ÿbers, which is a necessary feature for the study of most ÿber-reinforced uids.
As is visible in the above theoretical exposition, many extensions to this work can be envisaged. Richer ÿber-ÿber interaction laws could ÿrst be introduced, as, for example, the case of Carreau type or viscoelastic interactions, or dry friction between ÿbers. The great adaptability of the method would also enable one to account for the exibility of ÿbers, considering them, for example, as elastic beams.

Appendix A
As shown in Section 4, in the general case, the state of stress of the ÿber net is deÿned by the three tensors 0 , Ä 0 and 0 , the vector 0 being given by the relation according to deÿnition (32) of L k . Let us now multiply this last equation by a constant virtual ÿeld 0 . We get Terms like p b 1 ∧ 0 can be seen as a virtual ÿeld b , so forces f k 0 being solutions of the self-equilibrium equation (86), one obtains which gives an alternate deÿnition of vector 0 as This result then enables us to ÿnd the relation between 0 and 0 , using the deÿnition and properties of the permutation tensor e. e is the tensor whose components in e i ⊗ e j ⊗ e k are It is then easy to check that 0 = g −1 k∈CR (e · (G · k )) · f k 0 = e : k∈CR g −1 (G · k ⊗ f k 0 ) which is equivalent to 0 = e : 0 :