Multiscale modeling of fiber reinforced materials via non-matching immersed methods

TL;DR

This paper introduces a multiscale modeling approach for fiber reinforced materials using non-matching immersed methods, enabling complex fiber geometries and couplings without regularity assumptions, validated against classical models.

Contribution

It proposes a distributed Lagrange multiplier framework for coupling fibers and matrix, allowing flexible geometries and intractable configurations, with rigorous analysis and validation.

Findings

The method handles complex fiber geometries without alignment constraints.

Well-posedness and convergence are established for continuous and discrete models.

Validation confirms accuracy against classical homogenization and rule of mixtures models.

Abstract

Fiber reinforced materials (FRMs) can be modeled as bi-phasic materials, where different constitutive behaviors are associated with different phases. The numerical study of FRMs through a full geometrical resolution of the two phases is often computationally infeasible, and therefore most works on the subject resort to homogenization theory, and exploit strong regularity assumptions on the fibers distribution. Both approaches fall short in intermediate regimes where lack of regularity does not justify a homogenized approach, and when the fiber geometry or their numerosity render the fully resolved problem numerically intractable. In this paper, we propose a distributed Lagrange multiplier approach, where the effect of the fibers is superimposed on a background isotropic material through an independent description of the fibers. The two phases are coupled through a constraint condition, opening the way for intricate fiber-bulk couplings as well as allowing complex geometries with no alignment requirements between the discretisation of the background elastic matrix and the fibers. We analyze both a full order coupling, where the elastic matrix is coupled with fibers that have a finite thickness, as well as a reduced order model, where the position of their centerline uniquely determines the fibers. Well posedness, existence, and uniquess of solutions are shown both for the continuous models, and for the finite element discretizations. We validate our approach against the models derived by the rule of mixtures, and by the Halpin-Tsai formulation.

Tables

Table 1. Table 1 : Fiber Parameters

Fiber length	0.6	0.4	0.3	0.25	0.225	0.2	0.18
Fiber radius	0.03	0.02	0.015	0.0125	0.01125	0.01	0.009
Number of Fibers	79	268	637	1100	1509	2149	2947

Equations

$$S[u]=\mathbb{C}\nabla u,$$

$$\mathbb{C}=\{\begin{array}{ll}{\mathbb{C}}_{\Omega }& \text{in }\hat{\Omega },\\ {\mathbb{C}}_f& \text{in }{\Omega }_f.\end{array}$$

$$-\mathrm{div}(\mathbb{C}u)$$

$$u$$

$$H^{1}_{0,D}(\Omega)^{d}\coloneqq\{v\in H^{1}(\Omega)^{d}\colon{\left.\kern-1.2ptv\vphantom{\big{|}}\right|_{D}}=0\},$$

$$V$$

$$(\mathbb{C}\nabla u,\nabla v{)}_{\Omega }=(b,v{)}_{\Omega }\forall v\in V.$$

$$(\mathbb{C}\nabla u,\nabla v{)}_{\Omega }$$

$$=({\mathbb{C}}_f\nabla u,\nabla v{)}_{{\Omega }_f}+\underset{\Omega }{\underbrace{({\mathbb{C}}_{\Omega }\nabla u,\nabla v{)}_{\hat{\Omega }}+({\mathbb{C}}_{\Omega }\nabla u,\nabla v{)}_{{\Omega }_f}}}-({\mathbb{C}}_{\Omega }\nabla u,\nabla v{)}_{{\Omega }_f}$$

$$=({\mathbb{C}}_{\Omega }\nabla u,\nabla v{)}_{\Omega }+(\delta {\mathbb{C}}_f\nabla u,\nabla v{)}_{{\Omega }_f},$$

$$({\mathbb{C}}_{\Omega }\nabla u,\nabla v{)}_{\Omega }+(\delta {\mathbb{C}}_f\nabla u,\nabla v{)}_{{\Omega }_f}=(b,v{)}_{\Omega }\forall v\in V.$$

$$B_i$$

$$B_e$$

$$W$$

$${\left.\kern-1.2ptu\vphantom{\big{|}}\right|_{\Omega_{f}}}=w.$$

$$u,w=\arg \underset{\begin{array}{c}(u,v)\in V\times W\\ u{|}_{{\Omega }_f}=w\end{array}}{\inf }\psi (u,w),$$

$$\psi (u,w)=\frac{1}{2}({\mathbb{C}}_{\Omega }\nabla u,\nabla u{)}_{\Omega }+\frac{1}{2}(\delta {\mathbb{C}}_f\nabla w,\nabla w{)}_{{\Omega }_f}-(b,u{)}_{\Omega }.$$

$$\langle q,{\left.\kern-1.2ptu\vphantom{\big{|}}\right|_{\Omega_{f}}}-w\rangle_{Q\times W}=0\qquad\forall q\in Q.$$

$$u,w,\lambda =\arg \underset{\begin{array}{c}u\in V\\ w\in W\end{array}}{\inf }\left(\arg \underset{\begin{array}{c}\lambda \in Q\end{array}}{\sup }\psi (u,w,\lambda )\right),$$

$$\psi(u,w,\lambda):=\frac{1}{2}(\mathbb{C}_{\Omega}\nabla u,\nabla u)_{\Omega}+\frac{1}{2}(\delta\mathbb{C}_{f}\nabla w,\nabla w)_{\Omega_{f}}+\langle\lambda,{\left.\kern-1.2ptu\vphantom{\big{|}}\right|_{\Omega_{f}}}-w\rangle-(b,u)_{\Omega}.$$

$$⟨D_u\psi ,v⟩+⟨D_w\psi ,y⟩+⟨D_{\lambda }\psi ,q⟩=0\forall v\in V,\forall y\in W,\forall q\in Q,$$

$$\displaystyle(\mathbb{C}_{\Omega}\nabla u,\nabla v)_{\Omega}+\langle\lambda,{\left.\kern-1.2ptv\vphantom{\big{|}}\right|_{\Omega_{f}}}\rangle$$

$$(\delta {\mathbb{C}}_f\nabla w,\nabla y{)}_{{\Omega }_f}-⟨\lambda ,y⟩$$

$$\displaystyle\langle q,{\left.\kern-1.2ptu\vphantom{\big{|}}\right|_{\Omega_{f}}}-w\rangle$$

$${\mathcal{K}}_{\Omega }u$$

$${\mathcal{K}}_{f}w$$

$$\mathcal{B}u$$

$${\mathcal{K}}_{\Omega }$$

$${\mathcal{K}}_f$$

$$\mathcal{B}$$

$$\mathcal{M}$$

$$\mathbb{V}$$

$$\Vert (u,w){\Vert }_{\mathbb{V}}^2$$

$$\begin{array}{rl}\mathbb{F}:& \mathbb{V}\times \mathbb{V}⟶\mathbb{R}\\ & (\mathbf{u},\mathbf{v})⟼⟨{\mathcal{K}}_{\Omega }u,v⟩+⟨{\mathcal{K}}_{f}w,y⟩=({\mathbb{C}}_{\Omega }\nabla u,\nabla v{)}_{\Omega }+(\delta {\mathbb{C}}_f\nabla w,\nabla y{)}_{{\Omega }_f},\end{array}$$

$$\displaystyle\begin{aligned} \mathbb{E}\colon&\mathbb{V}\times Q\longrightarrow\mathbb{R}\\ &(\boldsymbol{u},q)\longmapsto\langle\mathcal{B}u,q\rangle-\langle\mathcal{M}w,q\rangle=\langle q,{\left.\kern-1.2ptv\vphantom{\big{|}}\right|_{\Omega_{f}}}-w\rangle,\end{aligned}$$

$$\mathbb{F}(\mathbf{u},\mathbf{v})+\mathbb{E}(\mathbf{v},\lambda )$$

$$\mathbb{E}(\mathbf{u},q)$$

Full text

In this paper, we propose a distributed Lagrange multiplier approach, where the effect of the fibers is superimposed on a background isotropic material through an independent description of the fibers. The two phases are coupled through a constraint condition, opening the way for intricate fiber-bulk couplings as well as allowing complex geometries with no alignment requirements between the discretisation of the background elastic matrix and the fibers.

We analyze both a full order coupling, where the elastic matrix is coupled with fibers that have a finite thickness, as well as a reduced order model, where the position of their centerline uniquely determines the fibers. Well posedness, existence, and uniquess of solutions are shown both for the continuous models, and for the finite element discretizations. We validate our approach against the models derived by the rule of mixtures, and by the Halpin-Tsai formulation.

1 Introduction

Numerous engineering applications require the efficient solution of partial differential equations involving multiple, complex geometries on different phases; composite materials are the prototypical example of such problems. During the past fifty years, the interest in composite materials flourished multiple times; it began for their applications to new materials in multiple fields, such as aerospace engineering hashin1972theory , civil engineering mehta1986concrete , and materials science behzad2007finite ; mallick2007fiber .

During the nineties, the increasing importance of biomechanics in life sciences lead to the development of numerous models describing, e.g., arterial walls holzapfel2000new , soft tissues holzapfel2001biomechanics , and muscle fibers huijing1999muscle .

Recent years saw the rise of new application fields, such as the study of natural fiber composites pickering2016review , and engineering methods to accurately recover the three-dimensional structure of a material sample, e.g., iizuka2019development ; konopczynski2019fully .

From the first studies on composites, it has been clear that their properties are strongly dependent on their internal structure: the volume ratio between each component, the orientation, the shape, all contribute substantially to the material’s properties hashin1972theory ; hashin1983analysis . One of the most wide-spread and significant example of composites is that of Fiber Reinforced Materials (FRMs), where thin, elongated structures (the fibers) are immersed in an underlying isotropic material (the elastic matrix).

We may separate the approaches used to study FRMs into two broad groups: i) “homogenization methods”, which study a complex inhomogeneous body by approximating it with a fictitious homogeneous body that behaves globally in the same way zaoui2002continuum , and ii) “fully resolved” methods, which use separate geometrical and constitutive descriptions for the elastic matrix and the individual fibers.

