Homogenization of composite materials reinforced with unidirectional fibres with complex curvilinear cross section: a virtual element approach

1.   Introduction

The use of composite materials for engineering applications is a continuously broadening field, due to their high performances under several points of view. Many are in fact the advantages offered by such materials, as for instance lighter weight, the ability to tailor layup for optimum strength and stiffness, improved fatigue life, corrosion resistance, and, through good design practice, reduced assembly costs due to fewer detail parts and fasteners. For these reasons these composites are used in several fields of engineering, for example in civil engineering are used for strengthening of reinforced concrete columns, in mechanical engineering in high performance racing cars, aerospace and biomechanical applications being as well a growing network in this regard. A special class in this realm is represented by fibre-reinforced composite materials. The specific strength of superior quality material fibres (especially carbon, glass, metal alloy, just to mention a few types) are higher than those of embedding matrices making these materials an interesting option for the construction of numerous technical devices. Being intrinsically heterogeneous, deriving homogenized equivalent material properties is mandatory in design practice. Amongst the many different methods developed in recent years, asymptotic homogenization has become widely known and quite popular due to the versatility in application with respect to material configuration of fibre and matrix arrangements. The method nonetheless requires a numerical approximation whenever the geometry of fibres or their disposition within the surrounding matrix becomes complex, being classical FE methods the usual tool in this context.

A powerful alternative to Finite Element Methods (FEM) and inherent limitations, is offered by the virtual element method (VEM), a recently introduced numerical method for approximation of PDEs ^[5,6] which can be viewed as an extension of FEM to general polygonal and polyhedral elements. The strongest aspects in favor of the VEM are its firm mathematical foundations, simplicity in implementation, and efficiency and accuracy in computations, as well as mesh adaptivity together with the possibility of having curvilinear polygonal meshes. The latter aspect in particular results quite efficient in the study of the overall mechanical behavior of fibre-reinforced composite materials through homogenization since the ensuing boundary value (cell) problem is posed, in general, on curvilinear domains bounded by material interfaces identified on the cross section orthogonal plane ^[3].

This work presents a curvilinear VEM approach for the computational homogenization of in-plane shear moduli for unidirectional fibre reinforced materials having inclusions with complex shape profiles. A key point of the procedure relies on an efficient 2D cubature algorithm for NURBS boundary polygons which grants positive weights and interior nodes. An extensive campaign of numerical applications illustrate accuracy and convergence patterns of the method, both for doubly periodic regular arrangements of fibres and for randomly distributed fibres within the matrix phase.

The paper is organized as follows: Section 2 presents the governing equations of the considered homogenization problem, Section 3 introduces the curvilinear VEM approximation space and the discretized form of the problem. Sections 4–7 are devoted to the cubature formula for NURBS-shaped curvilinear polygons. Last, Section 8 presents a set of selected numerical tests supporting the accuracy and efficiency of the proposed methodology. Conclusions are drawn in Section 9.

2.   Asymptotic homogenization of fibre-reinforced composite: the anti-plane problem

This section is devoted to a unified compact presentation of the so called computational asymptotic homogenization of antiplane shear moduli for both a doubly-period or a random fibre-reinforced composite. We consider a composite material, reinforced with long, parallel fibres, distributed in the material with a statistically homogeneous microstructure, given either by a doubly periodic or a random spatial arrangement. Fibres have all same cross section with a possibly complex shape: In the present context we assume that the curvilinear curve defining the boundary for a given cross section is described through a sigle NURBS or a $C^0$ regular blend of subsequent NURBS curves.

In either case, at microscale, the section orthogonal to fibres is represented by a doubly-periodic arrangement of repeating unit cells (RUC). A RUC is a parallelogram, having edges $L_1$, $L_2$, and an angle $\varphi$, containing a given number $F$ of fibres, as represented in Figure 1 for the two exemplary cases of doubly periodic and random type of composite. A given material setup is characterized by the so called volume fraction, denoted by $f = \sum_{j = 1}^{F}f_j$ which rapresents the ratio between fibre material area and RUC total area.

The case studied in this communication refers to effective in-plane elastic shear moduli, which are computed applying computational asymptotic homogenization through VEM, due to a lack of closed-form solutions of the aforementioned problem for fibres with complex cross section boundary. Statement of the problem is given by the following set of boundary value problems of equilibrium over the composite cross section domain:

indexed by a parameter $\varepsilon$. In the above, $w_{ \varepsilon}$ denotes the displacement field in the fibre axes direction, and $\cup_j\Omega ^{ {\rm f}}_{j \varepsilon}$ and $\Omega ^{ {\rm m}}_{ \varepsilon}$ indicate fibres and matrix domains, respectively, while$\cup_j\Gamma_{j \varepsilon}$ is defined as the set of all fibre/matrix interfaces, and $\boldsymbol{\nu}$ represents the unit normal to $\cup_j\Gamma_{j \varepsilon}$ towards $\Omega ^{ {\rm m}}_{ \varepsilon}$. The double square brackets operator ${[\kern-0.15em[}\cdot {]\kern-0.15em]}$ denotes jump across the fibre interface, defined as extra-minus-intra difference. Parameter $\varepsilon$ is a scaling factor for the microstructure, such that $\varepsilon = 1$ corresponds to the actual composite material, while homogenization limit is reached for the limit of $\varepsilon$ approaching zero.

The physical meaning of the above set of governing equations is as follows: Eq (2.1) represents translational equilibrium along fibre axis direction; Eq (2.2) represents in-plane equilibrium at fibre/matrix interface i.e., continuity of the normal-to-interface shear stress component; Eq (2.3) is actually a constitutive equation, also known as linear spring-layer model ^[19,21], ruling the zero-thickness imperfect interface behavior in terms of displacement jump ${[\kern-0.15em[} w_{ \varepsilon} {]\kern-0.15em]}$ and normal traction $\boldsymbol{G} \nabla w_{ \varepsilon} \cdot \boldsymbol{\nu}$ at the interface in such a way that parameter(s) $D_j$ are selected to represent the level of interface degradation ^[10]. Lastly, factor $\varepsilon^{-1}$ providess the right scaling for the homogenization limit ^[21].

