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

Steven Le Corre, Denis Caillerie, Laurent Orgéas, Denis Favier

## ► To cite this version:

Steven Le Corre, Denis Caillerie, Laurent Orgéas, Denis Favier. Behavior of a net of fibers linked by viscous interactions: theory and mechanical properties. Journal of the Mechanics and Physics of Solids, 2004, 52 (2), pp.395-421. 10.1016/S0022-5096(03)00090-5. hal-00531571

HAL Id: hal-00531571

https://minesparis-psl.hal.science/hal-00531571v1

Submitted on 3 Nov 2010

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Behavior of a net of fibers linked by viscous interactions: theory and mechanical properties

S. Le Corre<sup>a,b,*</sup>, D. Caillerie<sup>a</sup>, L. Orgéas<sup>a</sup>, D. Favier<sup>a</sup>

This paper presents an investigation of the macroscopic mechanical behavior of highly con-
centrated fiber suspensions for which the mechanical behavior is governed by local fiber–fiber
interactions.

The problem is approached by considering the case of a net of rigid fibers 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 fibers, 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 fiber net is regarded as a periodic assembly of identi-
cal 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 effective Cosserat medium, the other ones being usual Cauchy media. Lastly, formulations that enable an effective computation of constitutive equations are detailed. They show that the equivalent continuum behaves like an anisotropic power-law fluid.

Keywords: Homogenization; Lattice; Constitutive model; Non-linear behavior; Fiber-reinforced composites

<footnote><sup>a</sup>Laboratoire Sols, Solides, Structures, CNRS-Institut National Polytechnique de Grenoble-Université
Joseph Fourier, Domaine Universitaire BP 53, 38041 Grenoble cedex 9, France</footnote>

<footnote><sup>b</sup>CEMEF, Ecole des Mines de Paris, 1 rue Claude Daunesse, BP 207, 06904 Sophia Antipolis cedex,
France</footnote>
# 1. Introduction