As examples of analytical “homogenization methods” we recall the rule of mixtures (hashin1965elastic, ) and the empirical Halpin-Tsai equations halpin1976halpin , used to study a transversely isotropic unidirectional composite, where fibers are uniformly distributed and share the same orientation. The development of homogenization theory led, in recent years, to more complex models, e.g., holzapfel2001viscoelastic ; barrage2019modelling .

More intricate homogeneization approaches rely on numerical methods to provide a “cell” behaviour, which is then replicated using periodicity, using, e.g., the Finite Element Method adams1967transverse ; nakamura1993effects ; lu20143d , Fourier transforms moulinec1998numerical ; moulinec2003comparison ; eloh2019development , or Stochastic Methods lucas2015stochastic .

The fundamental limit of all “homogenization methods” approaches is the impossibility of adapting them to study composites with little regularity. In these cases, the different phases are typically modeled separately, as a continuum. This approach began with Pipkin Pipkin1968 on two dimensional membranes, and was then expanded to three dimensional examples by others (see, e.g., Kyriacou1996 for a detailed bibliography).

Fully resolved methods allow richer structures, but require a high numerical resolution, especially when material phases have different scales. The complex meshing and coupling often result in an unbearable computational cost, limiting the use of these methods.

The purpose of this paper is to introduce an approach which is fit for materials that have intermediate properties, i.e., they possess no particular regularity, and are made by a relatively high number of fiber components. Similarly to radtke2010computational , our model is inspired by the Immersed Boundary Method (IBM) peskin2002 , and by its variational counterparts boffi2008hyper ; Heltai-2008-a ; HeltaiCostanzo-2012-a ; RoyHeltaiCostanzo-2015-a ; HeltaiRotundo-2018-a , where the elastic matrix and the fibers are modeled independently, and coupled through a non-slip condition. The model we present can be interepreted as a variation of the embedded reinforced method (see goudarzi2019discrete and the references therein), where we aim at providing an efficient numerical method for FRMs that allows the modeling of complex networks of fibers, where one may also be interested in the elastic properties of single fibers, without requiring the resolution of the single fibers in the background elastic matrix. From the computational point of view, this approach allows the use of two independent discretizations: one describing the fibers, and one describing the whole domain, i.e., both the elastic matrix and the fibers. A distributed Lagrange multiplier is used to couple the independent grids, following the same spirit of the finite element immersed boundary method boffi2003finite ; boffi2008hyper , separating the Cauchy stress of the whole material into a background uniform behavior and into an excess elastic behavior on the fibers.

Section 2 introduces the classical fully resolved model of a collection of fibers immersed in an elastic matrix. For simplicity, we do not include dissipative terms, and restrict our study to linearly elastic materials. The problem is then reformulated exploiting classical results of mixed methods (see Chapter 4 of boffi2013mixed ), following ideas similar to those found in boffi2017fictitious , proving that both the continuous and discrete formulations we propose are well-posed with a unique solution.

The use of a full three dimensional model for the fibers still results in high computational costs; the obvious simplification would be to approximate the fibers with one-dimensional structures. This approach is non-trivial because the theoretical solution of the fully resolved variational problem does not posses enough Sobolev regularity to allow the restriction of a three-dimensional field to a one-dimensional domain. A possible workaround involves the use of weighted Sobolev spaces, combinded with graded meshes d2008coupling ; d2012finite . This approach remains unfeasible when the number of fibers is large. In Section 3, we propose and analyze an alternative strategy, where under some additional assumptions we construct a $3D-1D$ coupling that relies on local averaging techniques. A similar procedure is used in heltaimultiscale to model vascularized tissues. To conclude, we validate our thin fiber model in Section 4, and draw some conclusions in Section 5.