Assuming homogeneous isotropic linear elastic fibre and matrix material behavior, the in-plane effective shear moduli, which are collected by the constitutive tensor $\boldsymbol{G}$, are given by:

where ${\bf I}$ is the unit tensor. An interesting generalization for fibres to incorporate cylindrical orthotropy, with material grading along the radial direction, could be also taken into account for the particular case of circular fibres (see ^[2] for more details), but it is here omitted for simplicity.

Well posedness of the homogenization problem requires the following conditions to hold:

where $G ^{ {\rm m}}$, $G^{\it f}$ are shear moduli, for the matrix and the fibre, respectively, $D_j$ is the fibre/matrix interface elastic stiffness parameter (see Eq (2.3)).

2.1.   Homogenization of in-plane shear moduli

To derive the homogenized or effective in-plane shear moduli of the composite material, asymptotic homogenization shall be employed herein. With reference to a given RUC, the procedure identifies two characteristic separate length scales for the spatial domain under consideration, i.e., two different space variables: a macroscopic one, $x$, and a microscopic one, $y = x/ \varepsilon$, $y \in {\mathsf{D}}$, being ${\mathsf{D}}$ the RUC (see Figure 1), whose extra-fibre space, intra-fibre space and fibre-matrix interface are denoted by ${\mathsf{D}} ^{ {\rm m}}$, ${\mathsf{D}}^{ {\rm f}}_j$ and $\Gamma_j$, for $j = 1\ldots F$, respectively. Thence, an asymptotic expansion of the primary variable, the axial displacement field, is considered with respect to powers of $\varepsilon$:

where $w_{0}$, $w_{1}$, $w_{2}$ are doubly periodic functions over the RUC domain, with $w_{1}$, $w_{2}$ having zero integral average over ${\mathsf{D}}$. Substituting this expansion in the equilibrium form Eqs (2.1)–(2.3) and equating the power-like terms of $\varepsilon$, a set of three differential problems for $w_{0}$, $w_{1}$ and $w_{2}$, respectively are obtained, which in turn lead (see ^[9,25]) to the homogenized equilibrium equation for the so called macroscopic displacement $w_{0}$:

where $\nabla_{x} w_{0}$ is a macroscopic shear strain, and

is the effective constitutive tensor. In the above, the superscript $^ {\rm t}$ stands for a transpose, ${ {\rm d}} a$ is the area element on ${\mathsf{D}}$, $|\cdot|$ the Lebesgue measure. Function $\boldsymbol{\chi}(y)$ which is introduced as the cell function has components $\chi_s, s = 1, 2$, which represents the unique, null average, ${\mathsf{D}}$-periodic solutions of the cell problem:

where $\boldsymbol{e}_{s}$ is the unit vector parallel to the $y_{s}$ axis.

Using Gauss-Green Lemma and introducing an auxiliary cell function:

the following expression for the effective material moduli is obtained:

with ${ {\rm d}} l$ line element on $\Gamma_j$. Equation (2.14) is indeed applied to compute $\boldsymbol{G} ^{\#}$ in terms of $\tilde{ \boldsymbol{\chi}}$.

2.2.   Weak form

The cell problems (2.10)–(2.12) are usually rephrased in weak form through virtual work fundamental identity:

where $\tilde{{\bf{V}}}: = H^1_{sp}({\mathsf{D}})$ is the space of admissible auxiliary cell functions $\tilde {\boldsymbol{\chi}}$ which are RUC-shifted ${\mathsf{D}}$-periodic vector valued functions satisfying: Eq (2.13)

In a functional space setting:

Indicating ${{\bf{V}}}$ the space of the admissible ${\mathsf{D}}$-periodic variations of $\tilde{{\bf{V}}}$, the bilinear form associated with the stress divergence term results:

which, applying Gauss-Green lemma, using the interface elastic law (2.12) and observing the outward normals to fibre and matrix domains at their interface are mutually opposite, becomes:

or, in a more general fashion:

The form $a(\cdot, \cdot)$ is symmetric, continuous and coercive on $\tilde{\bf{V}}$, so that problem (2.15) is well posed.

3.   Lowest order augmented curvilinear virtual element method

A virtual element discretization of problem (2.15) for curvilinear polygons is here presented, much along the general idea outlined in ^[8]. Denote ${\mathcal{T}}_h$ as a simple polygonal mesh on ${\mathsf{D}}$, i.e., any decomposition of ${\mathsf{D}}$ in a finite set of simple polygons $E$, without holes and with boundary given by a finite number of edges. In this realm, and for the forthcoming application, we will primarily focus on meshes of quadrilaterals, given their simplicity, efficiency and wide possibility of generation through conventional mesh generators. For a curvilinear element with one or more edges lying on a fibre/matrix interface $\Gamma_j$, we shall describe its geometry in exact form as long as it can be recovered by a NURBS representation (cf. Section 4). Hence, the proposed methodology offers a wide feasibility in terms of fibre cross section boundary shapes allowing for complex profiles in this regard. Any given interface $\Gamma_j$ is then intended as a NURBS parametrization i.e., an invertible $C^1$ mapping with proper weights and control points or a $C^1$ blend of NURBS curves such that any portion is explicitly computable. Dropping index $j$ to lighten up notation, we simply indicate

to indicate a generic curved part of the fibre/matrix interface $\Gamma$ with the relevant NURBS representation.

Remark 3.1. In what follows we denote with $e$ any edge of the mesh and with $\nu$ a generic vertex. The symbol $h$ will be associated with a length quantity, hence $h_E$ denotes the diameter of element $E$ and $h_e$ the length of a (possibly curved) edge $e$. As usual, the maximum mesh element size is indicated by $h$ with no subscripts.

3.1.   The augmented virtual element space

For the discretization of the cell problem which is required to compute the effective shear moduli, we propose a simple, computationally efficient enrichment of the VEM space proposed in ^[8], for lowest order discretization $k = 1$, coupled with a NURBS like representation of curvilinear polygons abutting the material interfaces between fibres and matrix. Given a (curvilinear) polygon $E \in {\mathcal{T}}_h$ with some edge laying on a curved interface $\Gamma_j$ ($j \in \{1, 2, .., F\}$), for any such curved edge $e$, we denote with ${{\boldsymbol \gamma}}_e : [a, b] \rightarrow e$ the restriction of the parametrization for $\Gamma_j$ to edge $e$. Then we indicate the standard space of $\mathbb{R}^2$ polynomials of degree $k$ restrictions on edge $e$ as

The local virtual element space on $E$ is then defined introducing the space of traces on a curved edge as follows. For a given integer $k \geq 1$, on a given element $E$ with a curved edge $\gamma$, we consider the trace space

Then, for every integer $k \geq 1$ the local virtual element space is defined as

For $k = 1$, which is the primarily investigated case of the present contribution, the selected degrees of freedom of the local space are (cf. ^[8])

● pointwise evaluation at each vertex of $E$;

● one extra single degree of freedom for any given curved edge.

Note that the investigated case $k = 1$ amounts to a number of local degrees of freedom which equates the classical curvilinear case (cf. ^[3]) plus only one extra dof for the curved edge on the polygon. Figure 2 illustrates the choice of the additional node on the curved edge, associated to the extra degree of freedom in the augmented form for the relevant case $k = 1$ which simply amounts to adding an extra node for any given curvilinear element edge. In passing we note that, owing to interelement continuity, this node is shared by both curved elements abutting such an edge.

The global conforming space is obtained by a standard identification of degrees of freedom, i.e., as the unique values at the interelement, gluing local spaces with $C^0$ regularity:

The same holds for the spaces of discrete variations:

3.2.   Discretization of the problem

The discretization of the problem is an augmented combination of the scheme proposed in ^[8] for the case with straight and curvilinear edges and the original method introduced in ^[7] for curvilinear VEM; the projection space for virtual cell functions is here augmented with an additional deformation mode at no extra cost in terms of local degrees of freedom.

We propose here an augmented projection operator, local to element $E$, which represents an approximated cell function component $\chi_s$ ($s = 1, 2$) which is explicitly computable. Denote with $(\tilde x, \tilde y)$ the usual shifted centroidal Cartesian coordinates at element level. Let $[{{\mathcal P}}_{{\rm aug}}(E)]$ be the set of polynomials spanned by the monomials $\{1, \tilde x, \tilde y, \tilde x\tilde y\}$ (hence linears augmented by the skew-symmetric monomial $\tilde x\tilde y$) restricted to $E$. Given $E \in {\mathcal{T}}_h$ and any $v_h \in {\bf{V}}_h(E)$, the operator $\Pi^{\nabla} : {\bf{V}}_h(E) \rightarrow [{{\mathcal P}}_{{\rm aug}}(E)]$ is defined by

where $\nabla(v_h)$ denotes the gradient of $v_h$ with respect to $y_s$ ($s = 1, 2$) at the microscale. Hence, the augmented projector $\Pi^{\nabla}(v_h)$ is an extension of the standard $L^2$ projection of $v_h$ on $[{{\mathcal P}}_{1}(E)]$ since, at the price of the same number of local degrees of freedom of the standard $k = 1$ case, the approximated cell function belongs to a polynomial space which is larger than the standard one, encompassing also the quadratic skew-symmetric monomial term $\tilde x \tilde y$. This case is indeed quite relevant from the point of view of accuracy, since many RUC geometrical and material arrangements lead to a skew-symmetric solution of the homogenized equilibrium problem in terms of cell function $\chi_s$. It is immediate to check that the above operator $\Pi^{\nabla}(v_h)$ is readily computable through integration by parts and using the adopted degrees of freedom, see ^[5,6,8] for the relevant standard derivations.

The aforementioned integrals and all the quantities relevant to post processing of the solution (namely the homogenized material moduli) are computed for NURBS curvilinear polygons with the efficient quadrature formula outlined in Section 4.

Once the projector is defined, the development of the VEM method proceeds by deriving the local discrete counterpart of the bulk term in the bilinear form (2.19) as follows. For an $E \in {\mathcal{T}}_h$, for all $v_h, w_h \in {\bf{V}}_h(E)$, we define:

where the first term is a direct approximation of $\int_E \nabla w_h \cdot \boldsymbol{G} (\nabla v_h)$ by substituting $\nabla$ with $\nabla \Pi^{\nabla}$, and the second term is a stabilization term of the $dof_i-dof_i$ type (see ^[22] for a thorough discussion of the subject and possible variants). To this end, an additional operator $\pi:{\bf{V}}_h(E) \rightarrow {{\mathcal P}}_1(E)$ on linear monomials is introduced (given the relevant simplicity in coding), defined as the unique minimizer of the euclidean norm distance of the degrees of freedom values with respect to such polynomial space canonical basis. The stabilization form is then taken as:

where the ${\rm dof}_i$ symbol denotes evaluation at the $i^{th}$ local degree of freedom and the positive scalar $\alpha_E$ is introduced to ensure physical dimension consistency with the previous term. The present choice is indeed $\alpha_E = {\rm trace}(\boldsymbol{G})/2$. More details on the stabilization can be found for instance in ^[22].

The global discrete equilibrium equation at microscale over a RUC domain then reads:

for any $v_h, w_h$ in $\tilde{{\bf{V}}}_h$ or ${\bf{V}}_h$. Note that the jumps in the above expression do not imply any computational issue since the virtual cell function and their variations are explicit at the element boundaries.

The proposed VEM then reads

The above construction and the ensuing post-processing heavily relies on integrating known functions over the polygonal domain which can present NURBS-like curved edges. In the following section we detail the proposed efficient quadrature formula adopted throughout our numerical simulation campaign.

4.   Numerical cubature over rational spline curvilinear polygons

In this section we sketch the details on the computation of numerical rules over piecewise rational spline curvilinear polygons, following the details introduced in ^[31]. This formulation includes the domains studied in ^[4]. Such a problem was also considered in ^[26] in the framework of FEM as well as in ^[35].

In particular we consider Jordan domains ${\cal{S}}: \subset {\mathbb{R}}^2$:

$(1)$ Whose boundary $\partial {\cal{S}}$ is described by parametric equations $x = {\tilde{x}}(t)$, $y = {\tilde{y}}(t)$, $t \in [a, b]$, ${\tilde{x}}, {\tilde{y}} \in C([a, b])$, ${\tilde{x}}(a) = {\tilde{x}}(b)$ and ${\tilde{y}}(a) = {\tilde{y}}(b)$;

$(2)$ For which there are partitions $\{I^{(k)}\}_{k = 1, \ldots, M}$ of $[a, b]$, and $\{I_j^{(k)}\}_{j = 1, \ldots, m_k}$ of each $I^{(k)}\equiv [t^{(k)}, t^{(k+1)}]$, such that the restrictions of ${\tilde{x}}$, ${\tilde{y}}$ on each closed interval $I^{(k)}$ are rational splines, w.r.t. the subintervals $\{I_j^{(k)}\}_{j = 1, \ldots, m_k}$.

We adopt as notation

being $u_{k, 1}$, $u_{k, 2}$, $v_{k, 1}$, $v_{k, 2}$ splines on $I^{(k)}$, sharing the same knots and having degree, respectively, $\eta_{k, 1}$, $\eta_{k, 2}$, $\delta_{k, 1}$, $\delta_{k, 2}$, $k = 1, \ldots, M$.

Notice, that since ${\tilde{x}}$, ${\tilde{y}} \in C([a, b])$, we are assuming that the denominators $v_{k, 1}$, $v_{k, 2}$, $k = 1, \ldots, M$ are everywhere not null in the closed interval $I^{(k)}$.

In what follows we intend to show some examples of domains that fulfill these requests. First, this is the case of a spline curvilinear region as those presented in ^[32]. Indeed, let $V_k = ({\tilde{x}}(t_k), {\tilde{y}}(t_k)) \in {\mathbb{R}}^2$, $k = 1, \ldots, M$, $V_{M+1} = V_1$, be the vertices of such a Jordan domain ${\cal S}$, then $\partial {\cal S}: = \cup_{k = 1}^M V_k \frown V_{k+1}$ and each curvilinear side $V_k \frown V_{k+1}$ can be tracked by a parametric spline of degree $\delta_k$, interpolating an ordered subsequence of knots $P_{1, k} = V_k, P_{2, k}, \ldots, P_{m_k-1, k}, P_{m_k, k} = V_{k+1}$ with a suitable parametrization determining each $I_j^{(k)}$ (and thus each $I^{(k)}$).

Another large family of domains $\cal S$ belonging to the class defined above is when $\partial \cal S$ is a composite Bezier closed curve, whose $k$-th component is of the form

where ${\tilde{t}} = \frac{t^{(k+1)} + t^{(k)}}{2}+\frac{t^{(k+1)} - t^{(k)}}{2} t: = \omega_k(t), t \in [0, 1]$ and

are the Bernstein polynomials.

The framework also includes domains whose boundary $\partial \cal S$ is locally a $p$-th degree NURBS curve [24, p.117], where the $k$-th curvilinear side $V_k \frown V_{k+1}$ takes the form

where

● $\{P_{i, k}\}_{i = 1}^{m_k}$ are the control points,

● $\{\lambda_{i, k}\}_{i = 1}^{m_k}$ are the weights,

● $\{B_{i, p}\}_{i = 1}^{m_k}$ are the $p$-th degree B-spline basis functions [12, p.87] defined on the nonperiodic (and nonuniform) knot vector

with $t_{p+j}^{(k)} \leq t_{p+j+1}^{(k)}$, $j = 1, \ldots, m_k-1$.

Some classical examples are domains whose boundary is given by a polygon in which a side is substituted by a circular or ellipse arc (see Figure 3).

5.   The Tchakaloff-Davis-Wilhelmsen approach

In this paper we consider a technique based on an approach that dates back to Davis (1967) and Wilhelmsen (1976) in its general formulation (cf. ^[11,34]), giving a constructive proof of the well-known Tchakaloff existence theorem for positive quadratures (1957), cf. ^[33].

Our strategy is based on a theorem in ^[34] that says the following.

Theorem 5.1. Let be $\mathcal{F}\subset C({\cal{S}})$ a function space of dimension $k$ on a multidimensional compact set ${\cal{S}}\subset \mathbb{R}^d$ (such that $\mathcal{F}$ contains functions that do not vanish on ${\cal{S}}$), $X = \{x_i\}_{i\geq 1}\subset {\cal{S}}$ an everywhere dense point sequence, and $L$ a strictly positive linear functional on $\mathcal{F}$, i.e., $L(f) > 0$ for every $f\in \mathcal{F}$ not vanishing everywhere in ${\cal{S}}$.

Then, for sufficiently large $m$, the set $X_m = \{x_1, ..., x_m\}$ is a "Tchakaloff set" in ${\cal{S}}$, which means that, for every $f\in \mathcal{F}$, $L(f)$ can be represented as the integral with respect to a discrete positive measure with finite support in $X_m$ of cardinality not exceeding $k$ (i.e., as a linear combination with positive coefficients of at most $k$ values of $f$ in $X_m$).

Direct applications of the previous theorem regard the case in which $L(f)$ is an integral functional with respect to an absolutely continuous measure $w(x)dx$, i.e.,

as well as w.r.t. discrete measure with support of large cardinality (for example, a positive algebraic quadrature formula or a QMC formula),

where $Z_M = \{z_1, \dots, z_M\}\subset {\cal{S}}$ with $M > k$.

In both instances, Theorem 5.1 ensures that for sufficiently large $m$ there exist nodes $\{\xi_1, \dots, \xi_\nu\}\subset X_m\subset {\cal{S}}$ and corresponding positive weights $\{w_1, \dots, w_\nu\}$ such that

i.e., $X_m$ is a Tchakaloff set.

In our examples, having in mind to determine algebraic rules with a fixed degree of exactness $n$, we set $\mathcal{F} = \mathbb{P}_n({\cal{S}})$, the space of multivariate polynomials of total-degree not exceeding $n$ on ${\cal{S}}$. We also denote by $N = (n+1)(n+2)/2$ the dimension of the vector space $\mathbb{P}_n({\cal{S}})$.

Aiming to a practical implementation of Tchakaloff-Davis-Wilhelmsen theorem, we proceed as follows:

TDW measure compression algorithm: via moment-matching on the polynomial space $\mathbb{P}_n({\cal{S}})$

$(i)$ Set the moment residual tolerance $tol$, the starting cardinality $m$ and the cardinality increase factor $\theta > 1$;

$(ii)$ Select a basis $\{ \varphi_1, \dots, \varphi_N\}$ of $\mathbb{P}_n({\cal{S}})$ and compute the moments $b_j: = \int_{\cal{S}} \varphi_j dx$, $j = 1, ..., N$;

$(iii)$ Generate a (quasi)-uniformly distributed sequence $X_m = \{x_1, \dots, x_m\}$ in ${\cal{S}}$ (via an implementation of an in-domain routine over ${\cal{S}}$);

$(iv)$ Compute the Vandermonde-like matrix

$(v)$ Compute the factorization $V = QR$ with $Q\in \mathbb{R}^{m\times N}$, $R\in \mathbb{R}^{N\times N}$, and the modified moments $\beta$ such that $R^t\beta = b$;

$(vi)$ Solve the underdetermined system $Q^t u = \beta$ as a Non-Negative Least Squares problem

by Lawson-Hanson active-set method (see, e.g., ^[20] and its Matlab implementation via ${\tt{lsqnonneg}}$, as well as its alternatives ^{[14,15,16,27]}; for its application to numerical quadrature see, e.g., ^[36]);

$(vii)$ If $\|V^t u^\ast - b\|_2 > tol$ then goto $(iii)$ with $m: = \theta m$;

$(viii)$ Select the active weights and nodes:

A careful read of the algorithm shows that its implementation requires some key ingredients.

$(1)$ In ($ii$), once a suitable basis $\{ \varphi_1, \dots, \varphi_N\}$ of $\mathbb{P}_n({\cal{S}})$ is considered, we have to determine the moments $b_j: = \int_{\cal{S}} \varphi_j dx$, $j = 1, ..., N$;

$(2)$ In ($iii$), it is necessary to develop an in-domain routine for Jordan domains ${\cal{S}}$ taken in consideration.

We will develop these demands in the next sections, showing in Figure 4 the results obtained for two not trivial geometries. For practical aspects about the choice of $m$ and $\theta$ in (i) see ^[31]. The Matlab routines are available open-source at ^[29].

6.   On the moments computation

In the TDW measure compression algorithm, we have shown that once a suitable basis $\{ \varphi_1, \dots, \varphi_N\}$ of $\mathbb{P}_n({\cal{S}})$ is considered, we have to determine its moments $b_j: = \int_{\cal{S}} \varphi_j dx$, $j = 1, ..., N$.

In this section we develop these details, following the lines introduced in ^[32] and generalised in ^[31] to the Jordan domains mentioned above.

To this purpose, let ${\cal{R}}^* = [a_1, b_1] \times [a_2, b_2]$ be the bounding box of ${\cal S}$, i.e., the smallest cartesian rectangle ${\cal R} \supseteq {\cal S}$ with sides parallel to the axes. As polynomial basis $\{ \varphi_j\}_{1\leq j\leq N}$, we adopt the shifted lexicographically ordered total-degree product Chebyshev basis

(where $j$ is the position of $(h_1, h_2)$ in such ordering) where

● $T_h(\cdot) = \cos(h \arccos(\cdot))$ is the $h$-degree Chebyshev polynomial of first kind;

● $\alpha_i(s) = (2s-b_i-a_i)/(b_i-a_i)$, $s \in [a_i, b_i]$, $i = 1, 2$.

The use of this basis comes from the necessity of mitigating the possible extreme ill-conditioning of Vandermonde matrices in the standard monomial basis. By Gauss-Green theorem (see e.g., ^[1]),

where

It is easy to achieve that

If $P_{k, s}: = (\tilde{x}(t^{(k)}_s), \tilde{y}(t^{(k)}_s))$ and $P_{k, s} \frown P_{k, s+1}$ is the arc of $\partial {\cal S}$ joining $P_{k, s}$ with $P_{k, s+1}$, we obtain

The evaluation of the integrals on the r.h.s. of (6.2) is a delicate matter. Since $\tilde{x}$ and $\tilde{y}$ are in general rational functions one can use high order Gauss-Legendre rules ^[18], or adaptive routines as the MATLAB built-in routine ${\tt{integral}}$, or the Extended Rational Fejèr Quadrature Rules proposed in ^[13].

7.   The in-domain routine

In the previous section we have shown that the application of the TDW measure compression algorithm requires an in-domain algorithm for the Jordan domains $\cal S$ taken under consideration, i.e., an algorithm that determines if $P \in {\mathbb{R}}^2$ is inside (or not inside) such a ${\cal{S}}$. This problem has been analysed in ^[30,31] and for the sake of the reader, we sketch its details below.

A well-known strategy is based on the Jordan curve theorem that states that a point $P$ belongs to a Jordan domain ${\cal{S}}$ if and only if, having taken a point $P^* \notin {\cal{S}}$ then the segment ${\overline{P^*P}}$ crosses $\partial {\cal{S}}$ an odd number $c(P)$ of times.

In spite of its simplicity, there can be pathological cases, e.g., when ${\overline{P^*P}}$ crosses a vertex or when includes a segment of $\partial \Omega$.

Another problematic situation holds when the boundary $\partial \cal{S}$ has a critical point $S = ({\tilde{x}}(\gamma), {\tilde{y}}(\gamma))$ where

i.e., there is locally a vertical turn of boundary from left to right (or conversely from right to left). If we consider vertical segments ${\overline{P^*P}}$ then the Jordan theorem cannot be directly applied.

These special cases, illustrated in Figure 5, after some effort, can be classified by an algorithm, and thus we start from the most common situation that ${\overline{P^*P}}$ does not contain any critical point or vertical side.

Under these assumptions, let $\cal{R}$ be a rectangle, often called box, with sides parallel to the cartesian axes, whose interior contains $\cal S$, and suppose that we have to establish if $P = (P_x, P_y)$ belongs to $\cal S$.

Let $P^* = (P^*_x, P^*_y)$ be the point in $\cal{R}$ not internal to $\cal{S}$, such that $P^*_x = P_x$ and $P^*_y < P_y$. Geometrically it means that $P^*$ is not internal to $\cal{S}$, shares the same abscissa of $P$, but is vertically below $P$.

In what follows we compute the crossing number $c(P)$, i.e., the number of times in which the vertical segment ${\overline{P^*P}}$ crosses $\partial {\cal{S}}$.

First, we cover $\partial {\cal{S}}$ with the union of possibly overlapping rectangles, each one containing a portion of $\partial {\cal{S}}$ that has no vertical turning points and is parametrized by two rational functions, i.e., locally $({\tilde{x}}(\gamma)$, ${\tilde{y}}(\gamma))$ are the ratio of two polynomials (see Figure 6).

To this purpose, we observe that since ${\tilde{x}}$, ${\tilde{y}}$ are rational splines in each element $I^{(k)}$, $k = 1, \ldots, M$, that takes part to the partitioning of $[a, b]$, then there are

where the restriction of ${\tilde{x}}$, ${\tilde{y}}$ to $I_j^{(k)}$ are rational functions, i.e.,

with $u_{k, j, 1}$, $u_{k, j, 2}$, $v_{k, j, 1}$, $v_{k, j, 2}$ polynomials on the interval $I_j^{(k)}$, with degree, respectively, $\eta_{k, 1}$, $\eta_{k, 2}$, $\delta_{k, 1}$, $\delta_{k, 2}$ (notice that they do not depend on $j$ but just on the local degree of the splines $u_{k, 1}$, $u_{k, 2}$, $v_{k, 1}$, $v_{k, 2}$ in $I^{(k)}$).

If ${\tilde{x}} \equiv c$ in

we set

In other words, ${\cal{B}}_1^{(j, k)}$ is a box that actually consists of a vertical segment.

If ${\tilde{x}}^{\prime}$ has variable sign in $(t^{(k)}_j, t^{(k)}_{j+1})$, let ${\cal{N}}_j^{(k)} = \{t^{(j, k)}_{i}\}_{i = 1, \ldots, l_{j, k}}$ the set of $t^{(j, k)}_{i} \in (t^{(k)}_j, t^{(k)}_{j+1})$ such that ${\tilde{x}}^{\prime}(t^{(j, k)}_{i}) = 0$ (as observed before, the restriction of ${\tilde{x}}$ to $I^{(j, k)}$ is a rational function with the denominator nowhere null, and consequently ${\tilde{x}}^{\prime}$ exists), otherwise let ${\cal{N}}_j^{(k)} = \emptyset$. Next, set

where we suppose that its elements, say $T^{(j, k)}_i$, are in increasing order. Since

from

and $v^2_{k, j, 1}(t) \neq 0$ for each $t \in I_j^{(k)}$, we have that ${\tilde{x}}^{\prime}(t) = 0$ if and only if

and consequently ${\cal{N}}_j^{(k)}$ is available just solving a polynomial equation of degree $\eta_{k, 1}+\delta_{k, 1}-1$.

Now define as monotone boxes the rectangles ${\cal{B}}^{(j, k)}_i$,

where

By definition, if ${\cal{N}}_j^{(k)} = \emptyset$, necessarily there is only the monotone box ${\cal{B}}^{(j, k)}_1$. Since $\tilde{y}$ is a rational function in $[T^{(j, k)}_i, T^{(j, k)}_{i+1}]$, the evaluation of

can be explicitly determined computing the derivative of the polynomial ${\tilde{y}}^{\prime}$, its zeros in $[T^{(j, k)}_i, T^{(j, k)}_{i+1}]$ and the evaluation of ${\tilde{y}}$ at $T^{(j, k)}_i$ and $T^{(j, k)}_{i+1}$.

At this point, we have determined $I^{(j, k)}_i$, such that

● the restriction of ${\tilde{x}}$, ${\tilde{y}}$ to each interval $I^{(j, k)}_i \subseteq [a, b]$ are rational functions,

● ${\tilde{x}}$ is a monotone function (with no turning points of $\partial {\cal{S}}$ in the interior of each ${\cal{B}}^{(j, k)}_i$),

and we are ready to apply the crossing theorem to see if $P = (P_x, P_y)$ is inside the Jordan domain ${\cal{S}}$.

Let

be the set that contains all the monotone boxes ${\cal{B}}_l$ such that ${\overline{P^*P}} \cap {\cal{B}}_l \neq \emptyset$, and that consequently are the only ones that may contribute to the evaluation of the crossing number $c(P)$.

Consider one of these monotone boxes

If $P_y > \beta^{(l)}_2$ then the point is above the box, and thus the segment ${\overline{P^*P}}$ surely crosses ${\partial {\cal{S}}}$ once in ${\cal{B}}_l$ and below $P$, due to the monotonicity of $\tilde{x}$ in ${\cal{B}}_l$.

Otherwise, $P \in {\cal{B}}_l$. As assumed before, ${\overline{P^*P}}$ is free of critical points and vertical segments of the boundary, thus ${\cal{B}}_l$ includes a certain portion of $\partial {\cal{S}}$ described parametrically by two rational functions, say ${\tilde{x}}\vert_{{\cal{B}}_l}$, ${\tilde{y}}\vert_{{\cal{B}}_l}$, with arguments in the interval $I\vert_{{\cal{B}}_l} \subseteq [a, b]$, in which ${\tilde{x}}\vert_{{\cal{B}}_l}$ is monotone and such that $P_x \in {\tilde{x}}\vert_{{\cal{B}}_l}(I\vert_{{\cal{B}}_l}).$ This entails that necessarily there is a unique root $t^* \in I\vert_{{\cal{B}}_l}$ of the polynomial equation

Since ${\tilde{x}}\vert_{{\cal{B}}_l}(t) = \frac{ u(t) } { v(t) }$, for suitable polynomials $u$, $v$, then $t^*$ is the unique solution of the polynomial equation $u(t)-P_y \cdot v(t) = 0$ in $I\vert_{{\cal{B}}_l}$.

We also observe that,

● if ${\tilde{y}} (t^*) < P_y$ then the segment ${\overline{P^*P}}$ crosses the boundary ${\partial {\cal{S}}}$ once in the monotone box, below $P$;

● if ${\tilde{y}} (t^*) > P_y$ then the segment ${\overline{P^*P}}$ does not cross the boundary ${\partial {\cal{S}}}$ once in $B$, below $P$;

● if ${\tilde{y}} (t^*) = P_y$ then $P$ is on the boundary $\partial {\cal{S}}$.

After counting all the crossings, we determine whether a point $P$ is inside or not inside the domain ${\cal{S}}$, by means of Jordan theorem. In Figure 7 we illustrate by an example the contribution given by monotone boxes, while in Figure 8 we show the results given by our procedure when applied to two complicated geometries.

When the previous assumptions do not hold, i.e., the vertical segment ${\overline{P^*P}}$ contains a critical point or a portion of a vertical side of $\partial \cal{S}$ (see Figure 5), one can use an algorithm based on the well-known winding theorem to determine if $P$ belongs to $\cal{S}$. To this purpose, one can compute numerically the so called winding number ${{{wind}}}(P, {\tilde{x}}, {\tilde{y}}) \in {\mathbb{Z}}$,

If the quantity ${{{wind}}}(P, {\tilde{x}}, {\tilde{y}})$ is odd then the point belongs to $\cal{S}$ otherwise is not inside such domain.

In our numerical tests, on a cloud of points, this approach is slower than the evaluation of the crossing numbers. Our experience is that for general domains the winding number strategy will only seldomly called by the in-domain routine proposed here.

In ^[30] the authors take into account techniques to make the implementation efficient and safe from numerical issues. The numerical codes regarding this in-domain routine have been implemented in Matlab and are freely available at ^[28].

8.   Numerical tests

This section presents numerical tests to validate the proposed curvilinear VEM methodology for the homogenization of fibre-reinforced composite materials. In particular, in Section 8.1 we show a set of simple patch tests on a fictitious square domain with a minimal mesh comprising variously distorted curvilinear polygons. In Section 8.2 we perform asymptotic homogenization for the basic case of doubly periodic composite materials with complex fibre shapes. Finally, in Section 8.3 we address a real scale engineering problem i.e., computational homogenization of composite materials with randomly distributed highly complex unidirectional fibres, applying Monte Carlo simulations to obtain statistically averaged effective properties.

8.1.   Patch test

In this section we assess accuracy and robustness through representative boundary value problems on a unit square domain with known solution i.e., patch tests with a linear solution in terms of the cell function $\chi$. The domain is represented by a unit square and meshes with the lowest number of curvilinear quadrilaterals both convex and concave (see Figures 9 and 10) which are also progressively distorted by moving the inner curvilinear polygon which is characterized by having the whole boundary as a sequence of circular arches. The problem is split in two cases for each examined mesh configuration, namely applying a unit shear strain $\boldsymbol{e}_s$ ($s = 1, 2$) along the $y_1$ and $y_2$ direction (see Eq (2.10)), respectively, and solving the ensuing homogenized equilibrium equation, Eq (3.5), for each loading condition.

Numerical results in terms of $H^1$ error measure on the cell function are reported in Table 1, which shows that the proposed methodology exactly solve a general ${{\mathcal P}}_{1}$-type patch test with no major sensitivity to mesh distortion.

8.2.   Doubly periodic composites

In order to verify the capability of the method to tackle complex curvilinear fibre shape in efficient manner, we here study doubly periodic fibre-reinforced composites for different fibre arrangements and material setups. A given doubly periodic composite unit cell is identified through the geometrical features $\varphi$, $\kappa = L_2/L_1$, $f$, and the following dimensionless material parameters:

● fibre/matrix stiffness ratio (contrast factor) $\xi = G ^{ {\rm f}}/G ^{ {\rm m}}$;

● dimensionless interface parameter $\delta = D/(G ^{ {\rm m}} L_1)$;

The simulations refer to square RUC for simplicity i.e., $\varphi = \pi/2$, and $\kappa = 1$. The first benchmark corresponds to elliptic cross section inclusions for $f = 0.1, 0.2, 0.4$, respectively, as can be see in Figure 11, perfect interfaces i.e., $\delta = \infty$ and a contrast factor $\xi = 50$. The inclusions scale homothetically, with a constant ratio between semi-axes fixed at $2/3$. For each computational domain geometry, quadrilateral mesh discretizations are still adopted. Representative meshes of a quarter of the domain are portrayed in Figure 11, in view of the double symmetry offered by all three cases. Presented results are obtained for five mesh sizes applying uniform refinement. Lacking a closed form or semi-analytical solution for the case under investigation, as reference results we use standard Lagrangian quadrilateral quadratic finite elements with an overkilling mesh resolution.

In Figure 12 we report $h-$convergence plots for the cell function $\boldsymbol{\chi}(y)$ in the $H^1-$error norm for uniform mesh refinement, for the set of examined volume fractions.

The expected convergence rate is linear which is obtained for all examined patterns. We notice that the exact geometric representation of the curved interface produce accurate results in conjunction with the implemented NURBS quadrature formula. It is moreover observed that progressively higher volume fractions of fibres require higher computational cost to reach a given accuracy level. The second doubly periodic composite under investigation still regards a complex fibre boundary which, differing from the previous case, exhibits a re-entrant sharp corner. The RUC domain is represented in Figure 13 indicating still a square RUC which lodges a bilobe-type inclusion made of two circular fibres with the same radius fused together with an overlapping of a third of the radial length. Due to symmetry a quarter of the domain is meshed as can be inspected in Figure 13. We examine two sets of problems for a fixed volume fraction $f = 0.2$. First, we consider perfect interfaces ($\delta = \infty$) and values of stiffness ratio $\xi = 10,100, 1000$, respectively. Subsequently, we fix $\xi = 100$ and analyze the following values for the interface integrity indicator: $\delta = 10,100, 1000, \infty$.

Optimal (linear) convergence rate for all cases is observed in Figure 14 where $H^1-$error plots are obtained through uniform mesh refinement sequence. This further set of results confirms method consistency and capability of exactly representing the curvilinear interface geometry. An additional improvement which would furtherly enhance the computational performance of the methodology is adaptive mesh refinement, guided by an efficient and reliable error estimator procedure, which in principle would indicate localization of the error (hence local refinement) at fibre/matrix interface, and at unit cell boundary edges (see ^[3] for a thorough presentation of the method).

8.3.   Random composites with stocastically distributed inclusions

The engineering relevant case of random composites relies on a statistical homogenization approach which for instance amounts to solving for large number of RUC random realizations and subsequently infer homogenized material properties through statistical averaging of the results ^[17,23]. In this context, one key assumption is to consider statistically homogeneous randomly generated microstructures, yielding an isotropic effective behaviour ^[17]. It is known that with this aim, quantitative estimation of the RUC size plays an important role from accuracy and computational cost point of view, since the effective modulus $G ^{\#}$, obtained by Eq (2.14) is actually a random variable depending on the realizations of the RUC ${\mathsf{D}}$ examined throughout the statistical homogenization procedure. The main problem is then to determine a proper RUC size to guarantee a prescribed accuracy on the effective material moduli. The idea is then not to use too large RUCs requiring heavy computational effort, instead, sampling a higher number of smaller RUCs setups to get a prefixed accuracy ^[17]. To this end, in this context where we focus specifically on solving complex curvilinear fibre cross sections, we apply the statistical homogenization procedure presented recently in ^[3] in a simplified form, i.e., with standard uniform mesh refinement.

As a manufactured benchmark we present numerical simulations on four geometrical arrangements of square RUCs with equal, homogeneous isotropic trilobe-shaped fibres with volume fraction $f = 0.2$, stiffness ratio $\xi = 1000$, $\delta = \infty$. The RUC-to-fibre number ratio ranges from $2$ to $16$, see Figure 15 for a pictorial representation of such four typical realizations and relevant quadrilateral meshes. Fibre cross section profiles are no more axis-symmetric, and are obtained blending three elliptical arcs with as many circular fillets granting smoothness of the overall boundary. Each ellipse has semi-axes ratio $2/5$, while the ratio between circular arcs radius and minimal elliptical semi-axis is fixed in order to match the prescribed volume fraction of the composite.

The idea is then to solve a sufficiently large number $n$ of realizations with a standard uniform mesh refinement strategy and numerically assess the statistics of result distribution. We refer to ^[3] for a more detailed description of the computational homogenization procedure based on Monte Carlo simulations.

Figure 16 shows the normalized mean value $\mu_{G ^{\#}} / \Gamma^m$ and the dispersion of $G ^{\#}$ as a function of the number of trilobe inclusions per RUC. As it is expected, larger RUCs are associated to higher accuracy and lower dispersion and hence can be used to obtain a fair estimate of $\mu_{G ^{\#}}$.

9.   Conclusions

In this communication, we have addressed the computational asymptotic homogenization of reinforced composite materials with long parallel fibres having a complex curvilinear cross section. An augmented VEM formulation, based on the conjoined use of an efficient quadrature formula for NURBS-type curvilinear polygons and an augmented projection space, has been developed for the effective computational analysis of the above class of heterogeneous materials in the framework of antiplane deformation. VEM has been recently developed as a generalization of the FEM and it allows the use of curvilinear polygonal elements of general, including non-convex elements. The proposed method has been validated through en extensive numerical campaign showing its generality for modelling accurately multiphase complex material. Further research aims at extending the method to include material non-linear behaviour and damage.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Work partially supported by the DOR funds and the biennial project BIRD 192932 of the University of Padova, by the INdAM-GNCS 2022 Project "Methods and software for multivariate integral models" and by the INdAM-GNCS 2023 Project "Approximation and multivariate integration, with application to integral equations". This research has been accomplished within the RITA "Research ITalian network on Approximation", the "ANA&A SIMAI" group and within the UMI Group TAA "Approximation Theory and Applications" (A. Sommariva).

Conflict of interest

The authors declare no conflicts of interest.