Understanding the behavior of short fiber-reinforced fluids is of great interest in many industrial applications such as the processing of polymeric materials and other composites. In short fiber systems processing, fibers noticeably influence the flow of the fluid matrix and conversely, the flow of the matrix determines the spatial distribution and orientation of fibers, 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 fibers, and macroscopic mechanical behavior, under some restrictive assumptions. They generally apply to the case of rigid straight fibers immersed in a Newtonian fluid when short range interactions between fibers 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 effects of local interactions between fibers, that may be due to dry friction effects 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 effective calculation of the behavior possible only in the case of perfectly unidirectional or perfectly isotropic orientations of fibers (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 fluids reinforced with fibers. This is due to the difficulty of calculating velocity fields around fibers as in the Newtonian case (Jeffery, 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 fibers, the contribution to the total stress of the matrix and of the fiber-matrix interaction becoming then negligible. The critical fiber concentration for such an assumption to be valid cannot be established in a general way. It largely depends on the fiber'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 fibers 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 flows such as biaxial extension or simple shear flow.

As a first step towards a more general approach to the modelling of the behavior of short fiber systems, this paper exposes an homogenization method suitable for the modelling of highly concentrated fiber 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 fibers (microscopic scale), the proposed method enables finding 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 fibers. 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 fiber net to be transformed into continuous ones: this leads to the definition of the macroscopic stress tensors of the equivalent continuum (see Section 4). Depending on the fiber-fiber interaction laws at the microscopic level, three types of macroscopic constitutive equations are obtained, the effective 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 fluid.

For the sake of simplicity, the problem is exposed in the case of a planar fiber net, but, as it will be clear in the equations, the extension to a three dimensional problem would be straightforward.

## 1.1. 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

$$ \underline{\mathbf{a}} \cdot \underline{\mathbf{b}} = a_i b_i, \quad (\underline{\mathbf{a}} \cdot \underline{\mathbf{b}})_i = a_{ij} b_j, $$

$$ (\underline{\mathbf{a}} \cdot \underline{\mathbf{b}})_{ij} = a_{ik} b_{kj}, \quad (\underline{\mathbf{a}} : \underline{\mathbf{b}}) = a_{ij} b_{ji} $$

using Einstein's summation convention over repeated indices.

Tensorial product is denoted by the symbol $\otimes$, e.g.:

$$ (\underline{\mathbf{a}} \otimes \underline{\mathbf{b}})_{ij} = a_i b_j, \quad (\underline{\mathbf{a}} \otimes \underline{\mathbf{b}})_{ijk} = a_{ij} b_k. $$

The gradient of a vector $\underline{\mathbf{a}}$ with respect to space variables $x_i$ will be denoted $\nabla \underline{\mathbf{a}}$ and defined as

$$ (\nabla \underline{\mathbf{a}})_{ij} = \frac{\partial a_i}{\partial x_j}. $$

The same convention will be used for higher order spatial derivatives, e.g.:

$$ (\nabla \underline{\mathbf{a}})_{ijk} = \frac{\partial^2 a_i}{\partial x_j \partial x_k}. $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Discrete and periodic description of the fiber net, example of reference cell.</figcaption></figure>

## 2. Discrete description and microscopic modelling

### 2.1. Basic assumptions

The structure under consideration is a planar net made of rigid cylindrical fibers of uniform length $l$, that will be called bars. The plane is characterized by the Cartesian reference frame $\mathcal{R}_\phi = (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 first assumption is that the structure may be considered as periodic. Thus the fiber net is regarded as an assembly of identical cells as shown in Fig. 1, which are characterized by the two vectors of $\mathcal{R}_\phi$, $\underline{d}_1$ and $\underline{d}_2$, called periodicity vectors. The fiber net dimensions are therefore $L_1$ and $L_2$ which are such that

$$L_1 = N_1 \|\underline{d}_1\| \quad \text{and} \quad L_2 = N_2 \|\underline{d}_2\|, \qquad (1)$$

where $N_1$ is the number of cells in the $\underline{d}_1$ direction and $N_2$ the number of cells in the $\underline{d}_2$ direction.

The second assumption is that the fiber net is made of a huge number of cells $\mathcal{N}^c$, so that the scale separation parameter $\varepsilon$, defined as

$$\varepsilon = \frac{1}{\sqrt{\mathcal{N}^c}} \ll 1, \qquad (2)$$

may be considered as a very small parameter. Condition (2) is equivalent to the following, and so can have another meaning:

$$\varepsilon = \sqrt{\frac{s}{S}} = \frac{d}{L} \ll 1, \qquad (3)$$

where $s$ and $S$ are respectively the surface of one cell and the surface of the whole fiber net, $d = \sqrt{s}$ represents the characteristic length of the microstructure and $L = \sqrt{s}$
is the characteristic length of the whole net. Let us assume that $d$ has the same order of magnitude as the length of fibers. Condition (3) is therefore a condition of scale separation; it implies that the size of the fiber net should be large in comparison with the size of fibers.

## 2.2. Numbering system

The net of fibers 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 fibers. The geometry and topology of the net is therefore entirely defined by the spatial and angular position of each fiber, and by the connectivity of each fiber with the rest of the assembly.

Centers of cells are first located by a vector of integers $\underline{a}$ whose components are $(a_1, a_2)$. Their positions are located by vectors $\underline{p}_0(\underline{a})$ (see Fig. 1):

$$ \underline{p}_0(\underline{a}) = a_1\underline{d}_1 + a_2\underline{d}_2. \qquad (4) $$

The assumption of periodicity then suggests a system of numbering of bars and links reflecting the regularity of the microstructure. Each bar of the net is numbered by $\tilde{b} = (b, \underline{a}^b)$, which means that the bar $\tilde{b}$ is the $b^{th}$ bar of the cell $\underline{a}^b$. The set all bars of the net is denoted $\mathcal{B}$.

In the same way, the action of the bar $\tilde{c}$ on the bar $\tilde{b}$ is denoted either by $\tilde{k} = (\tilde{c}/\tilde{b})$ or by $\tilde{k} = (k, \underline{a}^k)$, the set of connections of the net being denoted $\mathcal{C}$. Bars $\tilde{b}$ and $\tilde{c}$ can respectively be considered as the interior and the exterior of the interaction $\tilde{k}$, so the following notation will be used:

$$ \tilde{b} = I(\tilde{k}), \quad \tilde{c} = E(\tilde{k}) \quad \text{and so} \quad \tilde{k} = (E(\tilde{k})/I(\tilde{k})). \qquad (5) $$

We will consider that the cell to which $\tilde{k}$ belongs is the belonging cell of the bar on which the action is exerted, i.e. the bar $\tilde{b} = I(\tilde{k})$, so that $\underline{a}^k = \underline{a}^b$. The reciprocal interaction $(\tilde{b}/\tilde{c})$ will be denoted $^k\tilde{k} = (^k k, \underline{a}^k)$ and will therefore belong to the cell $\underline{a}^c$. In the notations $\tilde{b} = (b, \underline{a}^b)$ and $\tilde{k} = (k, \underline{a}^k)$, integers $b$ and $k$ are the numbers of bars and connections in a reference cell the sets of which are denoted $\mathcal{B}_R$ and $\mathcal{C}_R$.

In this work, we will assume that the size of cells $d$ is about the same order as the length of fibers $l$, and that $d > l$. This implies that interactions can take place only between two fibers of the same cell or between fibers of two neighboring cells. The exterior of the connection $\tilde{k} = (k, \underline{a}^k)$ will therefore be located in the cell $\underline{a}^k$ or in a neighboring one. Anyway, it will belong to the cell $\underline{a}^k + \delta^k$, where $\delta^k$ is a vector of integers whose components take their values in $\{-1, 0, +1\}$. Notice that the periodicity assumption causes $\delta^k$ to be independent of the position of the cell.

## 2.3. Geometrical description

As illustrated in Fig. 2, the position of a bar $\tilde{b}$ is defined by the position of its center $P^b$, located by the vector $\underline{p}(\tilde{b}) = \underline{OP}^b$, and by its unit vector $\underline{e}^b$. The periodicity of the
<figure><img src="image_3.png" /><figcaption>Fig. 2. Locating and kinematics of fibers in the planar case.</figcaption></figure>

net implies that $\underline{p}(\tilde{b})$ can be partitioned as

$$ \underline{p}(\tilde{b}) = \underline{p}(b, \underline{a}^{\tilde{b}}) = \underline{p}_{0}(\underline{a}^{\tilde{b}}) + \varepsilon \underline{p}_{1}^{b} = a_{1}\underline{d}_{1} + a_{2}\underline{d}_{2} + \varepsilon \underline{p}_{1}^{b}, \qquad (6) $$

where $\underline{p}_0(\underline{a})$ is the macroscopic position of the belonging cell of $\tilde{b}$ in the fiber net, and $\varepsilon \underline{p}_1^b$ is the local position of bar $b$ in that cell. It is to be noted that $\underline{p}_1^b$ 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 fiber 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 $\underline{d}_i$ are small compared to the size $L$ of the net, so they can be written as $\underline{d}_i = \varepsilon \underline{\Lambda}_i$ where the vectors $\underline{\Lambda}_i$ are macroscopic vectors, independent of $\varepsilon$. Let us now introduce the pseudo continuous variable $\underline{\lambda}^b$ defined by $\underline{\lambda}^b = (\lambda_1, \lambda_2) = \varepsilon \underline{a}^{\tilde{b}}$. According to this notation, Eq. (6) becomes

$$ \underline{p}(\tilde{b}) = \lambda_1 \underline{\Lambda}_1 + \lambda_2 \underline{\Lambda}_2 + \varepsilon \underline{p}_1^b = \underline{p}_0(\underline{\lambda}) + \varepsilon \underline{p}_1^b. \qquad (7) $$

As illustrated in Fig. 3, macroscopic positions of bars are now parameterized by the vector of reals $\underline{\lambda}$ which is a variable of $\Lambda \subset [0, 1] \times [0, 1]$. $\Lambda$ is the reference parametric space, where the net is made of square cells of size $\varepsilon$. As $\varepsilon$ is assumed to tend to zero in the use of asymptotic expansion methods, $\underline{\lambda}$ tends to become, and will be used as a continuous variable, even if, strictly speaking, it only takes values such as $\underline{\lambda} = \varepsilon \underline{a}$ for a finite value of $\varepsilon$.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Parameterization of the geometry of the fiber net. Discrete description in the physical space; continuous description in the parametric space.</figcaption></figure>

Let us introduce $\underline{G}$, the gradient $\underline{\nabla} \underline{p}_0$ of the geometrical transformation $\underline{\lambda} \mapsto \underline{p}_0(\underline{\lambda})$. According to Eq. (7), the components of $\underline{G}$ relative to $\underline{e}_i \otimes \underline{e}_j$ can be expressed as

$$G_{ij} = (\underline{\Lambda}_i \cdot \underline{e}_j) \qquad (8)$$

and the Jacobian $g$ of the transformation is

$$g = \det \underline{G} = \| \underline{\Lambda}_1 \wedge \underline{\Lambda}_2 \|. \qquad (9)$$

## 2.4. Local kinematics and interactions

The problem being planar in $\mathcal{R}_\phi$ and fibers being considered as rigid bars, their motion is a rigid body planar motion. Kinematics of a bar $\tilde{b}$ can therefore be defined by:

* the velocity of its center $P^{\tilde{b}}$, denoted $v(\tilde{b})$

* its angular velocity $\omega(\tilde{b}) = \omega(\tilde{b})\underline{e}_3$, where $\underline{e}_3$ is the unit vector normal to the plane of fibers.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Locating of a contacting point along a bar $\tilde{b}$.</figcaption></figure>

The velocity of a point $Q$ of bar $\tilde{b}$ is then

$$ \underline{\mathbf{v}}^{\tilde{b}}(Q) = \underline{\mathbf{v}}(\tilde{b}) + \underline{\omega}(\tilde{b}) \wedge \underline{\mathbf{P}}^{\tilde{b}}\underline{\mathbf{Q}} = \underline{\mathbf{v}}(\tilde{b}) + \underline{\omega}(\tilde{b}) \wedge \rho \underline{\mathbf{e}}^{b}, \quad (10) $$

where $\rho$ is the curvilinear abscissa of $Q$ on the bar $\tilde{b}$ with respect to the center $P^{\tilde{b}}$ as shown in Fig. 4. $\rho$ is a small variable; its order of magnitude is about the size of the cell $d$ so it can be written as $\rho = d\xi = \varepsilon L\xi$, where $\xi$ is independent of the size of the problem $L$. In the following, it will be convenient to introduce the new variable $\underline{\phi} = L\underline{\omega}$, which enables us to rewrite Eq. (10) as

$$ \underline{\mathbf{v}}^{\tilde{b}}(Q) = \underline{\mathbf{v}}(\tilde{b}) + \varepsilon\xi\underline{\phi}(\tilde{b}) \wedge \underline{\mathbf{e}}^{b}. \quad (11) $$

The action $\tilde{k}$ of bar E($\tilde{k}$) on bar I($\tilde{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 interaction $\tilde{k}$ can therefore be partitioned into:

* a force $\underline{f}^{\tilde{k}}$ proportional to the difference of velocity of the two bars $\underline{\Delta v}^{\tilde{k}}$ at their contacting point, denoted $Q^{\tilde{k}}$:

$$ \underline{f}^{\tilde{k}} = \mu^k \|\underline{\Delta v}^{\tilde{k}}\|^m \underline{\Delta v}^{\tilde{k}}, \quad (12) $$

where the power law index $m$ is a real scalar ranging from 0 to 1. By definition, $\underline{\Delta v}^{\tilde{k}}$ equals the difference $\underline{v}^{E(\tilde{k})}(Q^{\tilde{k}}) - \underline{v}^{I(\tilde{k})}(Q^{\tilde{k}})$, so making use of Eq. (11), its expression reads

$$ \underline{\Delta v}^{\tilde{k}} = \underline{\mathbf{v}}(E(\tilde{k})) - \underline{\mathbf{v}}(I(\tilde{k})) + \varepsilon(\xi^{\tilde{k}}\underline{\phi}(E(\tilde{k})) \wedge \underline{\mathbf{e}}^{E(\tilde{k})} - \xi^{\tilde{k}}\underline{\phi}(I(\tilde{k})) \wedge \underline{\mathbf{e}}^{I(\tilde{k})}). \quad (13) $$

In the last equation, $\xi^{\tilde{k}}$ (resp. $\xi^{\tilde{k}}$) is the normalized local abscissa of $Q^{\tilde{k}}$ on the bar $\tilde{b} = I(\tilde{k})$ (resp. $\tilde{c} = E(\tilde{k})$).

* a moment $\underline{M}^{\tilde{k}}(Q^{\tilde{k}})$ relative to $Q^{\tilde{k}}$, proportional to the difference of angular velocities of both bars:

$$ \underline{M}^{\tilde{k}}(Q^{\tilde{k}}) = B^k \|\underline{\Delta \omega}^{\tilde{k}}\|^m \underline{\Delta \omega}^{\tilde{k}}, \quad (14) $$

where $\underline{\Delta \omega}^{\tilde{k}}$ denotes the difference of angular velocities $\underline{\omega}(E(\tilde{k})) - \underline{\omega}(I(\tilde{k}))$.
In order to normalize moments with respect to forces, vectors $\underline{m}^k(Q^k)$, such as $\underline{M}^k(Q^k)$ = $L\underline{m}^k(Q^k)$, 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

$$ \underline{m}^k(Q^k) = \beta^k \|\underline{\Delta\phi}^k\|^{m-1} \underline{\Delta\phi}^k, \quad (15) $$

where $\beta^k = B^k/L^{m+1}$. The moment of action $\tilde{k}$ expressed at point $P^k$, the center of bar $I(\tilde{k})$, will be denoted $\underline{m}^k$ and will have the following expression:

$$ \underline{m}^k = \beta^k \|\underline{\Delta\phi}^k\|^{m-1} \underline{\Delta\phi}^k + \epsilon \zeta^k \underline{e}^l(\tilde{k}) \wedge \underline{f}^k. \quad (16) $$

In Eqs. (12) and (15), the coefficients $\mu^k$ and $\beta^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

$$ \mu^k = \epsilon^{1-m} \mu_0^k \quad (17) $$

and

$$ \beta^k = \epsilon^{1+q-m} \beta_0^k, \quad (18) $$

where $\mu_0^k$ and $\beta_0^k$ 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 $\beta^k$ with respect to $\mu^k$. In the following, only cases where $q \ge 0$ will be examined, for situations where $\beta^k$ would be much greater than $\mu^k$ are not physically likely to happen in the considered suspensions. For simplicity, three cases will be discussed. They may represent three different situations one might expect when studying a specific 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 $\beta_0^k = 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.

### 3. 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 fiber suspensions (Batchelor, 1971; Gibson and Toll, 1999).
Such a process would require a further assumption on the velocity fields 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.

## 3.1. Asymptotic expansions

### 3.1.1. Velocity fields

Taking advantage of the smallness of the parameter $\epsilon$ and of the assumption of periodicity, velocity and angular velocity fields relative to a bar $\tilde{b}$ are expanded in discrete asymptotic series expansion in powers of $\epsilon$ (Moreau and Caillerie, 1995; Tollenaere and Caillerie, 1998). Thus, they are written as

$$ \mathbf{v}(\tilde{b}) = \mathbf{v}^b(\tilde{\lambda}) = \mathbf{v}_0^b(\tilde{\lambda}) + \epsilon \mathbf{v}_1^b(\tilde{\lambda}) + \epsilon^2 \mathbf{v}_2^b(\tilde{\lambda}) + \dots, \quad (19) $$

$$ \underline{\phi}(\tilde{b}) = \underline{\phi}^b(\underline{\lambda}) = \underline{\phi}_0^b(\underline{\lambda}) + \epsilon \underline{\phi}_1^b(\underline{\lambda}) + \epsilon^2 \underline{\phi}_2^b(\underline{\lambda}) + \dots, \quad (20) $$

where functions $\mathbf{v}_n^b$ and $\underline{\phi}_n^b$ are continuous functions of $\underline{\lambda}$ that generally depend on the bar b. $\underline{\lambda}$ reflects the macroscopic variation of those functions whereas index b is a local variable, depending on the corresponding bar of $\tilde{b}$ in the reference cell.

For a pair of connected bars $\tilde{k} = (\tilde{c}, \tilde{b})$, let us recall that $\tilde{b}$ is supposed to belong to the cell located by $\underline{\lambda}$, whereas $\tilde{c}$ may belong to a neighboring cell located by $\underline{\lambda} + \epsilon \underline{\delta}$. Thus $\mathbf{v}(\tilde{c}) = \mathbf{v}^{E(k)}(\underline{\lambda} + \epsilon \underline{\delta}^k)$ is expanded using a Taylor expansion around $\underline{\lambda}$, using the fact that $\epsilon \underline{\delta}^k$ is very small with respect to $\underline{\lambda}$:

$$ \begin{aligned} \mathbf{v}(\tilde{c}) &= \mathbf{v}^{E(k)}(\underline{\lambda} + \epsilon \underline{\delta}^k) \\ &= \mathbf{v}^{E(k)}(\underline{\lambda}) + \epsilon \underline{\nabla} \mathbf{v}^{E(k)}(\underline{\lambda}) \cdot \underline{\delta}^k + \epsilon^2 \frac{1}{2!} \underline{\nabla} \mathbf{v}^{E(k)}(\underline{\lambda}) : [\underline{\delta}^k \otimes \underline{\delta}^k] + \dots \end{aligned} \quad (21) $$

Carrying asymptotic expansion (19) into Eq. (21) and identifying each order of powers of $\epsilon$, the expansion of $\mathbf{v}(\tilde{c})$ reads

$$ \begin{aligned} \underline{\mathbf{v}}(\tilde{c}) = \underline{\mathbf{v}}^{E(k)}(\underline{\lambda} + \epsilon \underline{\delta}^k) = \underline{\mathbf{v}}_0^{E(k)}(\underline{\lambda}) + \epsilon(\underline{\mathbf{v}}_1^{E(k)}(\underline{\lambda}) + \underline{\nabla} \underline{\mathbf{v}}_0^{E(k)}(\underline{\lambda}) \cdot \underline{\delta}^k) \\ + \epsilon^2 \left( \underline{\mathbf{v}}_2^{E(k)}(\underline{\lambda}) + \underline{\nabla} \underline{\mathbf{v}}_1^{E(k)}(\underline{\lambda}) \cdot \underline{\delta}^k + \frac{1}{2!} \underline{\nabla} \underline{\mathbf{v}}_0^{E(k)}(\underline{\lambda}) : [\underline{\delta}^k \otimes \underline{\delta}^k] \right) + \dots \end{aligned} \quad (22) $$

The expansion of the angular velocity $\underline{\phi}(\tilde{c})$ is obtained exactly in the same manner.

### 3.1.2. Velocity difference fields

In the same way, velocity differences $\underline{\Delta} \mathbf{v}^{\tilde{k}}(13)$ and $\underline{\Delta} \underline{\phi}^{\tilde{k}} = \underline{\phi}^{E(\tilde{k})} - \underline{\phi}^{I(\tilde{k})}$ can be expanded:

$$ \begin{aligned} \underline{\Delta} \mathbf{v}^{\tilde{k}} &= \underline{\Delta} \mathbf{v}_0^{\tilde{k}} + \epsilon \underline{\Delta} \mathbf{v}_1^{\tilde{k}} + \epsilon^2 \underline{\Delta} \mathbf{v}_2^{\tilde{k}} + \dots, \\ \underline{\Delta} \underline{\phi}^{\tilde{k}} &= \underline{\Delta} \underline{\phi}_0^{\tilde{k}} + \epsilon \underline{\Delta} \underline{\phi}_1^{\tilde{k}} + \epsilon^2 \underline{\Delta} \underline{\phi}_0^{\tilde{k}} + \dots \end{aligned} \quad (23) $$
By substituting expansion (19) into Eq. (13), and making use of Taylor expansions (22), velocity differences $\underline{\Delta v}_i^k$ can be identified at each order of powers of $\epsilon$. For the two first orders ($\epsilon^0$ and $\epsilon^1$), we obtain

$$ \underline{\Delta v}_0^k = \underline{v}_0^{E(k)}(\underline{\lambda}) - \underline{v}_0^{I(k)}(\underline{\lambda}), \quad (24) $$

$$ \begin{aligned} \underline{\Delta v}_1^k = \underline{v}_1^{E(k)}(\underline{\lambda}) - \underline{v}_1^{I(k)}(\underline{\lambda}) &+ \xi^k \underline{\phi}_0^{E(k)}(\underline{\lambda}) \wedge \underline{\epsilon}^{E(k)} - \xi^k \underline{\phi}_0^{I(k)}(\underline{\lambda}) \wedge \underline{\epsilon}^{I(k)} \\ &+ \underline{\nabla} \underline{v}_0^{E(k)}(\underline{\lambda}) \cdot \underline{\delta}^k. \end{aligned} \quad (25) $$

For the two first orders of angular velocities difference, one obtains

$$ \underline{\Delta \phi}_0^k = \underline{\phi}_0^{E(k)}(\underline{\lambda}) - \underline{\phi}_0^{I(k)}(\underline{\lambda}), \quad (26) $$

$$ \underline{\Delta \phi}_1^k = \underline{\phi}_1^{E(k)}(\underline{\lambda}) - \underline{\phi}_1^{I(k)}(\underline{\lambda}) + \underline{\nabla} \underline{\phi}_0^{E(k)}(\underline{\lambda}) \cdot \underline{\delta}^k. \quad (27) $$

## 3.2. Equilibrium of the net

In the scope of this paper, inertia effects 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 bar $\tilde{b}$ can be described by the two following sets of equations:

$$ \forall \tilde{b}, \quad \sum_{\tilde{k} \in \mathcal{C}(\tilde{b})} \underline{f}^{\tilde{k}} = \underline{0} \quad \text{and} \quad \sum_{\tilde{k} \in \mathcal{C}(\tilde{b})} \underline{m}^{\tilde{k}} = \underline{0}, \quad (28) $$

where $\mathcal{C}(\tilde{b})$ is the set of connections of bar $\tilde{b}$.

The expansion of equilibrium equations of the whole fiber net can be simplified 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 $\epsilon$. Thus using two sets of virtual fields denoted $(\underline{v}$ and $\underline{\varphi})$, problem (28) is equivalent to the following virtual formulation:

$$ \forall (\underline{v}(\tilde{b}), \underline{\varphi}(\tilde{b})), \quad \sum_{\tilde{b}} \sum_{\tilde{k} \in \mathcal{C}(\tilde{b})} \underline{f}^{\tilde{k}} \cdot \underline{v}(\tilde{b}) + \underline{m}^{\tilde{k}} \cdot \underline{\varphi}(\tilde{b}) = 0. \quad (29) $$

Within such a summation, it can be noticed that if the connection $\tilde{k} = (\tilde{c}/\tilde{b})$ belongs to $\mathcal{C}(\tilde{b})$, connection $\tilde{k} = (\tilde{b}/\tilde{c})$ belongs to $\mathcal{C}(\tilde{c})$. Summation on the bars can therefore be transformed into a summation on the set of connected pairs $(\tilde{c}/\tilde{b})$ such that $\tilde{c} > \tilde{b}$. Doing so, formulation (29) can rewritten as

$$ \forall (\underline{v}(\tilde{b}), \underline{\varphi}(\tilde{b})), \quad \sum_{\tilde{k}=(\tilde{c}, \tilde{b}) \atop \tilde{c}>\tilde{b}} \left[ \underline{f}^{\tilde{k}} \cdot (\underline{v}(\underline{E}(\tilde{k})) - \underline{v}(\underline{I}(\tilde{k}))) + \underline{m}^{\tilde{k}} \cdot (\underline{\varphi}(\underline{E}(\tilde{k})) - \underline{\varphi}(\underline{I}(\tilde{k}))) - \epsilon \underline{L}^{\tilde{k}} \wedge \underline{f}^{\tilde{k}} \cdot \underline{\varphi}(\underline{E}(\tilde{k})) \right] = 0. \quad (30) $$
This last equation was obtained using the action and reaction theorem which in this case can be shown to imply

$$ \underline{f}^{\tilde{k}} + \underline{f}'^{\tilde{k}} = \underline{0} \quad \text{and} \quad \underline{m}^{\tilde{k}} + \underline{m}'^{\tilde{k}} - \varepsilon \underline{L}^{\tilde{k}} \wedge \underline{f}^{\tilde{k}} = \underline{0}, \qquad (31) $$

where vector $\varepsilon \underline{L}^k$ denotes the vector $\underline{P}^{\tilde{b}} \underline{P}^{\tilde{c}}$ and is therefore defined as

$$ \varepsilon \underline{L}^k = \underline{p}(\tilde{c}) - \underline{p}(\tilde{b}) = \varepsilon(\underline{G} \cdot \underline{\delta}^k + \underline{p}_1^c - \underline{p}_1^b). \qquad (32) $$

### 3.3. From discrete to continuous

Suppose that $\underline{u}^{\tilde{b}}$ is any vectorial field of the spatial variable $\underline{\lambda}$ with the following two properties:

$$ \begin{gathered} \sum_{\tilde{b} \in \mathcal{B}} \underline{u}^{\tilde{b}} = \sum_{\underline{a}} \sum_{b \in \mathcal{B}_R} \underline{u}^b, \\ \underline{u}^{\tilde{b}}(\underline{\lambda}) = \varepsilon^2 \underline{u}_0^b(\underline{\lambda}) + \mathcal{O}(\varepsilon^3). \end{gathered} $$

Then it can be shown that when $\varepsilon$ tends to zero, we have

$$ \lim_{\varepsilon \to 0} \sum_{\tilde{b} \in \mathcal{B}} \underline{u}^{\tilde{b}} = \int_{\Lambda} \sum_{b \in \mathcal{B}_R} \underline{u}_0^b(\underline{\lambda}), \qquad (33) $$

using the definition 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.

### 3.4. Self-equilibrium at the lower orders

In a first stage let us assume that both order zero velocities $\underline{v}_0^b$ and angular velocities $\underline{\phi}_0^b$ actually depend on the considered bar, so that order zero velocity differences $\underline{\Delta v}_0^k$ and $\underline{\Delta \phi}_0^k$ are non null functions.

#### 3.4.1. Self-equilibrium of forces

According to asymptotic expansions (23), expressions such as $||\underline{\Delta v}^{\tilde{k}}||^{m-1}$ have an asymptotic expansion of the following form:

$$ \begin{aligned} ||\underline{\Delta v}^{\tilde{k}}||^{m-1} &= ||\underline{\Delta v}_0^k + \varepsilon \underline{\Delta v}_1^k + \varepsilon^2 \underline{\Delta v}_2^k + \cdots||^{m-1} \\ &= ||\underline{\Delta v}_0^k||^{m-1} \left( 1 + \varepsilon(m-1) \frac{\underline{\Delta v}_0^k \cdot \underline{\Delta v}_1^k}{||\underline{\Delta v}_0^k||^2} + \cdots \right). \end{aligned} \qquad (34) $$
Subsequently, at this stage, the asymptotic expansions of interaction forces, defined by Eq. (12), can be written as

$$ \underline{f}^{\tilde{k}} = \varepsilon^{1-m} \mu_{0}^{k} \|\underline{\Delta} \underline{v}_{0}^{k}\|^{m-1} \left[ \underline{\Delta} \underline{v}_{0}^{k} + \varepsilon \left( \underline{\Delta} \underline{v}_{1}^{k} + (m-1) \frac{\underline{\Delta} \underline{v}_{0}^{k} \cdot \underline{\Delta} \underline{v}_{1}^{k}}{\|\underline{\Delta} \underline{v}_{0}^{k}\|^{2}} \underline{\Delta} \underline{v}_{0}^{k} \right) + \cdots \right]. \quad (35) $$

Their asymptotic expansion may therefore be written in the following way:

$$ \underline{f}^{\tilde{k}} = \varepsilon^{2-m} \left( \frac{1}{\varepsilon} \underline{f}_{-1}^{k} + \underline{f}_{0}^{k} + \varepsilon \underline{f}_{1}^{k} + \cdots \right) \quad (36) $$

with, at the lower order of $\varepsilon$ powers,

$$ \underline{f}_{-1}^{k} = \mu_{0}^{k} \|\underline{\Delta} \underline{v}_{0}^{k}\|^{m-1} \underline{\Delta} \underline{v}_{0}^{k}. \quad (37) $$

Taking into the virtual power formulation (30) virtual velocity fields $\underline{v}(\tilde{b})$ and $\underline{\varphi}(\tilde{b})$ such that

$$ \forall \tilde{b}, \quad \underline{v}(\tilde{b}) = \varepsilon^{1+m} \alpha_{1}(\underline{\lambda}) \underline{v}^{b} \quad \text{and} \quad \underline{\varphi}(\tilde{b}) = 0 \quad (38) $$

then carrying into Eq. (30) asymptotic expansions (36) and making use of property (33), we obtain, when $\varepsilon$ tends to 0:

$$ \forall (\alpha_1(\underline{\lambda}), \underline{v}^b), \quad \int_A \sum_{k \in \mathcal{C}_R} \underline{f}_{-1}^k \cdot (\underline{v}^{E(k)} - \underline{v}^{I(k)}) \alpha_1(\underline{\lambda}) \, d\lambda = 0. \quad (39) $$

This relation being verified for any continuous function $\alpha_1$, forces of the lower order are solution of the lower order self-equilibrium equation of the reference cell:

$$ \forall \underline{v}^b, \quad \sum_{k \in \mathcal{C}_R} \underline{f}_{-1}^k \cdot (\underline{v}^{E(k)} - \underline{v}^{I(k)}) = 0. \quad (40) $$

Then carrying into this last equation the constitutive relationship (37), we obtain a relation that must be satisfied by the first order velocities of bars in the reference cell. Velocities $\underline{v}_0^b$ are actually solutions of the problem

$$ \forall \underline{v}^b, \quad \sum_{k \in \mathcal{C}_R} \mu_0^k \|\underline{v}_0^{E(k)} - \underline{v}_0^{I(k)}\|^{m-1} (\underline{v}_0^{E(k)} - \underline{v}_0^{I(k)}) \cdot (\underline{v}^{E(k)} - \underline{v}^{I(k)}) = 0. \quad (41) $$

In the special case where $\underline{v}^b = \underline{v}_0^b$, one checks that order zero velocities have to satisfy the condition

$$ \sum_{k \in \mathcal{C}_R} \mu_0^k \|\underline{v}_0^{E(k)} - \underline{v}_0^{I(k)}\|^{m+1} = 0. \quad (42) $$

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:

$$ \forall b \in \mathcal{B}_R, \quad \underline{v}_0^b(\underline{\lambda}) = \underline{v}_0. \quad (43) $$

According to such a property, we also deduce that

$$ \forall k \in \mathcal{C}_R, \quad \underline{f}_{-1}^k = \underline{0}. \quad (44) $$
Such a result is however inconsistent with the asymptotic expansion (36) proposed for interaction forces which was based on the assumption that $\underline{\Delta v}_0^k \neq 0$. Actually, according to property (43), expressions such as $\|\underline{\Delta v}^k\|^{m-1}$ have an asymptotic expansion of the form:

$$ \|\underline{\Delta v}^k\|^{m-1} = \epsilon^m [\|\underline{\Delta v}_1^k\|^{m-1} + \epsilon(m-1)\|\underline{\Delta v}_1^k\|^{m-3}(\underline{\Delta v}_1^k \cdot \underline{\Delta v}_2^k) + \cdots]. \quad (45) $$

Subsequently, asymptotic expansions of interaction forces can be written as

$$ \underline{f}^k = \epsilon(\underline{f}_0^k + \epsilon\underline{f}_0^k + \epsilon^2\underline{f}_0^k + \cdots) \quad (46) $$

with, at order $\epsilon^0$,

$$ \underline{f}_0^k = \mu_0^k \|\underline{\Delta v}_1^k\|^{m-1} \underline{\Delta v}_1^k, \quad (47) $$

$$ \underline{\Delta v}_1^k = \underline{v}_1^{E(k)} - \underline{v}_1^{I(k)} + \underline{\nabla} \underline{v}_0 \cdot \underline{\delta}^k + \zeta^k \underline{\phi}_0^{E(k)} \wedge \underline{e}^{E(k)} - \zeta^k \underline{\phi}_0^{I(k)} \wedge \underline{e}^{I(k)}. \quad (48) $$

### 3.4.2. 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:

$$ \begin{aligned} \underline{m}^k ={}& \epsilon^{1+q-m} \beta_0^k [\|\underline{\Delta \phi}_0^k\|^{m-1} \underline{\Delta \phi}_0^k \\ & + \epsilon(\|\underline{\Delta \phi}_0^k\|^{m-1} \underline{\Delta \phi}_1^k + (m-1)\|\underline{\Delta \phi}_0^k\|^{m-3}(\underline{\Delta \phi}_0^k \cdot \underline{\Delta \phi}_1^k)\underline{\Delta \phi}_0^k) + \cdots] \\ & + \epsilon^2 \zeta^k \underline{e}^{I(k)} \wedge (\underline{f}_0^k + \epsilon \underline{f}_1^k + \epsilon^2 \underline{f}_0^k + \cdots). \end{aligned} \quad (49) $$

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 fields $\underline{v}(\tilde{b})$ and $\underline{\phi}(\tilde{b})$ such that

$$ \forall \tilde{b}, \quad \underline{v}(\tilde{b}) = 0 \quad \text{and} \quad \underline{\phi}(\tilde{b}) = \epsilon^{1+m-q} \alpha_2(\underline{\lambda}) \underline{\phi}^b \quad (50) $$

and then carrying into Eq. (30) asymptotic expansions (49), making use again of property (33), we obtain, when $\epsilon$ tends to 0

$$ \forall (\alpha_2(\underline{\lambda}), \underline{\phi}^b), \quad \int_{\Lambda} \sum_{k \in \mathcal{C}_R} \beta_0^k \|\underline{\Delta \phi}_0^k\|^{m-1} \underline{\Delta \phi}_0^k \cdot (\underline{\phi}^{E(k)} - \underline{\phi}^{I(k)}) \alpha_2(\underline{\lambda}) \, d\underline{\lambda} = 0. \quad (51) $$

According to property (33), all the higher order terms vanish when $\epsilon$ tends to zero. The last relation being verified for any function $\alpha_2$, order zero angular velocities are then shown to solve the following self-equilibrium problem on the reference cell:

$$ \forall \underline{\phi}^b, \quad \sum_{k \in \mathcal{C}_R} \beta_0^k \|\underline{\Delta \phi}_0^{E(k)} - \underline{\phi}_0^{I(k)}\|^{m-1} (\underline{\phi}_0^{E(k)} - \underline{\phi}_0^{I(k)}) \cdot (\underline{\phi}^{E(k)} - \underline{\phi}^{I(k)}) = 0. \quad (52) $$
This problem is strictly equivalent to the problem satisfied by order zero velocities so, for cases 1 and 2, we also obtain

$$ \forall b \in \mathcal{B}_R, \quad \underline{\phi}_0^b(\underline{\lambda}) = \underline{\phi}_0. \qquad (53) $$

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

$$ \underline{m}^k = \varepsilon(\underline{m}_0^k + \varepsilon\underline{m}_1^k + \varepsilon^2\underline{m}_2^k + \dots) \qquad (54) $$

with the following expressions of first order terms:

• Case 1:

$$ \underline{m}_0^k = \beta_0^k \|\underline{\Delta\phi}_1^k\|^m \underline{\Delta\phi}_1^k, \qquad (55) $$

$$ \underline{\Delta\phi}_1^k = \underline{\phi}_1^{\text{E}(k)} - \underline{\phi}_1^{\text{I}(k)} + \underline{\nabla}\underline{\phi}_0 \cdot \underline{\delta}^k. \qquad (56) $$

• Case 2:

$$ \underline{m}_0^k = 0, \quad \underline{m}_1^k = \beta_0^k \|\underline{\Delta\phi}_1^k\|^m \underline{\Delta\phi}_1^k, \qquad (57) $$

$$ \underline{\Delta\phi}_1^k = \underline{\phi}_1^{\text{E}(k)} - \underline{\phi}_1^{\text{I}(k)} + \underline{\nabla}\underline{\phi}_0 \cdot \underline{\delta}^k. \qquad (58) $$

• Case 3:

$$ \underline{m}_0^k = 0, \quad \underline{m}_1^k = \beta_0^k \|\underline{\Delta\phi}_0^k\|^m \underline{\Delta\phi}_0^k + \xi^k \underline{e}^{\text{I}(\tilde{k})} \wedge \underline{f}_0^k, \qquad (59) $$

$$ \underline{\Delta\phi}_0^k = \underline{\phi}_0^{\text{E}(k)} - \underline{\phi}_0^{\text{I}(k)}. \qquad (60) $$

# 4. Equilibrium of the equivalent continuum

Before detailing the obtaining of the equivalent continuum of the fiber 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.

## 4.1. Continuous media with internal rotation

A medium for which microscopic kinematics may imply rotational motion such as granular materials, foams or fiber 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 field $v(\underline{x})$
and an angular velocity field $\phi(\mathbf{x})$ and, as described by Eringen (1968) and Germain (1973), the virtual power formulation of its equilibrium reads

$$ 
\begin{aligned}
\forall (\underline{\mathbf{v}}^*, \underline{\phi}^*), \quad & - \int_{\Omega} (\underline{\boldsymbol{\sigma}} : \underline{\nabla} \underline{\mathbf{v}}^* + \underline{\zeta} \cdot \underline{\phi}^* + \underline{\mathbf{\kappa}} : \underline{\nabla} \underline{\phi}^*) d\Omega \\
& + \int_{\Omega} (\underline{\boldsymbol{f}} \cdot \underline{\mathbf{v}}^* + \underline{\boldsymbol{c}} \cdot \underline{\phi}^*) d\Omega = 0,
\end{aligned} 
\quad (61) $$

where $\underline{\mathbf{v}}^*$ and $\underline{\phi}^*$ are virtual velocity and angular velocity fields. In Eq. (61), $\underline{\boldsymbol{\sigma}}$ represents the stress tensor, $\underline{\mathbf{\kappa}}$ the couple-stress tensor (due to local rotations) and $\underline{\zeta}$ the micro-stresses vector. Vectors $\underline{\boldsymbol{f}}$ and $\underline{\boldsymbol{c}}$ respectively denote the external volume forces and moments. The principle of objectivity causes $\underline{\zeta}$ to be directly linked to $\underline{\boldsymbol{\sigma}}$ by the relation:

$$ \underline{\zeta} = \underline{\boldsymbol{e}} : \underline{\boldsymbol{\sigma}} = \underline{\boldsymbol{e}} : \underline{\boldsymbol{\sigma}}^A \quad (62) $$

tensor $\underline{\boldsymbol{e}}$ being the permutation tensor whose definition is given in Appendix A and $\underline{\boldsymbol{\sigma}}^A$ is the antisymmetric part of $\underline{\boldsymbol{\sigma}}$. Then by splitting $\underline{\nabla} \underline{\mathbf{v}}^*$ into its symmetric and antisymmetric parts $\underline{\mathbf{D}}^*$ and $\underline{\mathbf{R}}^*$ and accounting for relation (62), formulation (61) may also be rewritten as

$$ 
\begin{aligned}
\forall (\underline{\mathbf{v}}^*, \underline{\phi}^*), \quad & - \int_{\Omega} (\underline{\boldsymbol{\sigma}}^S : \underline{\mathbf{D}}^* + \underline{\boldsymbol{\sigma}}^A : (\underline{\mathbf{R}}^* - \underline{\boldsymbol{e}} \cdot \underline{\phi}^*) + \underline{\boldsymbol{\kappa}} : \underline{\nabla} \underline{\phi}^*) d\Omega \\
& + \int_{\Omega} \underline{\boldsymbol{f}} \cdot \underline{\mathbf{v}}^* + \underline{\boldsymbol{c}} \cdot \underline{\phi}^* d\Omega = 0.
\end{aligned} 
\quad (63) $$

Such a formulation leads to the following local equilibrium equations:

$$ \mathbf{div}(\underline{\boldsymbol{\sigma}}^S + \underline{\boldsymbol{\sigma}}^A) + \underline{\boldsymbol{f}} = \underline{\mathbf{0}}, \quad \mathbf{div}(\underline{\boldsymbol{\kappa}} + \underline{\boldsymbol{e}} : \underline{\boldsymbol{\sigma}}^A + \underline{\boldsymbol{c}} = \underline{\mathbf{0}} \quad (64) $$

and implies, for a viscous material, constitutive equations of the generic type:

$$ 
\begin{aligned}
\underline{\boldsymbol{\sigma}}^S &= \mathcal{S}(\underline{\mathbf{D}}, \underline{\mathbf{R}} - \underline{\boldsymbol{e}} \cdot \underline{\boldsymbol{\phi}}, \underline{\nabla} \underline{\boldsymbol{\phi}}), \\
\underline{\boldsymbol{\sigma}}^A &= \mathcal{A}(\underline{\mathbf{D}}, \underline{\mathbf{R}} - \underline{\boldsymbol{e}} \cdot \underline{\boldsymbol{\phi}}, \underline{\nabla} \underline{\boldsymbol{\phi}}), \\
\underline{\boldsymbol{\kappa}} &= \mathcal{K}(\underline{\mathbf{D}}, \underline{\mathbf{R}} - \underline{\boldsymbol{e}} \cdot \underline{\boldsymbol{\phi}}, \underline{\nabla} \underline{\boldsymbol{\phi}}).
\end{aligned} 
\quad (65) $$

In the special case where $\underline{\boldsymbol{\kappa}}$ and external moments $\underline{\boldsymbol{c}}$ actually happen to be null or negligible, $\underline{\nabla} \underline{\boldsymbol{\phi}}$ no longer influences the behavior. It is then clear from Eqs. (64) that vector $\underline{\boldsymbol{e}} \cdot \underline{\boldsymbol{\sigma}}^A$ is also null or negligible. In such a case, constitutive equations (65) simplify to

$$ \underline{\boldsymbol{\sigma}}^S = \mathcal{S}(\underline{\mathbf{D}}, \underline{\mathbf{R}} - \underline{\boldsymbol{e}} \cdot \underline{\boldsymbol{\phi}}), \quad (66) $$

$$ \underline{\boldsymbol{\sigma}}^A = \mathcal{A}(\underline{\mathbf{D}}, \underline{\mathbf{R}} - \underline{\boldsymbol{e}} \cdot \underline{\boldsymbol{\phi}}) = \underline{\mathbf{0}}. \quad (67) $$

The symmetry condition (67) then provides a new relation that theoretically enables $\underline{\mathbf{R}} - \underline{\boldsymbol{e}} \cdot \underline{\boldsymbol{\phi}}$ to be expressed as a function $\mathcal{F}$ of the strain rate tensor $\underline{\mathbf{D}}$. The stress tensor
may therefore be expressed as a function of $\underline{D}$ only:

$$ \underline{\sigma} = \underline{\sigma}^S = \mathcal{S}(\underline{D}, \mathcal{F}(\underline{D})) = \underline{\sigma}(\underline{D}). \qquad (68) $$

## 4.2. Equilibrium of the fiber net

For the fiber 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 fields $\underline{v}$ and $\underline{\varphi}$ such that

$$ \forall \tilde{b}, \quad \underline{v}(\tilde{b}) = \underline{v}_0(\underline{\lambda}) \quad \text{and} \quad \underline{\varphi}(\tilde{b}) = \underline{\varphi}_0(\underline{\lambda}). \qquad (69) $$

Then when $\varepsilon$ tends to zero, using Eq. (33), forces and moments can be shown to be solutions of the problem:

$$ \forall (\underline{v}_0, \underline{\varphi}_0), \quad \int_{\Lambda} \sum_{k \in \mathcal{R}} [\underline{f}_0^k \cdot (\underline{\nabla} \underline{v}_0 \cdot \underline{\delta}^k) + \underline{m}_0^k \cdot (\underline{\nabla} \underline{\varphi}_0 \cdot \underline{\delta}^k) - \underline{L}^k \wedge \underline{f}_0^k \cdot \underline{\varphi}_0] d\underline{\lambda} = 0 \qquad (70) $$

which, using the property $(\underline{a} \cdot \underline{b}) \cdot \underline{c} = \underline{b} \otimes \underline{c} : \underline{a}$, can also be written as

$$ \forall (\underline{v}_0, \underline{\varphi}_0), \quad \int_{\Lambda} [\underline{S}_0 : \underline{\nabla} \underline{v}_0 + \underline{M}_0 : \underline{\nabla} \underline{\varphi}_0 + \underline{Z}_0 \cdot \underline{\varphi}_0] d\underline{\lambda} = 0. \qquad (71) $$

In this last formulation, the following stress tensors, defined in the parametric space $\Lambda$, were introduced:

$$ \underline{S}_0 = \sum_{k \in \mathcal{R}} \underline{\delta}^k \otimes \underline{f}_0^k, \qquad (72) $$

$$ \underline{M}_0 = \sum_{k \in \mathcal{R}} \underline{\delta}^k \otimes \underline{m}_0^k, \qquad (73) $$

$$ \underline{Z}_0 = \sum_{k \in \mathcal{R}} \underline{f}_0^k \wedge \underline{L}^k. \qquad (74) $$

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 $\underline{\lambda} \mapsto \underline{p}_0(\underline{\lambda})$ in formulation (71). According to definition (8), the relation between the gradient of a vectorial field $\underline{u}$ with respect to $\underline{\lambda}$, denoted $\underline{\nabla} \underline{u}$, and its gradient with respect to $\underline{x} = \underline{p}_0$, denoted $\underline{\nabla}_x \underline{u}$, is given by

$$ \underline{\nabla} \underline{u} = \underline{\nabla}_x \underline{u} \cdot \underline{G}. \qquad (75) $$

The macroscopic equilibrium of the fiber net in the physical space then reads

$$ \forall (\underline{v}_0, \underline{\varphi}_0), \quad \int_{\Omega} [\underline{\sigma}_0 : \underline{\nabla}_x \underline{v}_0 + \underline{\kappa}_0 : \underline{\nabla}_x \underline{\varphi}_0 + \underline{\zeta}_0 \cdot \underline{\varphi}_0] dx = 0. \qquad (76) $$
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:

$$ 
\begin{aligned}
\underline{\underline{\mathbf{div}}} \underline{\underline{\underline{\boldsymbol{\sigma}}}}_0 &= \underline{\underline{0}}, \\
\underline{\underline{\underline{\mathbf{div}}}} \underline{\underline{\underline{\boldsymbol{\kappa}}}}_0 - \underline{\underline{\underline{\boldsymbol{\zeta}}}}_0 &= \underline{\underline{0}}.
\end{aligned}
\qquad (77) $$

Therefore, the state of stresses inside the fiber net is defined by the three following tensors:

$$ 
\begin{aligned}
\underline{\underline{\boldsymbol{\sigma}}}_0 &= g^{-1} \underline{\underline{\underline{\boldsymbol{G}}}} \cdot \underline{\underline{\underline{\boldsymbol{S}}}}_0: \text{ stress tensor,} \\
\underline{\underline{\underline{\boldsymbol{\kappa}}}}_0 &= g^{-1} \underline{\underline{\underline{\boldsymbol{G}}}} \cdot \underline{\underline{\underline{\boldsymbol{M}}}}_0: \text{ couple stress tensor,} \\
\underline{\underline{\underline{\boldsymbol{\zeta}}}}_0 &= g^{-1} \underline{\underline{\underline{\boldsymbol{Z}}}}_0: \text{ micro-stresses vector.}
\end{aligned}
\qquad (78) $$

In accordance with the general theory of Cosserat media exposed in Section 4.1, the macroscopic stress tensor $\underline{\underline{\sigma}}_0$ is non symmetric, and its antisymmetric part is directly linked to $\underline{\underline{\zeta}}_0$ by the relation

$$ \underline{\underline{\zeta}}_0 = \underline{\underline{e}} : \underline{\underline{\underline{\sigma}}}}_0 = \underline{\underline{e}} : \underline{\underline{\underline{\sigma}}}}_0^A. \qquad (79) $$

This general property of Cosserat media can easily be checked in the case of the fiber 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 first order interaction forces and moments. The stiffness of the fiber net will therefore be 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 fibers may be neglected as well as local moments.

## 5. Constitutive equations of the equivalent continuum

As Eqs. (72)–(74) show, the determination of the macroscopic state of stress $(\underline{\underline{\sigma}}_0, \underline{\underline{\kappa}}_0, \underline{\underline{\zeta}}_0)$ in the fiber net requires the determination of forces and moments of order $\varepsilon^0, \underline{\underline{f}}_0^k$ and $\underline{\underline{m}}_0^k$. Nevertheless, with the chosen local interaction relations, those quantities directly depend on the value of kinematic variables $\underline{\underline{v}}_1^b, \underline{\underline{\phi}}_0^b$ and $\underline{\underline{\phi}}_1^b$, as well as on the macroscopic velocity gradient $\underline{\underline{\nabla}} \underline{\underline{v}}_0$. To fully determine the constitutive equations corresponding to the three cases defined 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 $\beta^k$ with respect to translational one $\mu^k$, characterized by the parameter $q$ (see Section 2.4).
5.1. Case 1: $q=0$

According to results (43) and (53), order zero interaction laws become

$$ \underline{f}_0^k = \mu_0^k \|\underline{\Delta v}_1^k\|^{m-1} \underline{\Delta v}_1^k, \quad (80) $$

$$ \underline{m}_0^k = \beta_0^k \|\underline{\Delta \phi}_1^k\|^{m-1} \underline{\Delta \phi}_1^k \quad (81) $$

with

$$ \underline{\Delta v}_1^k = \underline{v}_1^{\text{E}(k)} - \underline{v}_1^{\text{I}(k)} + \underline{\phi}_0 \wedge (\xi^k \underline{\epsilon}^{\text{E}(k)} - \xi^k \underline{\epsilon}^{\text{I}(k)}) + \underline{\nabla} \underline{v}_0 \cdot \underline{\delta}^k, \quad (82) $$

$$ \underline{\Delta \phi}_1^k = \underline{\phi}_1^{\text{E}(k)} - \underline{\phi}_1^{\text{I}(k)} + \underline{\nabla} \underline{\phi}_0 \cdot \underline{\delta}^k. \quad (83) $$

Carrying into the general equilibrium formulation (30) the virtual functions

$$ \underline{v}(\tilde{b}) = \varepsilon \alpha_1(\underline{\lambda}) \underline{v}^b \quad \text{and} \quad \underline{\varphi}(\tilde{b}) = \varepsilon \alpha_2(\underline{\lambda}) \underline{\varphi}^b \quad (84) $$

and making $\varepsilon$ tend to zero, one gets

$$ \begin{gathered} \forall(\alpha_1, \underline{v}^b), \quad \forall(\alpha_2, \underline{\varphi}^b), \\ \int_{\Lambda} \sum_{k \in \mathcal{C}_R} \underline{f}_0^k \cdot (\underline{v}^{\text{E}(k)} - \underline{v}^{\text{I}(k)}) \alpha_1 \, d\lambda + \int_{\Lambda} \sum_{k \in \mathcal{C}_R} \underline{m}_0^k \cdot (\underline{\varphi}^{\text{E}(k)} - \underline{\varphi}^{\text{I}(k)}) \alpha_2 \, d\lambda = 0. \quad (85) \end{gathered} $$

This relation being satisfied for all fields $\alpha_1$ and $\alpha_2$, forces $\underline{f}_0^k$ and moments $\underline{m}_0^k$ are solutions of the problems:

$$ \forall \underline{v}^b, \quad \sum_{k \in \mathcal{C}_R} \underline{f}_0^k \cdot (\underline{v}^{\text{E}(k)} - \underline{v}^{\text{I}(k)}) = 0, \quad (86) $$

$$ \forall \underline{\varphi}^b, \quad \sum_{k \in \mathcal{C}_R} \underline{m}_0^k \cdot (\underline{\varphi}^{\text{E}(k)} - \underline{\varphi}^{\text{I}(k)}) = 0. \quad (87) $$

Formulations (86) and (87) are respectively force and moment order zero self-equilibrium equations. They are strictly equivalent to the following non-linear systems:

$$ \forall b \in \mathcal{B}_R, \quad \sum_{k \in \mathcal{C}(b)} \underline{f}_0^k(\underline{v}_1^{\text{E}(k)}, \underline{v}_1^{\text{I}(k)}, \underline{\phi}_0, \underline{\nabla} v_0) = \underline{0}, \quad (88) $$

$$ \forall b \in \mathcal{B}_R, \quad \sum_{k \in \mathcal{C}(b)} \underline{m}_0^k(\underline{\phi}_0^{\text{E}(k)}, \underline{\phi}_1^{\text{I}(k)}, \underline{\nabla} \phi_0) = \underline{0}. \quad (89) $$

Constitutive equations of the equivalent continuum can then be calculated considering as given the macroscopic fields $\underline{\nabla} \underline{v}_0$, $\underline{\nabla} \underline{\phi}_0$ and $\underline{\phi}_0$. The computation gives velocities $\underline{v}_1^b$ and $\underline{\phi}_1^b$ as functions of $\underline{\nabla} \underline{v}_0$, $\underline{\nabla} \underline{\phi}_0$ and $\underline{\phi}_0$, then Eqs. (72)-(74) and (78) give $\underline{\underline{\sigma}}_0$, $\underline{\underline{\kappa}}_0$ and $\underline{\underline{\zeta}}_0$ in terms of those fields, 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 fields 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:

$$ 
\begin{aligned}
\underline{\underline{\sigma}}_0 &= \underline{\underline{\sigma}}_0(\underline{\underline{\nabla}}\underline{\underline{v}}_0, \underline{\underline{\phi}}_0), \\
\underline{\underline{\kappa}}_0 &= \underline{\underline{\kappa}}_0(\underline{\underline{\nabla}}\underline{\underline{\phi}}_0), \\
\underline{\underline{\zeta}}_0 &= \underline{\underline{\epsilon}} : \underline{\underline{\sigma}}_0^A = \underline{\underline{\zeta}}_0(\underline{\underline{\nabla}}\underline{\underline{v}}_0, \underline{\underline{\phi}}_0)
\end{aligned} 
\quad (90) $$

which is consistent with the general formulation (65) obtained by continuous media theory (see Section 4.1).

## 5.2. Case 2: $q=1$

In this case, as shown in Section 3.4 order zero moments are null. The macroscopic couple-stress tensor $\underline{\underline{\kappa}}_0$ is therefore automatically null, and from Eqs. (77), the micro-stresses vector $\underline{\underline{\zeta}}_0$ is also null. Property (79) then immediately causes the anti-symmetric part of the macroscopic stress tensor to be null. Thus, the equilibrium of the equivalent continuum does not imply any local moment and $\underline{\underline{\sigma}}_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 $\underline{\underline{f}}_0^k$. Thanks to property (53), their expression is the same as in case 1; they are defined by Eq. (80), with $\underline{\underline{\Delta}}\underline{\underline{v}}_1^k$ given by Eq. (82). Then, adopting the same technique as in case 1, order zero forces $\underline{\underline{f}}_0^k$ can be proved to solve the self-equilibrium (86), which is strictly equivalent to the following non-linear system:

$$ \forall b \in \mathcal{B}_R, \quad \sum_{k \in \mathcal{C}(b)} \underline{\underline{f}}_0^k(\underline{\underline{v}}_1^{\mathrm{E}(k)}, \underline{\underline{v}}_1^{\mathrm{I}(k)}, \underline{\underline{\phi}}_0, \underline{\underline{\nabla}}\underline{\underline{v}}_0) = \underline{\underline{0}}. \quad (91) $$

Furthermore, property (79) implies an additional relation on forces $\underline{\underline{f}}_0^k$ which reads

$$ \sum_{k \in \mathcal{C}(b)} \underline{\underline{f}}_0^k(\underline{\underline{v}}_1^{\mathrm{E}(k)}, \underline{\underline{v}}_1^{\mathrm{I}(k)}, \underline{\underline{\phi}}_0, \underline{\underline{\nabla}}\underline{\underline{v}}_0) \wedge (\underline{\underline{G}} \cdot \underline{\underline{\delta}}^k) = \underline{\underline{0}}. \quad (92) $$

This last relation shows that $\underline{\underline{\phi}}_0$ can theoretically be expressed as a function of $\underline{\underline{\nabla}}\underline{\underline{v}}_0$. One then notices that Eqs. (91) and (92) enable $\underline{\underline{v}}_1^b$ and $\underline{\underline{\phi}}_0$ to be calculated in terms of the macroscopic field $\underline{\underline{\nabla}}\underline{\underline{v}}_0$. Here again, an explicit determination of the macroscopic stress tensor $\underline{\underline{\sigma}}_0$ can be achieved by the numerical solution of the non-linear system formed by both equations.

Finally, in this case, the fiber 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

$$ \underline{\underline{\mathbf{d}}\mathbf{v}} \underline{\underline{\sigma}}_0 = \underline{\underline{0}} \quad (93) $$
with a constitutive equation of the type

$$ \underline{\underline{\sigma}}_0 = \underline{\underline{\sigma}}_0(\underline{\underline{\nabla}}\underline{\underline{v}}_0, \underline{\underline{\phi}}_0(\underline{\underline{\nabla}}\underline{\underline{v}}_0)), \quad (94) $$

where $\underline{\underline{\sigma}}_0$ is a symmetric tensor. The relation between $\underline{\underline{\phi}}_0$ and $\underline{\underline{\nabla}}\underline{\underline{v}}_0$, is obtained in an implicit way by ensuring the symmetry condition.

Furthermore, in accordance with the general theory, because of the symmetry of $\underline{\underline{\sigma}}_0$, only the symmetric part of the macroscopic velocity gradient $\underline{\underline{\nabla}}\underline{\underline{v}}_0$ contributes to the total dissipated mechanical power. This causes $\underline{\underline{\sigma}}_0$ to depend only on the macroscopic strain rate $\underline{\underline{D}}_0$, defined as

$$ \underline{\underline{D}}_0 = \frac{1}{2}(\underline{\underline{\nabla}}\underline{\underline{v}}_0 + \iota\underline{\underline{\nabla}}\underline{\underline{v}}_0). \quad (95) $$

So in Eq. (94) $\underline{\underline{\nabla}}\underline{\underline{v}}_0$ may be replaced by $\underline{\underline{D}}_0$ and the equivalent continuum exhibits a general fluid-like behavior.

### 5.3. Case 3: $q = 1 + m$

In this case, moments of order 0 are immediately null, which causes $\underline{\underline{\kappa}}_0$ and $\underline{\underline{\zeta}}_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 defined by the single tensor $\underline{\underline{\sigma}}_0$. Its determination, according to Eq. (72), requires the determination of order zero interaction forces $\underline{\underline{f}}_0^k$.

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 $\underline{\underline{\phi}}_0^b$. Those variables now depend on the considered bar in the reference cell.

Forces $\underline{\underline{f}}_0^k$ are therefore defined by

$$ \underline{\underline{f}}_0^k = \mu_0^k \|\underline{\underline{\nabla}}\underline{\underline{v}}_1^k\|^{m-1} \underline{\underline{\nabla}}\underline{\underline{v}}_1^k \quad (96) $$

with

$$ \underline{\underline{\nabla}}\underline{\underline{v}}_1^k = \underline{\underline{v}}_1^{\mathrm{E}(k)} - \underline{\underline{v}}_1^{\mathrm{I}(k)} + \zeta_0^k \underline{\underline{\phi}}_0^{\mathrm{E}(k)} \wedge \underline{\underline{\epsilon}}^{\mathrm{E}(k)} - \zeta_0^k \underline{\underline{\phi}}_0^{\mathrm{I}(k)} \wedge \underline{\underline{\epsilon}}^{\mathrm{I}(k)} + \underline{\underline{\nabla}}\underline{\underline{v}}_0 \cdot \underline{\underline{\delta}}^k. \quad (97) $$

They solve Eq. (86), which is equivalent to the non-linear-system:

$$ \forall b \in \mathcal{B}_R, \quad \sum_{k \in \mathcal{C}(b)} \underline{\underline{f}}_0^k(\underline{\underline{v}}_1^{\mathrm{E}(k)}, \underline{\underline{v}}_1^{\mathrm{I}(k)}, \underline{\underline{\phi}}_0^{\mathrm{E}(k)}, \underline{\underline{\phi}}_0^{\mathrm{I}(k)}, \underline{\underline{\nabla}}\underline{\underline{v}}_0) = \underline{\underline{0}}, \quad (98) $$

the unknowns of which are now $\underline{\underline{v}}_1^b$ and $\underline{\underline{\phi}}_0^b$. It therefore brings 2 equations per bar whereas 3 unknowns per bar are to be determined, $\underline{\underline{\nabla}}\underline{\underline{v}}_0$ being considered as data.

It is to be noted that in this case, order one moments $\underline{\underline{m}}_1^k$, defined by

$$ \underline{\underline{m}}_1^k = \beta_0^k \|\underline{\underline{\phi}}_0^{\mathrm{E}(k)} - \underline{\underline{\phi}}_0^{\mathrm{I}(k)}\|^{m-1} (\underline{\underline{\phi}}_0^{\mathrm{E}(k)} - \underline{\underline{\phi}}_0^{\mathrm{I}(k)}) + \zeta_0^k \underline{\underline{\epsilon}}^{\mathrm{I}(k)} \wedge \underline{\underline{f}}_0^k \quad (99) $$

imply the same kinematic variables as $\underline{\underline{f}}_0^k$.

The missing equations thus can be obtained by taking in Eq. (30) the virtual functions $\underline{\underline{v}}(\tilde{b})=0$ and $\underline{\underline{\phi}}(\tilde{b})=\alpha_2(\underline{\underline{\lambda}})\underline{\underline{\phi}}_0^b$ leading to the moments order one self-equilibrium following
formulation:

$$ \forall \underline{\boldsymbol{\varphi}}^b, \sum_{k \in \mathcal{C}_R} \underline{\boldsymbol{m}}_1^k \cdot (\underline{\boldsymbol{\varphi}}^{\mathrm{E}(k)} - \underline{\boldsymbol{\varphi}}^{\mathrm{I}(k)}) + \underline{\boldsymbol{L}}^k \wedge \underline{\boldsymbol{f}}_0^k \cdot \underline{\boldsymbol{\varphi}}^{\mathrm{E}(k)} = 0. \quad (100) $$

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

$$ \forall b \in \mathcal{B}_R, \sum_{k \in \mathcal{C}(b)} \underline{\boldsymbol{m}}_1^k (\underline{\boldsymbol{v}}_1^{\mathrm{E}(k)}, \underline{\boldsymbol{v}}_1^{\mathrm{I}(k)}, \underline{\boldsymbol{\phi}}_0^{\mathrm{E}(k)}, \underline{\boldsymbol{\phi}}_0^{\mathrm{I}(k)}, \underline{\boldsymbol{\nabla}} \underline{\boldsymbol{v}}_0) = \underline{\boldsymbol{0}}. \quad (101) $$

Vectors $\underline{\boldsymbol{v}}_1^b$, and $\underline{\boldsymbol{\phi}}_0^b$ can therefore be computed by the simultaneous solving of Eqs. (98) and (101) in terms of $\underline{\boldsymbol{\nabla}} \underline{\boldsymbol{v}}_0$, what will provide $\underline{\boldsymbol{f}}_0^k$, and then $\underline{\boldsymbol{\sigma}}_0$ in terms of this macroscopic velocity gradient.

Finally, as in case 2, the fiber 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:

$$ \underline{\boldsymbol{\sigma}}_0 = \underline{\boldsymbol{\sigma}}_0(\underline{\boldsymbol{D}}_0), \quad (102) $$

where $\underline{\boldsymbol{\sigma}}_0$ is a symmetric tensor and $\underline{\boldsymbol{D}}_0$ is the macroscopic strain rate tensor. Cases 2 and 3 finally happen to lead to the same type of equivalent continuous medium, even if the calculation of their constitutive equations leads to somewhat different 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.

## 6. 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 first 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 $\underline{\boldsymbol{v}}_1^b$ and $\underline{\boldsymbol{\phi}}_0^b$ in terms of $\underline{\boldsymbol{\nabla}} \underline{\boldsymbol{v}}_0$. Such a system can be summarized in the following way:

$$ \underline{\boldsymbol{F}}(\underline{\boldsymbol{X}}, \underline{\boldsymbol{\nabla}} \underline{\boldsymbol{v}}_0) = \underline{\boldsymbol{0}}, \quad (103) $$

where $\underline{\boldsymbol{F}}$ is a vector that contains the force and moment equilibrium equations of each bar of the reference cell, and $\underline{\boldsymbol{X}}$ is a vector containing the kinematic unknowns $\underline{\boldsymbol{v}}_1^b$ and $\underline{\boldsymbol{\phi}}_0^b$ relative to each bar. $\underline{\boldsymbol{F}}$ is a block vector where blocks of components $[3b - 2:3b]$, relative to bar $b$, are denoted $\underline{\boldsymbol{F}}^b$ and read

$$ \underline{\boldsymbol{F}}^b = \sum_{k \in \mathcal{C}(b)} \left\{ \begin{array}{c} \underline{\boldsymbol{f}}_0^k \\ \underline{\boldsymbol{m}}_1^k \end{array} \right\} \quad (104) $$
whereas $\mathbf{X}$ may be assembled in blocks $\mathbf{X}^b$ such as

$$ \underline{\mathbf{X}}^b = \left\{ \begin{matrix} \underline{\mathbf{v}}_1^b \\ \underline{\boldsymbol{\phi}}_0^b \end{matrix} \right\}. \qquad (105) $$

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 sufficiently high for every bar to be connected with some of its neighbors.

Let us now consider that $\mathbf{X}$ is the solution of Eq. (103) and that $\hat{\mathbf{X}}$ is the solution of

$$ \underline{\mathbf{F}}(\hat{\mathbf{X}}, \alpha \underline{\nabla} \mathbf{v}_0) = 0, \qquad (106) $$

where $\alpha$ is a non-null real scalar, and study the relation between $\mathbf{X}$ and $\hat{\mathbf{X}}$.

If $\underline{\mathbf{v}}_1^b$, $\underline{\boldsymbol{\phi}}_0^b$ and $\underline{\hat{\mathbf{v}}}_1^b$, $\underline{\hat{\boldsymbol{\phi}}}_0^b$ are the respective solutions of problems (103) and (106), forces and moments solutions of Eq. (106) are then

$$ \underline{\hat{\mathbf{f}}}_0^k = \mu_0^k \|\underline{\Delta \hat{\mathbf{v}}}_1^k\|^{m-1} \underline{\Delta \hat{\mathbf{v}}}_1^k, \qquad (107) $$

$$ \underline{\hat{\mathbf{m}}}_1^k = \beta_0^k \|\underline{\Delta \hat{\boldsymbol{\phi}}}_0^k\|^{m-1} \underline{\Delta \hat{\boldsymbol{\phi}}}_0^k + \xi^k \underline{\boldsymbol{e}}^{l(\tilde{k})} \wedge \underline{\hat{\mathbf{f}}}_0^k \qquad (108) $$

with

$$ \underline{\Delta \hat{\mathbf{v}}}_1^k = \underline{\hat{\mathbf{v}}}_1^{\mathrm{E}(k)} - \underline{\hat{\mathbf{v}}}_1^{\mathrm{I}(k)} + \xi^k \underline{\hat{\boldsymbol{\phi}}}_0^{\mathrm{E}(k)} \wedge \underline{\boldsymbol{e}}^{\mathrm{E}(k)} - \xi^k \underline{\hat{\boldsymbol{\phi}}}_0^{\mathrm{I}(k)} \wedge \underline{\boldsymbol{e}}^{\mathrm{I}(k)} + \alpha \underline{\nabla} \underline{\mathbf{v}}_0 \cdot \underline{\boldsymbol{\delta}}^k, \qquad (109) $$

$$ \underline{\Delta \hat{\boldsymbol{\phi}}}_0^k = \underline{\hat{\boldsymbol{\phi}}}_0^{\mathrm{E}(k)} - \underline{\hat{\boldsymbol{\phi}}}_0^{\mathrm{I}(k)}. \qquad (110) $$

They can therefore be rewritten as

$$ \underline{\hat{\mathbf{f}}}_0^k = \alpha^m (\mu_0^k \|\alpha^{-1} \underline{\Delta \hat{\mathbf{v}}}_1^k\|^{m-1} \alpha^{-1} \underline{\Delta \hat{\mathbf{v}}}_1^k), \qquad (111) $$

$$ \underline{\hat{\mathbf{m}}}_1^k = \alpha^m (\beta_0^k \|\alpha^{-1} \underline{\Delta \hat{\boldsymbol{\phi}}}_0^k\|^{m-1} \alpha^{-1} \underline{\Delta \hat{\boldsymbol{\phi}}}_0^k + \xi^k \underline{\boldsymbol{e}}^{l(\tilde{k})} \wedge \mu_0^k \|\alpha^{-1} \underline{\Delta \hat{\mathbf{v}}}_1^k\|^{m-1} \alpha^{-1} \underline{\Delta \hat{\mathbf{v}}}_1^k). \qquad (112) $$

Carrying Eqs. (111) and (112) in Eq. (106) and according to the expression of vector $\underline{\mathbf{F}}$, one checks that $\hat{\mathbf{X}}$ is solution of the problem:

$$ \alpha^m \underline{\mathbf{F}}(\alpha^{-1} \underline{\hat{\mathbf{X}}}, \underline{\nabla} \underline{\mathbf{v}}_0) = 0 \Leftrightarrow \underline{\mathbf{F}}(\alpha^{-1} \underline{\hat{\mathbf{X}}}, \underline{\nabla} \underline{\mathbf{v}}_0) = 0. \qquad (113) $$

The uniqueness of the solution then immediately causes

$$ \underline{\hat{\mathbf{X}}} = \alpha \underline{\mathbf{X}}. \qquad (114) $$

Subsequently, forces and moments resulting from problem (106) are such that

$$ \underline{\hat{\mathbf{f}}}_0^k = \alpha^m \underline{\mathbf{f}}_0^k \quad \text{and} \quad \underline{\hat{\mathbf{m}}}_1^k = \alpha^m \underline{\mathbf{m}}_1^k, \qquad (115) $$

where $\underline{\mathbf{f}}_0^k$ and $\underline{\mathbf{m}}_1^k$ result from the initial problem (103).
Now forming the macroscopic stress tensors $\underline{\underline{\sigma}}_0$ and $\underline{\underline{\sigma}}_0$ defined by Eqs. (78) and (72) corresponding to both problems finally leads to

$$ \underline{\underline{\sigma}}_0 = \alpha^m \underline{\underline{\sigma}}_0. \qquad (116) $$

This result therefore enables one to deduce the following important property:

$$ \forall \alpha \in \mathbb{R}, \quad \underline{\underline{\sigma}}_0(\alpha \underline{\underline{D}}_0) = \alpha^m \underline{\underline{\sigma}}_0(\underline{\underline{D}}_0). \qquad (117) $$

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 fluid, with a strain rate sensitivity equal to the strain rate sensitivity postulated at the level of interactions between fibers. Actually, if one defines a norm $||\cdot||_{\text{eq}}$ in the space of second order tensors, property (117) enables us to write the following relation, characteristic of power law fluids:

$$ \underline{\underline{\sigma}}_0(\underline{\underline{D}}) = \underline{\underline{\sigma}}_0 \left( \frac{||\underline{\underline{D}}||_{\text{eq}} \underline{\underline{D}}}{||\underline{\underline{D}}||_{\text{eq}}} \right) = ||\underline{\underline{D}}||_{\text{eq}}^m \underline{\underline{\sigma}}_0 \left( \frac{\underline{\underline{D}}}{||\underline{\underline{D}}||_{\text{eq}}} \right). \qquad (118) $$

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 fields $\underline{\underline{\nabla}} v_0$ and $\underline{\underline{\phi}}_0$ (directly linked to $\underline{\underline{\nabla}} \underline{\underline{\phi}}_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:

$$ \begin{aligned} \forall \alpha \in \mathbb{R}, \quad & \underline{\underline{\sigma}}_0(\alpha \underline{\underline{\nabla}} v_0, \alpha \underline{\underline{\phi}}_0) = \alpha^m \underline{\underline{\sigma}}_0(\underline{\underline{\nabla}} v_0, \underline{\underline{\phi}}_0), \\ & \underline{\underline{\zeta}}_0(\alpha \underline{\underline{\nabla}} v_0, \alpha \underline{\underline{\phi}}_0) = \alpha^m \underline{\underline{\zeta}}_0(\underline{\underline{\nabla}} v_0, \underline{\underline{\phi}}_0), \\ & \underline{\underline{\kappa}}_0(\alpha \underline{\underline{\nabla}} \underline{\underline{\phi}}_0) = \alpha^m \underline{\underline{\kappa}}_0(\underline{\underline{\nabla}} \underline{\underline{\phi}}_0). \end{aligned} \qquad (119) $$

The equivalent Cosserat medium exhibits a degree $m$ homogeneity property in terms of the pairs $(\underline{\underline{\nabla}} v_0, \underline{\underline{\phi}}_0)$.

## 7. Conclusions

This theoretical work on the behavior of a net of rigid fibers linked by punctual power law fiber-fiber interactions shows several interesting results.

If the scale separation assumption (3) is satisfied 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 defined 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 different fiber-fiber interaction laws leads to 2 main different 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 fiber 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 fiber-fiber interactions. Such a case would probably be relevant for almost rigid interactions ($m \approx 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, defined by a single symmetric stress tensor.

Constitutive equations of the fiber 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 fiber-fiber 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 fluid with the same strain rate sensitivity $m$ as the one postulated at the scale of fibers. 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 fields.

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 specific net of fibers. It requires almost no restrictive physical assumptions and enables an easy, computer time efficient, analysis of the behavior for fiber nets eventually including a great number of fibers, which is a necessary feature for the study of most fiber-reinforced fluids.

As is visible in the above theoretical exposition, many extensions to this work can be envisaged. Richer fiber-fiber interaction laws could first be introduced, as, for example, the case of Carreau type or viscoelastic interactions, or dry friction between fibers. The great adaptability of the method would also enable one to account for the flexibility of fibers, considering them, for example, as elastic beams.

## Appendix A

As shown in Section 4, in the general case, the state of stress of the fiber net is defined by the three tensors $\underline{\underline{\sigma}}_0$, $\underline{\underline{\kappa}}_0$ and $\underline{\underline{\zeta}}_0$, the vector $\underline{\underline{\zeta}}_0$ being given by the relation

$$g\underline{\underline{\zeta}}_0 = \sum_{k \in \mathcal{C}_R} \underline{\underline{f}}_0^k \wedge \underline{\underline{L}}^k = \sum_{k \in \mathcal{C}_R} \underline{\underline{f}}_0^k \wedge (\underline{\underline{G}} \cdot \underline{\underline{\delta}}^k + \underline{\underline{p}}_1^c - \underline{\underline{p}}_1^b)$$
according to definition (32) of $\underline{L}^k$. Let us now multiply this last equation by a constant virtual field $\phi_0$. We get

$$ \forall \underline{\phi}_0, \quad g\underline{\zeta}_0 \cdot \underline{\phi}_0 = \sum_{k \in \mathcal{C}_R} \underline{f}_0^k \wedge (\underline{\underline{G}} \cdot \underline{\delta}^k) \cdot \underline{\phi}_0 + \sum_{k \in \mathcal{C}_R} \underline{f}_0^k \cdot (\underline{\underline{p}}_1^c \wedge \underline{\phi}_0 - \underline{\underline{p}}_1^b \wedge \underline{\phi}_0). $$

Terms like $\underline{\underline{p}}_1^b \wedge \underline{\phi}_0$ can be seen as a virtual field $\underline{v}^b$, so forces $\underline{\underline{f}}_0^k$ being solutions of the self-equilibrium equation (86), one obtains

$$ \forall \underline{\phi}_0, \quad g\underline{\zeta}_0 \cdot \underline{\phi}_0 = \sum_{k \in \mathcal{C}_R} \underline{f}_0^k \wedge (\underline{\underline{G}} \cdot \underline{\delta}^k) \cdot \underline{\phi}_0 $$

which gives an alternate definition of vector $\underline{\zeta}_0$ as

$$ \underline{\zeta}_0 = g^{-1} \sum_{k \in \mathcal{C}_R} \underline{f}_0^k \wedge (\underline{\underline{G}} \cdot \underline{\delta}^k). $$

This result then enables us to find the relation between $\underline{\underline{\sigma}}_0$ and $\underline{\underline{\zeta}}_0$, using the definition and properties of the permutation tensor $\underline{\underline{e}}$. $\underline{\underline{e}}$ is the tensor whose components in $\underline{\underline{e}}_i \otimes \underline{\underline{e}}_j \otimes \underline{\underline{e}}_k$ are

$$ e_{ijk} = \begin{cases} 1 & \text{if } ijk \text{ belongs to the set } \{123, 231, 312\}, \\ -1 & \text{if } ijk \text{ belongs to the set } \{132, 321, 213\}, \\ 0 & \text{in any other case.} \end{cases} $$

It is then easy to check that

$$ \underline{\underline{\zeta}}_0 = g^{-1} \sum_{k \in \mathcal{C}_R} (\underline{\underline{e}} \cdot (\underline{\underline{G}} \cdot \underline{\underline{\delta}}^k)) \cdot \underline{\underline{f}}_0^k = \underline{\underline{e}} : \sum_{k \in \mathcal{C}_R} g^{-1} (\underline{\underline{G}} \cdot \underline{\underline{\delta}}^k \otimes \underline{\underline{f}}_0^k) $$

which is equivalent to

$$ \underline{\underline{\zeta}}_0 = \underline{\underline{e}} : \underline{\underline{\sigma}}_0. $$

## References

Advani, S., 1994. Flow and rheology in polymer composites manufacturing. In: Composite Materials Series, Vol. 10. Elsevier, New York.

Bakhvalov, N.S., Panasenko, B., 1989. Homogenisation: averaging processes in periodic media—mathematical problems in the mechanics of composite materials. In: Mathematics and its Applications Soviet Series, Vol. 36. Kluwer Academic Publishers, Dordrecht-Boston-London.

Batchelor, G., 1971. The stress generated in a non-dilute suspension of elongated particles by pure straining motion. J. Fluid Mech. 46 (4), 813–829.

Bensoussan, A., Lions, J.-L., Papanicolaou, G., 1978. Asymptotic Analysis for Periodic Structures. North-Holland, Amsterdam.

Cosserat, E., Cosserat, F., 1909. Théorie Des Corps Déformables. Herman et fils edition, Paris.

Dumont, P., Orgéas, L., Le Corre, S., Favier, D., 2003. Anisotropic viscous behaviour of sheet molding compounds (SMC) during compression molding. Int. J. Plast. 19, 625–646.

Eringen, A.C., 1968. Mechanics of micromorphic continua. In: E. Kröner (ed), Mechanics of Generalized Continua, I.U.T.A.M. Symposium, Freudenstadt-Stuttgart, Germany, Springer, Berlin, pp. 18–35.

Fredrickson, G., Shaqfeh, E.S.G., 1989. Heat and mass transport in composites of aligned slender fibers. Phys. Fluids 1 (1), 3–20.
Germain, P., 1973. La méthode des puissances virtuelles en mécanique des milieux déformables. J. Méc. 12 (2), 235-274.

Gibson, A., Toll, S., 1999. Mechanics of the squeeze flow of planar fibre suspensions. J. Non-Newtonian Fluid Mech. 82, 1-24.

Jeffery, G.B., 1922. The motion of ellipsoidal particles immersed in a viscous fluid. Proc. Rheo. Soc. A(102), 161-179.

Le Corre, S., Orgéas, L., Favier, D., Tourabi, A., Maazouz, A., Venet, C., 2002. Shear and compression behaviour of sheet molding compounds. Composite Sci. Technol. 62, 571-577.

Moreau, G., 1996. Homogénéisation de Structures Discrètes En Elasticité et En Incrémental. Ph.D. Thesis, Institut National Polytechnique de Grenoble.

Moreau, G., Caillerie, D., 1995. Continuum modeling of lattice structures in large displacement framework. Dev. Comput. Tech. Struct. Eng. 20, 53-70.

Pradel, F., 1998. Homogénéisation Des Milieux Continus et Discrets Périodiques Orientés. Une Application Aux Mousses. Ph.D. Thesis, Ecole Nationale des Ponts et Chaussées, Paris, France.

Sanchez-Palencia, E., 1980. Non-Homogeneous Media and Vibration Theory. Lecture Notes in Physics. Springer, Berlin, Germany.

Servais, C., Manson, J.-A.E., Toll, S., 1999. Fiber-fiber interaction in concentrated suspensions: disperse fibers. J. Rheology 43 (4), 991-1004.

Shaqfeh, E.S.G., Fredrickson, G., 1990. The hydrodynamic stress in a suspension of rods. Phys. Fluids 2 (1), 7-24.

Toll, S., Manson, J.-A.E., 1994. Dynamics of a planar concentrated suspension with non-hydrodynamic interaction. J. Rheol. 38 (4), 985-997.

Tollenaere, H., Caillerie, D., 1998. Continuous modeling of lattice structures by homogenization. Adv. Eng. Software (special issue) 29 (7), 699-705.