2 Three-dimensional model

Many bi-phasic materials present a relatively simple fiber structure but result in a very intricate elastic matrix. Consider, for example, Figure 1: constructing a discretization grid for the fibers themselves maybe simple enough, but building a fully resolved grid for the surrounding elastic matrix, in this case, may require eccessive resolution, and result in a computationally hard problem to solve. We wish to describe a new approach, where we substitute the complex mesh needed for the elastic matrix with a simple one describing the whole domain, and overlap the fiber structure independently with respect to the background grid, and couple the two systems via distributed Lagrange multipliers.

2.1 Problem formulation

As a model bi-phasic material, we consider a linearly elastic fiber reinforced material, in the quasi-static small strain regime. We consider continuous fibers with a perfect bond between the two phases (see (callister2007materials, , Chapter 16)), leading to a no-slip condition between the fibers and the elastic matrix. The extension to finite strain elasticity, dynamic problems, or non-perfect bonds do not present additional difficulties and will be the subject of future investigations.

To describe the composite, we use a connected, bounded, Lipschitz domain $\Omega\subset\mathbb{R}^{d}$ of dimension $d=3$, composed of a fiber phase $\Omega_{f}\subset\Omega$, and an elastic matrix $\hat{\Omega}\coloneqq\Omega\setminus\Omega_{f}\subset\mathbb{R}^{d}$, which we assume to be a connected, Lipschitz domain. We describe each of the $n_{f}\in\mathbb{N}$ fibers with a connected, Lipschitz domain and $\Omega_{f}$ is obtained as the union of these (possibly overlapping) domains.

Remark**.**

The results of this section hold for a general domain $\Omega_{f}$, union of multiple components with the required regularity. The property of the fibers of being thin, elongated, structures, is only needed for the model of Section 3, and plays no role in this section.

Given a continuous displacement field $u\colon\Omega\rightarrow\mathbb{R}^{d}$, representing a deformation from the equilibrium configuration, the corresponding stress tensor on $\Omega$ can be expressed using the stress-strain law Gurtin1981 :

[TABLE]

where $\mathbb{C}$ is a symmetric $4^{th}$ order tensor that takes the form:

[TABLE]

Here $\mathbb{C}_{\Omega}$ and $\mathbb{C}_{f}$ are assumed to be constant over their respective domains, and represent the elasticity tensors of the elastic matrix and of the fibers. The assumption of perfect bonds between the fibers and the elastic matrix (see (callister2007materials, , Chapter 16)) allows one to define the full order problem as a single elasticity problem on the union of the two domains, with pecewise elastic properties. This assumption may be relaxed by reformulating the fiber elasticity equations and the elastic matrix elasticity equations separately, and coupling them with appropriate interface (or bonding) conditions.

The classical formulation of static linear elasticity can be thought of as a force balance equation (see, for example, Gurtin1981 ):

Problem 1** (Classic Strong Formulation).**

Given an external force density field $b$, find the displacement $u$ such that

[TABLE]

Due to the piecewise nature of $\mathbb{C}$, it is natural to reformulate Problem 1 into a variational or weak formulation. Given a subset $D$ of $\partial\Omega$, we introduce the notation for the subspace of the Sobolev space $H^{1}(\Omega)$ with functions vanishing on a subset $D$ of the boundary $\partial\Omega$:

[TABLE]

with norm $\|\cdot\|_{V}=\|\cdot\|_{H^{1}(\Omega)}\coloneqq\|\cdot\|_{\Omega}+\|\nabla\cdot\|_{\Omega}$, where the symbol $\|\cdot\|_{A}$ represents the $L^{2}(A)$ norm over the measurable set $A\subset\Omega$, and $(\cdot,\cdot)_{A}$ represents the $L^{2}$ scalar product on the given domain $A$. We define the following space:

[TABLE]

The standard weak formulation reads:

Problem 2** (Classic Weak Formulation).**

Given $b\in L^{2}(\Omega)^{d}$, find $u\in V$ such that:

[TABLE]

The main idea behind our reformulation is to rewrite Problem 2 into an equivalent form, where we define two independent functional spaces. The novelty we introduce is to define the functional spaces on $\Omega$ and $\Omega_{f}$, and not on $\hat{\Omega}$ and $\Omega_{f}$. To achieved this, we define two fictitious materials: one with the same properties of the elastic matrix, occupying the full space $\Omega$, and one describing the “excess elasticity” of the fibers separately, defined on $\Omega_{f}$ only. The first step in this direction is to split the left-hand side of Equation (3) on the two domains:

[TABLE]

where \delta\mathbb{C}_{f}\coloneqq\mathbb{C}_{f}-{\left.\kern-1.2pt(\mathbb{C}_{\Omega})\vphantom{\big{|}}\right|_{\Omega_{f}}}.

For simplicity, we improperly use the expression “elastic matrix equation” and “fiber equation”, even though they should be really considered as the “whole domain equation”, and a “delta fiber equation”. This formal separation does not change the original variational problem, which can still be stated explicitly: given $b\in L^{2}(\Omega)^{d}$, find $u\in V$ such that:

[TABLE]

The boundary of the domain $\Omega$ induces a natural splitting on the boundary of the fibers: we define the following partition of $\partial\Omega_{f}$:

[TABLE]

where $B_{i}$ is the interface between the fibers and the elastic matrix, while $B_{e}$ is the interface between the fibers the exterior part of $\Omega$, that lies on the boundary $\partial\Omega$. Next we define the restriction of $H_{0}^{1}(\Omega)^{d}$ on the fibers:

[TABLE]

and separate the solution of Problem 2 in two components, one describing the whole material, the other describing only the fibers.

When reformulating the problem with two independent variables, the perfect bond condition between the fibers and the elastic matrix must be imposed via a volumetric non-slip constraint on the solution $(u,w)\in V\times W$:

[TABLE]

The described setting is similar to the distributed Lagrange multiplier method used to model fluid structure interaction problems with non-matching discretisations, as described in boffi2003finite ; auricchio2015fictitious ; boffi2014mixed .

The result is a constrained minimization problem:

[TABLE]

where we defined the total elastic energy of the system as

[TABLE]

To impose the non-slip constraint of Equation (8), we use the duality product $W^{\prime}\times W$ as in boffi2017fictitious . We define $Q\coloneqq W^{\prime}$, and enforce the perfect bond on the fibers by asking that

[TABLE]

For simplicity we will omit the subscript $Q\times W$ from now on.

The constrained minimization expressed in Equation (9) is equivalent to the saddle point problem:

[TABLE]

where the constraint is imposed weakly as in 11 through a Lagrange multiplier:

[TABLE]

A solution to Equation (12) is obtained by solving the Euler-Lagrange equation:

[TABLE]

that is:

Problem 3** (Saddle Point Weak Formulation).**

Given $b\in L^{2}(\Omega)^{d}$, find $u\in V,w\in W,\lambda\in Q$ such that:

[TABLE]

or, equivalently,

[TABLE]

where

[TABLE]

2.2 Well-posedness, existence and uniqueness

Existence and uniqueness follow from standard saddle point theory boffi2013mixed . We introduce the product Hilbert space, with its norm:

[TABLE]

and we indicate with $\boldsymbol{u}\coloneqq(u,w),\boldsymbol{v}\coloneqq(v,y)$ the elements of $\mathbb{V}$. We define the bilinear forms

[TABLE]

with which we reformulate the problem as: find $\boldsymbol{u}\in\mathbb{V},\lambda\in Q$ such that

[TABLE]

Proposition 2.1**.**

There exists a constant $\alpha_{1}>0$ such that:

[TABLE]

Moreover $\alpha_{1}=1$.

The proof for this proposition, and its discrete version, are variations on the one found in boffi2017fictitious .

Proof.

The non slip condition is given by the duality pairing between $Q=W^{\prime}$ and $W$; by definition of the norm in the dual space $Q$:

[TABLE]

where the last inequality can be proven fixing $v=0$. The final statement is found dividing by $\|q\|_{Q}$ and taking the $\inf_{q\in Q}$. ∎

Notice that $\alpha_{1}=1$, and does not depend on the two domains or on the two spaces.

Proposition 2.2**.**

Assume $\mathbb{C}_{\Omega}$ and $\mathbb{C}_{f}$ to be strongly elliptic with constants $c_{\Omega}$ and $c_{f}$ respectively, such that $c_{f}>c_{\Omega}>0$; then there exists a constant $\alpha_{2}>0$ such that:

[TABLE]

Proof.

An immediate consequence of the hypotheses is that $\delta\mathbb{C}_{f}$ is elliptic of constant $c_{f}-c_{\Omega}$. Following the proof for a similar statement found in auricchio2015fictitious ; boffi2014mixed , given an element $(v,y)\in\ker(\mathbb{E})$, the fact that {\left.\kern-1.2ptv\vphantom{\big{|}}\right|_{\Omega_{a}}}=y allows to use the Poincaré inequality on $v$ to control the norm of $y$; for every $(u,w)\in\ker(\mathbb{E})$:

[TABLE]

where we used the Poincaré inequality with its positive constant $c_{p}$ on $u\in V$, with $\alpha_{2}~{}\coloneqq~{}\frac{c_{p}\min(c_{\Omega},c_{f}-c_{\Omega})}{2}$, and we used the scalar product of $H^{1}(\Omega)^{d}$:

[TABLE]

and the analogous one for $H^{1}(\Omega_{f})^{d}$. The final statement is obtained dividing by $\|(u,w)\|_{\mathbb{V}}$ and considering the $\inf_{(u,w)\in\ker{\mathbb{E}}}$. ∎

Remark**.**

This paper does not intend to focus on the choice of elastic tensors. Strong ellipticity is a common property among them, and holds in the case of linearly elastic materials (see e.g. marsden1994mathematical ): let $u\in V$, $w\in W$

[TABLE]

where $\mu_{\delta}\coloneqq\mu_{f}-\mu_{\Omega}$ and $\lambda_{\delta}\coloneqq\lambda_{f}-\lambda_{\Omega}$, and $Eu=\frac{\nabla u+\nabla u^{T}}{2}$ is the symmetric gradient. These are the elastic tensors we use in Section 4, for our numerical tests.

Propositions 2.2 and 2.1 imply that the inf-sup conditions are satisfied, and that Problem 3 is well-posed, and has a unique solution.

2.3 Finite element discretization

The distributed Lagrange formulation of Problem 3 makes it possible to use independent triangulations for its numerical solution. Consider the family $\mathcal{T}_{h}(\Omega)$ of regular meshes in $\Omega$, and a family $\mathcal{T}_{h}(\Omega_{f})$ of regular meshes in $\Omega_{f}$, where we denote by $h$ the maximum diameter of the elements of the two triangulations. We assume that no geometrical error is committed when meshing, i.e., $\Omega=\bigcup_{T_{h}\in\mathcal{T}_{h}(\Omega)}T_{h}$, and $\Omega_{f}=\bigcup_{S_{h}\in\mathcal{T}_{h}(\Omega_{f})}S_{h}$. We consider two independent finite element spaces $V_{h}\subset V$, and $W_{h}\subset W$, and we take $Q_{h}=W_{h}$. All finite dimensional spaces are endowed with the norms of their continuous counterparts.

Problem 4** (Discrete Weak Formulation).**

Given $b\in L^{2}(\Omega)^{d}$, find $u_{h}\in V_{h},~{}w_{h}\in W_{h},~{}\lambda_{h}\in Q_{h}$ such that:

[TABLE]

Existence, well posedness, and convergence are guaranteed if the classical discrete inf-sup conditions are satisfied:

Proposition 2.3**.**

Assume that the $L^{2}$ projection

[TABLE]

is continuous and $H^{1}$-stable, i.e., there exists a positive constant $c_{W}$ such that for all $w\in W$:

[TABLE]

Then there exists a constant $\alpha_{3}>0$, independent of $h$, such that:

[TABLE]

Proof.

For every $q_{h}\in Q_{h}$, by definition of the $Q$ norm, and by property (26), there exists $\hat{w}\in W$ such that:

[TABLE]

Therefore,

[TABLE]

and we conclude as in Proposition 2.1. ∎

Proposition 2.4**.**

Assume that $\mathbb{C}_{\Omega}$ and $\mathbb{C}_{f}$ are strongly elliptic with constants $c_{\Omega}$ and $c_{f}$ respectively, such that $c_{f}>c_{\Omega}>0$; then there exists a constant $\alpha_{4}>0$, independent of $h$, such that:

[TABLE]

where

[TABLE]

The proof follows the one of Proposition 2.2.

Error estimate

Proposition 2.1 and 2.2 allow us to apply the theory from Chapter $5$ of boffi2013mixed , obtaining the following error estimate:

Theorem 2.1**.**

Consider $\mathbb{C}_{\Omega}$ and $\mathbb{C}_{f}$, elastic stress tensors satisfying the hypothesis of Proposition 2.2, the domains $\Omega$ and $\Omega_{f}$ with the regularity required in Section 2, and $b\in L^{2}(\Omega)^{d}$. Then the following error estimate holds for $(u,w,\lambda)$, solution to Problem 3, and $(u_{h},w_{h},\lambda_{h})$, solution of Problem 4:

[TABLE]

where $C_{e}>0$, and depends on $\alpha_{3},\alpha_{4},c_{\Omega},c_{f}$ and the norm of the operators $\mathcal{K}_{\Omega}$ and $\mathcal{K}_{f}$.

We remark how the constant $C_{e}$ is affected by the coupling between the two meshes; as intuition suggests, the quality of the solution does not depend only the on the ability of $V$ and $W$ to individually describe it, but also on the coupling between them.

Non-matching meshes

One of the basic assumptions made in the continuous case is that, for every element $v\in V$, we have that {\left.\kern-1.2ptv\vphantom{\big{|}}\right|_{\Omega_{f}}}\in W; similarly every element $w\in W$ can be extended to an element of $V$.

With an independent discretization of the two meshes, the inclusion $W_{h}\subset V_{h}$ is in general false. In this case, the two requirements for a good approximation of the problem are that the projection $P_{W}\colon W\rightarrow W_{h}$ is $H^{1}$-stable, and that the kernel of $\mathbb{E}_{h}$ is “rich enough”.

The $H^{1}$ stability condition can be easily obtained by inverse inequalities on quasi-uniform meshes. For more general meshes, the stability of the $L^{2}$ projection has been investigated, for example, in Crouzeix2001 ; bramble2002stability ; Bank2013 . The non-matching nature of the discretization does not deteriorate the approximation properties of the underlying spaces, provided that the two discretizations have comparable local mesh sizes.

In particular, the triangle inequality and Bramble-Hilbert lemma imply readily that for any $u$ in $V$ we can write

[TABLE]

i.e., even with non-matching meshes, the passage through the two non-matching discretizations still allows one to control in a straight forward way the $L^{2}$ norm of the error for any $H^{1}$ function on $\Omega$. Provided that the discretizations on $\Omega$ and $\Omega_{f}$ have comparable mesh sizes (respectively $h_{V}$ and $h_{W}$ in the equation above), then the passage through non-matching meshes keeps the error in the same order.

Moreover, if we integrate exactly on the non-matching grids, it is possible to guarantee that globally constant and linear functions are included in the kernel, ensuring that $\ker(\mathbb{E}_{h})~{}\neq~{}\{(0,0)\}$.

3 Thin fibers

The computational cost of discretizing numerous three-dimensional fibers might render Problem 4 too computationally demanding: a possible simplification is to approximate the fibers with one-dimensional structures. This construction is a non-trivial because the restriction (or trace) of a three-dimensional function to a one-dimensional domain is not well defined in the $H^{1}$ Sobolev space, which is the natural space where the elasticity problem is well posed.

Instead of resorting to weighted Sobolev spaces and graded meshes, as done in d2008coupling ; d2012finite , we define the coupling between the one-dimensional fibers and the three-dimensional elastic matrix through an averaging technique that renders the problem well posed.

To simplify the exposition, we consider a single fiber with constant radius $a$; the same results hold for a finite collection of fibers. Given a one-dimensional curve $\Gamma$ immersed in $\mathbb{R}^{3}$, we call fiber its tubular extension of radius $a>0$:

[TABLE]

such that $\partial\Omega_{f}$ is non-intersecting.

Given a point $x$ on $\Gamma$, we call $D_{a}(x)$ the disk of radius $a$ orthogonal to $\Gamma$, and we assume that the fibers can be written as the image of a diffeomorphism

[TABLE]

where $\Omega_{a}$ is a straight cylinder of radius $a$ and of axis $\Upsilon$ aligned with the coordinate $x_{1}$, and $\Phi$ satisfies the following hypotheses:

1. [i)]

$\Phi(\Omega_{a})=\Omega_{f}$ 3. 2.

$\Phi\left(\Upsilon\right)=\Gamma$ 4. 3.

$D_{a}(\Phi(x))=\Phi\left(D_{a}(x)\right)$, for every $x\in\Upsilon$.

Given two vectors $u,v\in\mathbb{R}^{n}$, their tensor product is a matrix of size $n\times n$, denoted as $u\otimes v$; defined for each index $i,j\in\mathbb{N}$, $1\leq i,j\leq n$ as:

[TABLE]

Let $W\coloneqq H^{1}(\Omega_{f})^{d}$, and $W_{\Gamma}\coloneqq H^{1}(\Gamma)^{d}$. For $w\in W_{\Gamma}$, the surface gradient along the curve is defined as:

[TABLE]

where $t$ is the unitary tangent vector to the curve $\Gamma$ and $\varphi_{w}$ is a smooth extension of $w$ in a tubular neighborhood of $\Gamma$.

The following result can be used to define the coupling between the three dimensional elastic matrix and the one dimensional fibers:

Theorem 3.1**.**

Let $\Omega_{f}$ be a fiber of radius $a$ with center line $\Gamma$.

The operator $\mathcal{R}\colon W\rightarrow W_{\Gamma}$, defined as

[TABLE]

is bounded and continuous.

Proof.

We start by considering straight cylinders, and smooth functions. Let $u\in C^{\infty}(\Omega_{a})$; the following inequalities hold:

[TABLE]

where $\Omega_{a}$ is a straight cylinder of axis $\Upsilon$.

Consider the coordinates $(x_{1},x_{2},x_{3})$, where $x_{1}$ is aligned with the cylinder’s axis and $x_{2},x_{3}$ are normal to the axis. The first inequality can be proven directly, using either Jensen’s or Hölder inequality:

[TABLE]

The second inequality can be proven in a similar fashion, by splitting the integral used to compute $\mathcal{R}u(x)$ on $\Gamma$ and on the domain $D\coloneqq\{x_{2},x_{3}\text{ s.t. }x_{2}^{2}+x_{3}^{2}\leq a^{2}\}$, which does not depend on $x_{1}$, and observing that $|\nabla_{\Upsilon}\mathcal{R}u(x)|^{2}=|\partial_{x_{1}}\mathcal{R}u(x)|^{2}$.

From (33) and (34), we conclude with a density argument that

[TABLE]

The generalization to a generic fiber $\Omega_{f}=\Phi(\Omega_{a})$ follows applying a change of coordinate transformation through $\Phi$. In particular, for any $u\in H^{1}(\Omega_{a})$:

[TABLE]

where we defined the following positive constants:

[TABLE]

For the first inequality, we begin with the $L^{2}$ part of the norm:

[TABLE]

Similarly, for $\|\nabla(u\circ\Phi)\|_{\Omega_{a}}^{2}$:

[TABLE]

To prove the second inequality we follow a similar procedure, splitting the integral over the centerline of the tubular neighbourhood and the restrictions of $\Phi$ to the perpendicular disks. The combined inequalities imply the thesis:

[TABLE]

∎

A possible right inverse of the restriction operator is the extension operator $\mathcal{E}$:

[TABLE]

where $P_{\Gamma}$ is the geometric projection from the domain $\Omega_{f}$ to $\Gamma$, i.e., for every $y\in\Omega_{f}$, $P_{\Gamma}y$ is the element of $\Gamma$ such that:

[TABLE]

Equation (30) guarantees that $P_{\Gamma}$ is well-defined. The operator $\mathcal{E}$ is clearly bounded in $L^{2}$; to prove that its gradient is bounded it is sufficient to adapt the proof of Theorem 3.2.

The following result allows us to define a function $g$ on $\Gamma$, describing the geometry of $\Omega_{f}$, and reduce our computations to a linear integration.

Theorem 3.2**.**

For a general tubular neighbourhood $\Omega_{f}$ of a curve $\Gamma$, for which the curvature $\kappa$ and the torsion $\tau$ are defined a.e., there exists a function $g$ defined on $\Gamma$ satisfying:

[TABLE]

for every $w\in H^{1}(\Gamma)^{d}$, with $c_{\Gamma}\coloneqq\pi a^{2}$.

Proof.

Let $I$ be an interval, and let $\omega\colon I\rightarrow\Gamma$ be the arclength parametrization of $\Gamma$. From the hypotheses, it is possible to define an orthonormal basis on $\Gamma$: the tangent $t(s)$, the normal $n(s)$, and the binormal $b(s)$, along with the curvature $\kappa(s)$ and the torsion $\tau(s)$ (see, e.g., (kreyszig, )).

We consider a change of coordinates w.r.t. $t$, $n$, and $b$, which is typically used in physics to study wave propagation, optics and particle trajectories lee2018accelerator ; tang1970orthogonal ; ryder1993nonlinear , based on the coordinates $(r,\vartheta,s)$, with the orthogonal metric:

[TABLE]

and the explicit formula for the gradient:

[TABLE]

For a detailed description of these formulas, and the exact definition of the angles $\vartheta$ and $\hat{\theta}$ see (ryder1993nonlinear, , Chapter 3).

Given $w\in H^{1}(\Gamma)^{d}$ the extension $\mathcal{E}w$ is constant on each orthogonal disk, i.e., when computing the gradient with Equation (39) the components $e_{r}$ and $e_{\theta}$ are [math].

Thus:

[TABLE]

Where we defined the function:

[TABLE]

∎

The function $g$ encodes a geometrical stiffness information, which reflects the fact that rods with finite thickness store more energy (i.e., the function $g$ grows) when they are bent with high curvature w.r.t. to their thickness (i.e., when $\kappa r$ is close to one). In the case of a straight cylinder $\kappa(s)=0$ a.e., and $g$ reduces to the constant one.

The condition $r\leq a<\max_{s}1/\kappa(s)$, which is necessary to guarantee that $\Omega_{f}$ is non self-intersecting (duistermaat2004multidimensional, ), plays a major role also in defining the geometrical stiffness $g$ on the whole disk. Under these hypotheses it is also possible to prove that there exists a constant $C_{g}>0$ such that $g(s)\geq C_{g}>0$.

3.1 Problem Formulation

When the fibers are thin compared to the domain size (i.e., $a\ll\text{diam}(\Omega))$, it is reasonable to make the following assumptions:

1. [i)]

during the deformation, the radius of the fibers does not change; 3. 2.

the deformation field inside the fiber is similar to the deformation on the centerline of the fiber itself.

If we consider then a solution $u$ to Problem (3), we expect that

[TABLE]

jusifying the replacement of the space $W\subseteq H^{1}(\Omega_{f})^{d}$ with the much smaller space $\mathcal{E}W_{\Gamma}\subset W\subseteq H^{1}(\Omega_{f})^{d}$.

If we restrict functions in $W$ in Problem (4) to the subspace $\mathcal{E}W_{\Gamma}$, and exploit Theorem 3.2 and the fact that $\mathcal{R}\mathcal{E}w=w$ for any $w\in W_{\Gamma}$, we obtain the following energy functional for thin fibers

[TABLE]

where the non-slip condition on $\Omega_{f}$ has been replaced by an average non-slip condition:

[TABLE]

where $Q_{\Gamma}\coloneqq W_{\Gamma}^{\prime}$, and here and below the notation $\langle\cdot,\cdot\rangle$ is used to indicate the duality product between $Q_{\Gamma}^{\prime}=W_{\Gamma}$ and $Q_{\Gamma}:=W_{\Gamma}^{\prime}$.

We obtain the following formulation for the coupling of thin fibers with a thick elastic matrix:

Problem 5** (1D-3D Weak Formulation).**

Given $b\in L^{2}(\Omega)^{d}$, find $u\in V,w\in W_{\Gamma},\lambda\in Q_{\Gamma}$ such that:

[TABLE]

Similarly to what was done in the full order problem, we define the space $\mathbb{V}\coloneqq V\times W_{\Gamma}$, its norm $\|(\cdot,\cdot)\|_{\mathbb{V}}\coloneqq\|\cdot\|_{V}+\|\cdot\|_{W_{\Gamma}}$, and the operators:

[TABLE]

3.2 Well-posedness, existence and uniqueness

Proposition 3.1**.**

There exists a constant $\alpha_{5}>0$ such that:

[TABLE]

where:

[TABLE]

Moreover $\alpha_{5}=1$.

Proof.

The non slip condition is given by the duality pairing between $Q_{\Gamma}=W_{\Gamma}^{\prime}$ and $W$; by definition of the norm in the dual space $Q_{\Gamma}$:

[TABLE]

where the last inequality can be proven fixing $v=0$.

The final statement is found dividing by $\|q\|_{Q_{\Gamma}}$ and taking the $\inf_{q\in Q}$. ∎

Proposition 3.2**.**

Assume $\mathbb{C}_{\Omega}$ and $\mathbb{C}_{f}$ are strongly elliptic, with constants $c_{\Omega}$ and $c_{f}$ respectively such that $c_{\Omega}>c_{f}>0$; there exists a constant $\alpha_{6}>0$ such that:

[TABLE]

Proof.

For every $(u,w)\in\ker\mathbb{E}_{T}$:

[TABLE]

with $\alpha_{6}\coloneqq\frac{c_{p}c_{\Gamma}c_{b}\min(c_{\Omega},C_{g}(c_{f}-c_{\Omega}))}{2}$.

The result is obtained dividing and considering the $\inf_{(u,w)\in\ker{\mathbb{E}_{T}}}$. ∎

Using the saddle point theory we conclude that Problem 5 is well-posed, and has a unique solution.

3.3 Finite element discretization

The discretization of Problem 5 for thin fibers follows the steps of Section 2.3, on the domains $\Omega$ and $\Gamma$ (the fiber’s one-dimensional core). Consider two independent discretizations for these domains; the family $\mathcal{T}_{h}(\Omega)$ of regular meshes in $\Omega$, and a family $\mathcal{T}_{h}(\Gamma)$ of regular meshes in $\Gamma$. We assume no geometrical error is committed when meshing and consider two independent finite element discretizations $V_{h}\subset V,\ W_{h}\subset W_{\Gamma}$, and let $Q_{h}=W_{h}$.

Problem 6** (1D-3D Discretized Weak Formulation).**

Given $b\in L^{2}(\Omega)^{d}$, find $u_{h}~{}\in~{}V_{h},w_{h}\in W_{h},\lambda_{h}\in Q_{h}$ such that:

[TABLE]

Again, existence, well posedness and convergence are guaranteed if the classical inf-sup conditions are satisfied:

Proposition 3.3**.**

Assume the $L^{2}$ projection $P_{W_{\Gamma}}\colon W_{\Gamma}\rightarrow W_{h}$ is continuous and $H^{1}$-stable, as in Proposition 2.3; then there exists a constant $\alpha_{7}>0$, independent of $h$, such that:

[TABLE]

Proposition 3.4**.**

Assume $\mathbb{C}_{\Omega}$ and $\mathbb{C}_{f}$ are strongly elliptic, with constants $c_{\Omega}$ and $c_{f}$ respectively such that $c_{f}>c_{\Omega}>0$; then there exists a constant $\alpha_{8}>0$, independent of $h$, such that:

[TABLE]

where:

[TABLE]

The proofs of these propositions are variations on the proof Propositions 2.3 and 2.4. It is possible to prove an estimate analogous to the one of Theorem2.1:

Theorem 3.3**.**

Consider $\mathbb{C}_{\Omega}$ and $\mathbb{C}_{f}$, elastic stress tensors satisfying the hypothesis of Proposition 3.2, the domains $\Omega$ and $\Omega_{f}$ with the regularity required in this section, and $b\in L^{2}(\Omega)^{d}$. Then the following error estimate holds for $(u,w,\lambda)$, solution to Problem5, and $(u_{h},w_{h},\lambda_{h})$, solution of Problem 6:

[TABLE]

where $C_{e}>0$, and depends on $\alpha_{5},\alpha_{6},c_{\Omega},c_{f}$ and the norms of the operators $\mathcal{K}_{\Omega}\colon V\rightarrow V^{\prime}$ defined as $\langle\mathcal{K}_{\Omega}u,\cdot\rangle\coloneqq(\mathbb{C}_{\Omega}\nabla u,\nabla\cdot)_{\Omega}$, and $\mathcal{K}_{f}\colon W_{\Gamma}\rightarrow Q_{\Gamma}$ defined as $\langle\mathcal{K}_{f}w,\cdot\rangle\coloneqq(\delta\mathbb{C}_{f}\nabla w,\nabla\cdot)_{\Omega_{f}}$.

4 Numerical validation

The analytical solution of Problems 3 and 5, even for simple configurations, is non-trivial: we chose some FRMs structures which are studied in literature, and used the known approximated solutions as a comparison for our model.

Using the deal.II library dealii-2020 ; dealII90 ; BangerthHartmannKanschat2007 ; MaierBardelloniHeltai-2016-a ; SartoriGiulianiBardelloni-2018-a , and the deal.II step-60 tutorial step60 we developed a model for thin fibers proposed in Section 3, and compared it with the Rule of Mixtures and the Halpin-Tsai configurations in some pull and push tests.

Numerical Setting

To solve Problem 6 on a collection of fibers we define:

1. [i)]

$V_{h}=\text{span}\{v_{i}\}_{i=1}^{N}\subset H^{1}(\Omega)$, finite element space of dimension $N\in\mathbb{N}$ defined on $\Omega$, the elastic matrix; 3. 2.

$W_{h}=\text{span}\{w_{a}\}_{i=a}^{M}\subset H^{1}(\Gamma)$ the finite element discretization of dimension $M\in\mathbb{N}$ defined on the collection of $n_{f}\in\mathbb{N}$ fibers $\Gamma\coloneqq\bigcup_{k=1}^{n_{f}}\Gamma_{k}\subset\Omega$; 4. 3.

$Q_{h}$: the space of the Lagrange multiplier, discretized as $W_{h}$.

We use $i,j$ as indices for the space $V_{h}$ and $a,b$ as indices for the space $W_{h}$, and assume all hypotheses on spaces and meshes of Section 3.3 are satisfied. For each fiber $\Gamma_{k}$ we define its tubular neighborhood $\Omega_{a,k}$, and define

[TABLE]

In our implementation, we compute integrals over $\Omega_{f}$ by quadrature formulas, computing the integration over the orthogonal disks $D_{a}(x)$ to $\Gamma$ using the mid point rule.

The restriction of the operators of Problem 6 to the discrete finite element spaces produce the following sparse matrices:

[TABLE]

Here $B$ is the coupling matrix from $V_{h}$ to $W_{h}$, $M$ the mass-matrix of $W_{h}$. The discretization of Problem 6 becomes: find $(u,w,\lambda_{h})\in V_{h}\times W_{h}\times Q_{h}$ such that

[TABLE]

where $g_{i}\coloneqq(b,v_{i})_{\Omega}$. We reduce the system through the solution of the following linear problems:

[TABLE]

obtaining:

[TABLE]

where $P_{\Gamma}\coloneqq L^{-1}B$. Boundary conditions are imposed weakly, using Nitsche method as done in RotundoKimJiang-2016-a .

4.1 Model description

For our simulations we use linear finite elements and uniformly refined hexaedral meshes. The elastic matrix is the unitary cube $\Omega\coloneqq[0,1]^{3}$, where for boundary conditions we refer to Figure 3b. The used elastic tensors are described in Equations 20-22; for the model description we use the following parameters:

• •

$r_{\Omega},r_{\Gamma}$: global refinements of the meshes describing $\Omega$ and $\Gamma$ respectively.

• •

$\lambda_{\Omega},\mu_{\Omega},\lambda_{f},\mu_{f}$: Lamé parameters for the elastic matrix and the fibers respectively.

• •

$\beta$: fiber volume ratio or representative volume element (RVE), i.e., $\beta=|\Omega_{f}|/|\Omega|$.

• •

$a$: the radius of the fibers.

4.2 Homogeneous fibers

We consider a unidirectional composite, where fibers are uniform in properties and diameter, continuous, and parallel throughout the composite $\Omega$ (see Figure 3a). We compare the results obtained with our model with the ones obtained using the Rule of Mixtures (hashin1965elastic, ; agarwal2017analysis, ) to approximate the stress-strain equation:

[TABLE]

This approximation agrees with experimental tests especially for tensile loads, and when the fiber ratio $\beta$ is small.

Multiple tests were run, keeping $\beta$ constant, while increasing the fiber density and reducing the fiber diameter; we expect this process to render the coupled model solution increasingly close to the homogenized one.

Comparing solutions

Figures 4a and 4b illustrate the influence of $\Omega$’s refinement on the final result, when using few fibers on a pull test. Stiff fibers oppose being stretched, deforming the elastic matrix $\Omega$ through the non-slip condition: near each fiber, the deformation of $\Omega$ should be symmetrical, resembling a cone. This effect is better described in Figure 4b, where the higher value of $r_{\Omega}$ results in greater geometrical flexibility of the elastic matrix, allowing a better description of the effect of each fiber. Lower values of $r_{\Omega}$ result in a non symmetrical solution, as in Figure 4a. The lower geometrical flexibility results in an “averaged” solution which, in the case of few fibers, is closer to the homogenized model.

The slopes in the plots have been computed by taking the ratio of the error in two subsequent simulation, i.e.,

[TABLE]

where $\mathbf{e}_{1},\mathbf{e}_{2}$ be the $L^{2}$ errors and $n_{f,1},n_{f,2}$ be the number of fibers for two subsequent simulations.

Pull test along fibers

Dirichlet homogeneous conditions is applied to face [math], Neumann homogeneous conditions is applied to faces $2,3,4,5$. In the Push Test, the Neumann condition $0.05$ is applied to face $1$. For the Pull Test, the value $-0.05$ is applied to the same face. Boundary conditions are applied only to $\partial\Omega$; the fibers interact through the coupling with the elastic matrix. We report here only data from pull tests, as push tests gave comparable results.

The use of the projection matrix $P_{\Gamma}\colon V_{h}\rightarrow W_{h}$, and the error estimate for the fully three-dimensional case (Inequality 28), both suggest that the solution quality on the elastic matrix depends on both $V_{h}$ and $W_{h}$. This is apparent in Figure 5a, where for $r_{\Gamma}=1$ the mesh of $\Gamma$ is unable to describe the stretch of the material, resulting in the error remaining approximately constant after a certain fiber density is reached. A similar behaviour emerges in the case $r_{\Gamma}=2$. Figure 5b shows that refining only the elastic matrix does not improve the solution quality: as the number of fibers increases, the error converges to approximately the same value, which is limited by $r_{\Gamma}$.

Figure 6a shows an error comparison as the value of $\mu_{F}$ varies: as expected our model is better suited for stiff fibers.

4.3 Random fibers

We consider a pull test on a random chopped fiber reinforced composite: we distribute small fibers at a random points of $\Omega$, with a random direction parallel to the $\langle x,y\rangle$ plane; the fibers share the same size and properties. If a fiber surpasses the edge of $\Omega$, it is cut.

For more details on the algorithm used to distribute the fibers see the Random Sequential Adsorption algorithm pan2008analysis ; our implementation generates only the plane angle, and does not implement an intersection-avoidance mechanism.

As a comparison model, we estimate the material parameters using the empirical Halpin-Tsai halpin1976halpin ; which we report here for longitudinal moduli, as described in agarwal2017analysis . The fibers have length $l$, diameter $d$, the fiber and the matrix Young moduli are $E_{f}$ and $E_{m}$ respectively, $\beta$ is the volume fraction occupied by the fibers. We define two empirical constants:

[TABLE]

This allows to compute the longitudinal and transverse moduli for aligned short fibers:

[TABLE]

If fibers are randomly oriented in a plane the following equations can be used to predict the elastic modulus:

[TABLE]

Since a random fiber composite is considered isotropic in its plane, the Poisson’s ratio can be calculated as:

[TABLE]

According to this approximation, the properties of this composite do not depend directly on the fiber length or radius, but on the aspect ratio $l/d$. Our test setting runs on the unitary cube, with a fiber ratio $\beta=0.135$ and a fiber aspect ratio of $\frac{l}{2r}\approx 10$, where $l$ is the fiber’s length and $r$ is the fiber’s radius; the values used are described in Table 1.

We could not find an exact estimate of the error convergence, but we expect the solution to improve as the number of fibers increases because:

• •

the fiber radius $a$ reduces, improving of the average non-slip condition,

• •

more fibers result in a more homogeneous material on the planes parallel to the $\langle x,y\rangle$ plane.

Following pan2008analysis , we consider a short fiber E-glass/urethane composite: the fiber and matrix Young’s modulus are, respectively, $E_{f}~{}=~{}70GPa$ and $E_{m}~{}=~{}3GPa$, while the Poisson ratios are $\nu_{f}~{}=~{}0.2$ and $\nu_{m}~{}=~{}0.38$. These values werex converted to the Lamé parameters using the classic formulas for hyper elastic materials.

The predicted parameters for the composite are: $E_{C}~{}=~{}2.20GPa$ and $\nu_{C}~{}=~{}0.38GPa$; these are slightly different from pan2008analysis because, in the Halpin-Tsai equations, $l/d$ was used instead of $2l/d$ . The boundary conditions used for the pull tests are the same of Paragraph 4.2.

We limit the global refinements of $\Gamma$, in order to obtain cells of approximately the same size on both $\Omega$ and the fibers. The results are shown in Figure 6b: as the number of fiber increases the error reduces, but because the random fiber model is more complex than the homogeneous one, the final error achieved is higher than the one reached in the previous test.

5 Conclusions

Starting from a linearly elastic description of bi-phasic materials, we derive a new formulation for fiber reinforced materials where independent meshes are used to discretize the fibers and the elastic matrix, and the copling between the two phases is obtained via a distributed Lagrange multiplier.

We prove existence and uniqueness of a solution for the final saddle point problem, and we analyze a simplified model where fibers are discretised as one-dimensional.

The model is validated against the Rule of Mixtures and the Hapin-Tsai equations, where we test our discretisation with uniform and random distributions of fibers. The true benefit of our model, however, lies in the possibility to tackle complex and intricated fiber structures independently from the background elastic matrix discretization, opening the way to the efficient simulation of complex multi-phase materials.

From the numerical analysis point of view, there are some issues that deserve further development, such as finding better preconditioners for the final system, and exploring different solutions for the one dimensional coupling operator.

The formulation of our method makes it particularly suited for extensions to more complex situations, e.g., three-phasic materials, or materials where perfect-bond is replaced by other, more realistic conditions.

Acknowledgements

The authors would like to thank Prof. Lucia Gastaldi, for the fruitful conversations and suggestions that greatly improved this manuscript, and to thank Dr. Giovanni Noselli, and Dr. Maicol Caponi for their priceless suggestions and insights in their fields of expertise.