Numerical upscaling of discrete network models

Abstract

In this paper a numerical multiscale method for discrete networks is presented. The method gives an accurate coarse scale representation of the full network by solving sub-network problems. The method is used to solve problems with highly varying connectivity or random network structure, showing optimal order convergence rates with respect to the mesh size of the coarse representation. Moreover, a network model for paper-based materials is presented. The numerical multiscale method is applied to solve problems governed by the presented network model.

1 Introduction

Network structures are used to model a wide variety of phenomena, such as flow in porous media, traffic flows, elasticity of materials, body deformation in computer graphics, molecular dynamics, and fiber materials. In these applications, the microscale behaviour determines the macroscale properties of the system. Often a full microscale model is difficult or impossible to work with because of the vast computational complexity. Therefore, there is an interest in constructing coarser, but still accurate, representations of the entire system. Such a procedure is sometimes referred to as upscaling or homogenization. In this work a numerical upscaling method for discrete networks is presented.

There exist several numerical upscaling methods for partial differential equations (PDE) based on the idea of homogenization, such as the heterogeneous multiscale method (HMM) [20], the multiscale finite element method (MsFEM) [8], and the more recent works [3, 17]. The upscaling approach presented in this paper is based on the localized orthogonal decomposition method (LOD) [4, 15], which in turn is inspired by the variational multiscale method (VMM) [9]. Multiscale methods applied to network problems are for instance investigated by Ewing [5] and Ilev et al. [11] who study the heat conductivity of network materials and develop an upscaling method by solving the heat equation locally over small sub-domains. These local solutions are used to compute an effective global thermal conductivity tensor. Della Rossa et al. [2] investigate network models of traffic flows and derive a governing PDE for the macroscale by formulating traffic flow equations for single network nodes and interpreting the relations as finite difference approximations. The macroscale parameters are resolved using a two-scale averaging technique. Chu et al. [1] develop a multiscale method for networks representing flows in a porous medium. The medium is modelled as a network where nodes represent pores and edges represent throats. The conductance of each throat is assumed to be given by Hagen–Poiseuille equation, and using mass conservation equations for the flow through the network, a model for the microscale is attained.

The numerical upscaling method proposed in this work is developed for general unstructured networks. The network is supposed to represent the microscale, and the macroscale is represented by a finite element mesh which is coarse in comparison to the fine scale network. The coarse grid does not have to be related to the network in any way except that both cover the same computational domain, and therefore the method can be applied to arbitrary network geometries. The coarse FEM grid is used to define a macroscale solution space spanned by basis functions defined at each coarse grid node as in standard FEM. The upscaling idea is to modify the coarse basis functions to account for the microscale features of the network. This is accomplished by solving local sub-network problems at each coarse basis function. The modified basis functions are thereafter used to solve a global low-dimensional system resulting in an accurate macroscale solution. The method leads to modified basis functions that decay exponentially, and hence localization of the local sub-network problems can be utilized, reducing the computational cost considerably while preserving optimal convergence rates.

Moreover, this paper includes a two-dimensional network model, which can be used to model paper-based materials in form of fiber networks. The macroscale mechanical properties of paper-based materials are of great interest. Paper is a heterogeneous material built up of fibers bonded together into a network structure. The mechanical properties of paper depend primarily on the properties of the fibers and the bonds between them. In [12, 14, 19], computational fluid dynamics and advanced contact modeling are used to simulate the paper forming process. One future aim is to utilize that framework together with the proposed multiscale method to create virtual fiber networks and investigate the macroscale mechanical properties. A network representation including fibers and bonds is a suitable methodology to study the mechanical properties of paper [10, 13, 18]. Moreover, the varying properties of single fibers and bonds, as well as an interest for fracture propagation simulations, call for an upscaling approach. The presented network model is based on forces arising at the nodes when the network is displaced, acting to restore the initial configuration. The network model is similar to lattices models like [16, 21] where edges are represented by springs. Moreover, angle springs between pair of edges are included. A novelty of the network model in this work is a third type of force phenomenon resulting in an effect similar to the Poisson effect. Force equilibrium equations at each node result in a matrix equation which can be very large. For a regular network, the model converges to the linear elasticity equation when the length of the network edges tends to zero. The numerical upscaling method is applied to the network model and numerical examples are solved to demonstrate the convergence rates of the method. The examples show how the proposed numerical upscaling method resolves fine scale features which the standard FEM cannot.

The outline of this text is as follows. In Sect. 2, the general problem formulation is stated. Thereafter, in Sect. 3, the theory of the numerical upscaling method is presented. Sect. 4 contains error analysis, and in Sect. 5, the two-dimensional network model is described. In Sect. 6, numerical examples are presented, showing the convergence rates of the proposed method. Lastly, in Sect. 7, conclusions and future work are discussed.

2 Problem formulation

Consider a problem modelled by a network with properties governed by a connectivity matrix $K\in \mathbb {R}^{n\times n}$. The matrix K can for instance be the discrete Poisson operator describing heat conduction, the finite difference discretization of the linear elasticity operator, or represent a more complex model, such as of the mechanics of a fiber network. Let $F\in \mathbb {R}^{n}$ denote the load vector and let the solution vector be denoted u, belonging to a vector space $V\subset \mathbb {R}^n$. The network problem can be stated in two equivalent ways, either:

$$\begin{aligned} \text {Find }u:\quad {\bar{K}} u = {\bar{F}}, \end{aligned}$$

(2.1)

or:

$$\begin{aligned} \text {Find }u\in V:\quad v^T Ku=v^T F, \quad \forall v\in V. \end{aligned}$$

(2.2)

In the first formulation, (2.1), ${\bar{K}}$ and $\bar{F}$ denotes modifications of K and F by explicitly including the restriction of u to the space V, for instance by holding some nodes fixed. To ensure existence and uniqueness of the second formulation, (2.2), it is assumed that K is symmetric and positive definite on V. A matrix $K\in \mathbb {R}^{n\times n}$ is positive definite on a subset $V\subset \mathbb {R}^n$ if $v^TKv>0$ for all nonzero $v\in V$. Moreover, a symmetric positive definite matrix K constitutes a scalar product $\langle u,v\rangle =u^TKv$ on V, a property that will be used later.

In Fig. 1 three examples of networks are shown. The network in Fig. 1a exemplifies a finite difference grid for the unit square, with K as the resulting discretization of the linear elasticity operator. This problem setup can be used to find the node displacements u under applied node forces F. To attain a solvable system $Ku=F$, some degrees of freedom have to be prescribed, resulting in the restricted solution space V. The network in Fig. 1b can represent a conductive medium, governed by the discrete Poisson equation. The temperature at each node is contained in u. The network in Fig. 1c is a fiber network building up a paper sheet. The fibers are modelled as chains of edges connected at nodes with bonds between fibers at common network nodes.

Fig. 1

Three examples of networks: a regular square network (a), a regular square network with randomly perturbated nodes (b), and a fiber network (c) (generated as in [12])

The objective of this paper is to develop a numerical upscaling method for networks, circumventing the computational issues arising when materials of macrosize are considered. The idea is to reduce the size of the system by introducing a subspace $V_\mathrm{ms}\subset V$, as a coarse representation of the network. This space is called the multiscale space and it should fulfil the condition that $\dim V_\mathrm{ms}$ is much lower than $\dim V$. The multiscale solution is attained from the problem

$$\begin{aligned} \text {Find }u\in V_\mathrm{ms}:\quad v^T Ku_\mathrm{ms}=v^T F, \quad \forall v\in V_\mathrm{ms}. \end{aligned}$$

The aim is to construct a multiscale space such that an error $\Vert u-u_\mathrm{ms}\Vert$ is small. To achieve this, a FE-type coarse space is first introduced, which does not have the desired approximation properties. This coarse space is then modified by solving local sub-networks problems, resulting in the desired multiscale space. In the following section such a numerical homogenization method is presented.

3 Numerical homogenization of networks

Consider a network with N nodes and properties governed by a symmetric and positive semi-definite matrix $K\in \mathbb {R}^{n\times n}$, where $n=d\cdot N$ is the number of degrees of freedom of the network, and d denotes the number of degrees of freedom at each node. For instance, for an elastic network, where node displacements are to be solved, the number of degrees of freedom at each node will be two or three, depending on if the network is two- or three-dimensional. In the following presentation the space is assumed to be two-dimensional, but the method works analogously for three dimensions. Denote by $p_i\in \mathbb {R}^2$ the position of the node corresponding to degree of freedom $i = 1, \dots , n$. Note that groups of d degrees of freedom correspond to the same position.

Let the solution vector be denoted $u\in \mathbb {R}^{n}$. The ordering of nodes and their degrees of freedom is arranged such that if $d=2$, u(1) and u(2) correspond to the first and second degree of freedom of node 1, u(3) and u(4) correspond to the first and second degree of freedom of node 2, and so on, with analogous ordering if d is larger. Here u(i) denotes the i:th component of vector u. Let $F\in \mathbb {R}^{n}$ denote the load vector. The system $Ku=F$ is not necessarily solvable without prescribing some degrees of freedom. Consider fixed constraints with zero displacement (non-zero displacement is treated in Sect. 3.4) and let $\mathscr {N}_D\subset \{1,\dots , n\}$ be the set of indices corresponding to the fixed degrees of freedom. Let $\mathscr {N}=\{1,\dots , n\}{\setminus }\mathscr {N}_D$. Denote by $V \subset \mathbb {R}^{n}$ the restricted solution space defined by

$$\begin{aligned} V = \{v\in \mathbb {R}^n : v(i)=0, \quad i\in \mathscr {N}_D\}. \end{aligned}$$

The variational formulation of the network displacement problem reads:

$$\begin{aligned} \text {Find }u\in V:\quad v^T Ku=v^T F, \quad \forall v\in V. \end{aligned}$$

(3.1)

For the problem to be solvable it is assumed that K in addition to being symmetric, also is positive definite on the restricted solution space V.

3.1 Coarse grid representation

The overall idea of the upscaling method is to introduce a coarse grid, representing the network at the macroscale. See Fig. 2a for an illustration of a network with a coarse grid representation. At each coarse node, d number of basis functions are defined similarly as in the finite element method. These basis functions span a low dimensional solution space which gives an insufficient description of the fine scale features. To include the fine scale information, the basis functions are modified by solving local sub-network systems. Thereafter the modified basis functions are used to solve a global system, smaller than the full system including all nodes, resulting in an upscaled approximation of the original problem. In what follows, the details of this procedure are described.

Fig. 2

Example of a square network with a FEM quadrilateration representation

Let the coarse grid be denoted $\mathscr {T}$, containing M coarse nodes and let $m=d\cdot M$ be the degrees of freedom of the coarse grid. One choice of coarse grid is a quadrilateration as in Fig. 2a. It is assumed that the coarse grid constitutes a good approximation of the computational domain of the network, and that each coarse element contains at least one network node and that $N> M$. Let $\varLambda _i:\mathbb {R}^2\rightarrow \mathbb {R}, \,i=1,\dots , m$, denote the coarse nodal basis functions of the grid $\mathscr {T}$. For a quadrilateration, bilinear basis functions are suitable, illustrated in Fig. 2b.

Let $\mathscr {M}_D\subset \{1,\dots , m\}$ be the set of indices corresponding to fixed coarse degrees of freedom. The fixation of coarse grid nodes is determined from the set of fixed network nodes, $\mathscr {N}_D$. Consider a coarse node with basis function $\varLambda _i$, describing for instance the x-displacement of that node. If there exists a network node with fixed degree of freedom j such that $p_j$ lies in the support of $\varLambda _i$, and j also describes x-displacement, then the coarse degree of freedom i should be fixed. This is illustrated in Fig. 3. The fixation condition is equivalently stated as:

$$\begin{aligned} i\in \mathscr {M}_D \quad \text {if} \quad \exists j\in \mathscr {N}_D :\quad \Big (\varLambda _i(p_j)\ne 0 \quad \text {and}\quad i\equiv j \!\!\!\!\!\!\pmod {d}\Big ). \end{aligned}$$

(3.2)

The modulo operation is used to check if two degrees of freedom, here i and j, correspond to the same displacement direction (e.g. x or y for $d=2$). For networks with more complex boundary geometry, the coarse grid has to be refined at the boundary to attain a proper representation of the fixed boundary conditions.

Fig. 3

Example illustrating the fixation of coarse network nodes for $d=1$. The coarse grid nodes are marked with squares and the network nodes with dots. Fixed nodes are encircled. To the left, fixed network nodes are marked with circles, illustrating the set $\mathscr {N}_D$. Based on $\mathscr {N}_D$ and the condition (3.2), the coarse grid nodes which should be fixed, $\mathscr {M}_D$, have been marked with larger circles in the right plot

Let $\mathscr {M}=\{1, \dots , m\}{\setminus } \mathscr {M}_D$ denote the set of nonprescribed coarse degrees of freedom. The positions of the coarse nodes, $\{P_i\}_{i=1}^m$, defined similarly as the positions of the network nodes, are a subset of $\mathbb {R}^2$, likewise as the nodes of the network. However, these two subsets do not have to be related, but as already noted, it is assumed that each coarse element contains at least one network node.

Next, two vector spaces are introduced, the coarse space $V_H$, and the detail space W. The coarse space is defined from the coarse basis functions in the following way. Let $\lambda _i\in \mathbb {R}^n, i=1,\dots , m$, be the interpolation of the coarse nodal basis functions to the network nodes given by

$$\begin{aligned} \lambda _i(j) = {\left\{ \begin{array}{ll} \varLambda _i(p_j), &{}\quad \text {if} \quad i\equiv j \pmod {d},\\ 0, &{} \quad \text {else}. \end{array}\right. } \end{aligned}$$

See Fig. 2c for an illustration of the interpolated vector $\lambda _i$ of the coarse nodal basis function $\varLambda _i$. Note that $\left( \sum _{i=1}^m \lambda _i\right) (j) = 1, \,\,\forall j=1,\dots ,n$.

The coarse space is defined as the span of the interpolated non-fixed basis functions, that is

$$\begin{aligned} V_H = \text {span}(\{\lambda _i\}_{i\in \mathscr {M}}), \end{aligned}$$

with dimension $\dim V_H = m_H := |\mathscr {M}|$. Note that $\lambda _i$ is defined for all $i=1,\dots , m$, but $\lambda _i\in V$ only if $i\in \mathscr {M}$. Let the matrix $B_H = [\{\lambda _i\}_{i\in \mathscr {M}}]\in \mathbb {R}^{n\times m_H}$ contain the basis vectors of the coarse space $V_H$ as its columns. It is assumed that the columns are linearly independent. The matrix $B_H$ is called the prolongation matrix and acts as a map $B_H:\mathbb {R} ^{m_H}\rightarrow V_H$.

To define the detail space, a restriction matrix $C_H\in \mathbb {R}^{m_H\times n}$ is introduced acting as a map $C_H:\mathbb {R}^n\rightarrow \mathbb {R}^{m_H}$. In this work, the restriction matrix is chosen as $C_H=B_H^T$ (for examples of other choices of restriction operator, see [4]). With $C_H=B_H^T$, an equivalent definition of the coarse space $V_H$ is as the range of the map $B_HC_H:V\rightarrow V_H$, that is

$$\begin{aligned} V_H= \{B_HC_Hv:\,\, v\in V\}. \end{aligned}$$

The detail space is defined as the null space of the restriction matrix:

$$\begin{aligned} W= \{v\in V: \,\,C_Hv = 0\}. \end{aligned}$$

With $C_H=B_H^T$ it holds that $v\in W$ if the bilinear weighted average of v is zero for each interpolated bilinear basis function $\lambda _i$, i.e. $\lambda _i^T v= 0,\,\forall i\in \mathscr {M}$.

The coarse space and the detail space constitute a splitting of V such that each $v\in V$ can be uniquely decomposed as $v= v_H+w$ where $v_H\in V_H$ and $w\in W$. Before proving this fact, a lemma is stated showing the relation between the spaces $\mathbb {R}^{m_H}$, $\mathbb {R}^{n}$ and $V_H$, and the maps in-between, illustrated in Fig. 4.

Fig. 4

A sketch of the two spaces $\mathbb {R}^{m_H}$ and $\mathbb {R}^n$, and the subspace $V_H\subset \mathbb {R}^n$. The four mappings $B_H:\mathbb {R}^{m_H}\rightarrow V_H$, $C_H:\mathbb {R}^{n}\rightarrow \mathbb {R}^{m_H}$, $C_HB_H:\mathbb {R}^{m_H}\rightarrow \mathbb {R}^{m_H}$ and $B_HC_H: \mathbb {R}^{n}\rightarrow V_H$ are also shown

Lemma 3.1

If $B_H$ has linearly independent columns and $C_H=B_H^T$, then for each $v_H\in V_H$ there exists ${\bar{v}}_H\in V_H$ such that $B_HC_H{\bar{v}}_H = v_H$. Moreover, if $v_1, v_2\in V_H$ such that $B_HC_Hv_1 = B_HC_H v_2$, then $v_1=v_2$.

Proof

The map $B_H:\mathbb {R}^{m_H}\rightarrow V_H$ is one-to-one since $B_H$ has linearly independent columns. Moreover, the map $C_HB_H:\mathbb {R}^{m_H}\rightarrow \mathbb {R}^{m_H}$ is invertible since $C_HB_H=B_H^TB_H$ is symmetric and positive definite due to the fact that $x^TC_HB_Hx=|B_Hx|^2 \ge 0$ and $|B_Hx|=0$ implies $x=0$. Given $v_H\in V_H$, it exists $a\in \mathbb {R}^{m_H}$ such that $B_Ha=v_H$, and since $C_HB_H$ is invertible it exists $b\in \mathbb {R}^{m_H}$ such that $C_HB_Hb=a$. Therefore $v_H=B_HC_HB_Hb$ leading to ${\bar{v}}_H = B_Hb\in V_H$.

To prove the second part, assume $v_1\ne v_2$. Then there exists $a_1,a_2\in \mathbb {R}^{m_H}$ with $a_1\ne a_2$ such that $v_1=B_Ha_1$ and $v_2=B_Ha_2$. Since $C_HB_H$ is invertible, $C_HB_Ha_1\ne C_HB_Ha_2$, contradicting the fact that $B_HC_HB_Ha_1= B_HC_HB_Ha_2$, hence $v_1=v_2$. $\square$

Proposition 3.1

If $B_H$ has linearly independent columns and $C_H=B_H^T$, then $V = V_H\oplus W$ uniquely.

Proof

For $v\in V$, let $v_H=B_HC_Hv$. Lemma 3.1 states the existence of ${\tilde{v}}_H\in V_H$ such that $v_H=B_HC_H\tilde{v}_H$. Since $B_H$ has linearly independent columns the relation $0=v_H-v_H=B_HC_Hv-B_HC_H{\tilde{v}}_H = B_H(C_Hv-C_H{\tilde{v}}_H)$ implies that $C_Hv-C_H{\tilde{v}}_H = 0$ with conclusion that $v-\tilde{v}_H\in W$. Therefore $v={\tilde{v}}_H+(v-{\tilde{v}}_H)$ is a desired decomposition. To show uniqueness, consider two decompositions $v=v_{H,1} + w_1$ and $v= v_{H,2}+w_2$. Then $v_{H,1} + w_1= v_{H,2}+w_2$, and applying $B_HC_H$ on both sides gives $B_HC_Hv_{H,1}=B_HC_Hv_{H,2}$. From the last part of Lemma 3.1 it follows that $v_{H,1}=v_{H,2}$. $\square$

Using the detail space W, together with the connectivity matrix K, the multiscale space $V_\mathrm{ms}$ is defined as the K-orthogonal complement of W:

$$\begin{aligned} V_\mathrm{ms} = \{v\in V: w^TKv = 0,\quad \forall w\in W\}. \end{aligned}$$

The spaces W and $V_\mathrm{ms}$ constitute another splitting of V implying that every $v\in V$ can be decomposed uniquely as $v = v_\mathrm{ms}+w$ where $v_\mathrm{ms}\in V_\mathrm{ms}$ and $w\in W$.

Proposition 3.2

Assume $B_H$ has linearly independent columns and $C_H=B_H^T$. If K is symmetric and positive definite on V, then $V = V_\mathrm{ms}\oplus W$ uniquely.

Proof

Consider $v\in V$. From Proposition 3.1 it is known that $v=v_H+ {\tilde{w}}$ with $v_H\in V_H$ and ${\tilde{w}}\in W$. Let $z\in W: \, w^TKz=w^TKv_H,\, \forall w\in W$, which has a unique solution since K is symmetric and positive definite on V. Define $v_\mathrm{ms}=v_H-z$ and $w={\tilde{w}} + z$, where the second sum is in W. Since $x^TKv_\mathrm{ms}=x^TKv_H - x^TKz=0, \,\forall x\in W$, it is true that $v_\mathrm{ms}\in V_\mathrm{ms}$, giving the desired decomposition as $v = v_\mathrm{ms} + w$. To prove uniqueness, consider $v=v_{\text {ms}, 1} + w_ 1$ and $v=v_{\text {ms}, 2} + w_ 2$. Then $0=x^TK(v-v) = x^TK(w_1 + v_{\text {ms}, 1} - w_2 - v_{\text {ms}, 2})= x^TK(w_1-w_2), \, \forall x\in W$, implying that $w_1=w_2$. $\square$

The multiscale solution, $u_\mathrm{ms}\in V_\mathrm{ms}$, to the original problem (3.1), is defined as the solution to the problem

$$\begin{aligned} \text {Find }u_\mathrm{ms}\in V_\mathrm{ms}:\quad v^TKu_\mathrm{ms} = v^TF,\quad \quad \forall v\in V_\mathrm{ms}. \end{aligned}$$

(3.3)

Proposition 3.3

Let K be symmetric and positive definite on V, and $F\in V$. Then there exists a unique solution to problem (3.3).

Proposition 3.4

Let $u_f\in W$ be such that $w^TKu_f=w^TF, \,\, \forall w \in W$. Then the sum $u = u_\mathrm{ms} + u_f$, where $u_\mathrm{ms}$ is the solution to the multiscale problem (3.3), solves the original problem (3.1).

Proof

Same arguments as used in Proposition 3.3 show that there exists a unique such $u_f$. According to Proposition 3.2, $v\in V$ can be decomposed as $v=v_\mathrm{ms} + w$, where $v_\mathrm{ms}\in V_\mathrm{ms}$ and $w\in W$. Using orthogonality, and that $u_f$ and $u_\mathrm{ms}$ are solutions to their respective problem, it can be derived that

$$\begin{aligned} \begin{aligned} v^TKu&= (v_\mathrm{ms} + w)^TK(u_\mathrm{ms} + u_f) \\&= v_\mathrm{ms}^TKu_\mathrm{ms} + v_\mathrm{ms}^TK u_f+ w^T K u_\mathrm{ms}+ w^T Ku_f\\&= v_\mathrm{ms}^TF +0+ 0+ w^T F\\&= v^TF. \end{aligned} \end{aligned}$$

$\square$

To solve the multiscale problem (3.3), it is convenient to construct a basis for the multiscale space $V_\mathrm{ms}$, which can be used to simplify the problem to a matrix equation. A basis for $V_\mathrm{ms}$ is constructed using the vectors $\lambda _i$, by defining modification vectors $\phi _i\in V$, $i\in \mathscr {M}$, as solutions to the problems

$$\begin{aligned} \phi _i\in W : \quad w^TK(\lambda _i-\phi _i) = 0,\quad \forall w\in W. \end{aligned}$$

(3.4)

Proposition 3.5

If K is symmetric and positive definite on V, $B_H$ has linearly independent columns and $C_H=B_H^T$, then the vectors $\{\lambda _i-\phi _i\}_{i\in \mathscr {M}}$ constitute a basis for $V_\mathrm{ms}$.

Proof

The problem to find $\phi _i$ such that $w^TK\phi _i=w^TK\lambda _i$, $\forall w\in W$, has a unique solution $\phi _i\in W$ since K is symmetric and positive definite on V. By construction, it is also true that $\lambda _i-\phi _i\in V_\mathrm{ms}$. To prove linear independence, consider $\sum a_i (\lambda _i-\phi _i)=0$ and apply $B_HC_H$ to both sides. Using that $C_H\phi _i=0$ gives $B_HC_H\sum a_i\lambda _i = 0 = B_HC_H0$, and by the second part of Lemma 3.1 it follows that $\sum a_i\lambda _i = 0$ implying $a_i=0$. $\square$

Using the above constructed basis for $V_\mathrm{ms}$, $\{\lambda _i - \phi _i\}_{i\in \mathscr {M}}$, to assemble the matrix

$$\begin{aligned} B_\mathrm{ms}= [\{\lambda _i - \phi _i\}_{i\in \mathscr {M}}]\in \mathbb {R}^{n\times m_H}, \end{aligned}$$

reduces the variational form of the multiscale problem (3.3) to the equivalent matrix problem

$$\begin{aligned} B_\mathrm{ms}^TKB_\mathrm{ms}U_\mathrm{ms} = B_\mathrm{ms}^TF, \end{aligned}$$

(3.5)

where $U_\mathrm{ms}\in \mathbb {R}^{m_H}$ and $u_\mathrm{ms} = B_\mathrm{ms}U_\mathrm{ms}$.

At this point, a multiscale space with low dimension compared to the solution space V has been constructed as was desired in the problem formulation. Moreover, it has been shown that the problem can be solved through a matrix equation after constructing a basis for the multiscale space. However, the problems of solving the modified basis functions have the same size as the original problem. It turns out that this can be circumvented, utilizing that the modifications $\phi _i$ decay exponentially, implying that the problems can be localized. This is presented in the following section.

3.2 Localization

The described method requires systems to be solved which are as large as K. However, as will be demonstrated by numerical examples in Sect. 6, the modifications $\phi _i$ decay fast away from its coarse node. Therefore the problems (3.4), of calculating $\phi _i$, can be localized with preserved convergence rates. The localization is accomplished by solving each problem (3.4) on a restricted domain, called patch.

As in FEM, it is natural to assemble the stiffness matrix elementwise. In this work, it is suitable to assemble the connectivity matrix K over each coarse element E, such that $K=\sum _{E\in \mathscr {T}} K_E$. The local element connectivity matrices $K_E:V\rightarrow V$ are assembled for each element $E\in \mathscr {T}$ by only considering edges in each element. See Fig. 5 for an illustration. For unstructured networks, edges may intersect the elements. Such a situation is resolved by temporarily dividing the edges at intersection points. Using the decomposition of the connectivity matrix into local element matrices, the modifications $\phi _i$ can analogously be assembled as $\sum _{E\in \mathscr {T}} \phi _i^E$, where $\phi _i^E$ is the solution to the problem

$$\begin{aligned} \text {Find }\phi _i^E\in W: \quad w^TK\phi _i^E=w^TK_E\lambda _i,\quad \forall w\in W. \end{aligned}$$

Fig. 5

Element E is marked with a square and the edges which are included in the assemble of $K_E$ is marked with thick lines

Proposition 3.6

The sum of the elementwise modifications, $\sum _{E\in \mathscr {T}} \phi _i^E$, solves the original problem (3.4).

Proof

It follows that

$$\begin{aligned} w^TK\sum _{E\in \mathscr {T}} \phi _i^E = \sum _{E\in \mathscr {T}} w^TK \phi _i^E = \sum _{E\in \mathscr {T}} w^TK_E \lambda _i = w^T\sum _{E\in \mathscr {T}}K_E \lambda _i = w^TK \lambda _i. \end{aligned}$$

$\square$

With elementwise assembly, it is convenient to use patches centred at each element. For an element $E\in \mathscr {T}$, let its patch be denoted $\omega _E\subset \mathbb {R}^2$. One suitable choice of patch geometry is a circle with center coinciding with the element center, as illustrated in Fig. 6a. Another choice is to use a fixed number of layers of coarse elements surrounding the considered element E, as depicted in Fig. 6b. Let the patch size be described by the parameter $\rho$ such that $\rho H$ is the radius of the patch, where H denotes the coarse element size. In the two examples shown in Fig. 6, the value of $\rho$ corresponds to 1.5. Let $\mathscr {N}_E \subset \mathscr {N}$, denote the degrees of freedom of network nodes that are in the patch $\omega _E$. Similarly, let $\mathscr {M}_E \subset \mathscr {M}$, denote the degrees of freedom of coarse nodes that are in the patch of element E.

Fig. 6

Two square networks with different types of patch geometries

For each element E, its localized subspace of the detail space is defined according to

$$\begin{aligned} {\tilde{W}}_E= \{w\in W: \quad w(i) = 0, \quad \forall i\notin \mathscr {N}_E\}. \end{aligned}$$

The localized modification ${{\tilde{\phi }}}_i$ is attained by solving the problems

$$\begin{aligned} \text {Find }{{\tilde{\phi }}}_i^E\in {\tilde{W}}_E: \quad w^TK \tilde{\phi }_i^E=w^TK_E\lambda _i,\quad \forall w\in {\tilde{W}}_E, \end{aligned}$$

for each element E, and taking the sum ${\tilde{\phi }}_i = \sum _{E\in \mathscr {T}} {\tilde{\phi }}_i^E$. The resulting localized multiscale space is denoted ${\tilde{V}}_\mathrm{ms}$, and is defined as

$$\begin{aligned} {\tilde{V}}_\mathrm{ms} =\text {span}(\{\lambda _i-{\tilde{\phi }}_i\}_{i\in \mathscr {M}}). \end{aligned}$$

3.3 Algebraic formulation

In this section the algebraic formulation of the numerical multiscale method is presented. The multiscale method consists of two main steps. First, the basis for the multiscale space $V^\text {ms}$ is constructed by calculating each $\phi _i$ from (3.4). Secondly, the attained modified basis is used to solve the global multiscale problem (3.3). In this section, elementwise assembling and localization, described in Sect. 3.2, is employed.

First, a useful matrix notation is introduced. Consider a matrix $A\in \mathbb {R}^{a\times b}$. Let $\mathscr {A}\subset \{1,2,\dots ,a\}$ and $\mathscr {B}\subset \{1,2,\dots ,b\}$ be subsets of the matrix row and column indices respectively. A new matrix $A(\mathscr {A}, \mathscr {B})\in \mathbb {R}^{|\mathscr {A}|\times |\mathscr {B}|}$ is extracted from A by only considering rows corresponding to indices in $\mathscr {A}$ and columns corresponding to indices in $\mathscr {B}$, that is $A(\mathscr {A}, \mathscr {B})= (a_{ij})_{(i,j)\in \mathscr {A}\times \mathscr {B}}$ .

Using the introduced matrix notation, the following matrices and vectors are defined for $E\in \mathscr {T}$ and $i\in \mathscr {M}$:

$$\begin{aligned} \begin{aligned} K^E&= K(\mathscr {N}_E, \mathscr {N}_E),\\ C_H^E&= C_H(\mathscr {M}_E, \mathscr {N}_E),\\ r^E_i&= K_E(\mathscr {N}_E, :) \lambda _i. \end{aligned} \end{aligned}$$

The multiscale method can now be formulated as solving several matrix systems. For each element $E\in \mathscr {T}$, and each coarse degree of freedom $i\in \mathscr {M}$, the following system is solved

$$\begin{aligned} \begin{bmatrix} K^E&\quad {C_H^E}^T\\ C_H^E&\quad 0 \end{bmatrix} \begin{bmatrix} {\tilde{\varphi }}_i^E \\ \eta _i^E \end{bmatrix} = \begin{bmatrix} r_ i^E \\ 0 \end{bmatrix}. \end{aligned}$$

The correction vectors ${\tilde{\phi }}_i^E$ is attained from ${\tilde{\varphi }}_i^E$ as

$$\begin{aligned} \begin{aligned} {\tilde{\phi }}_i^E(\mathscr {N}_E)&= {\tilde{\varphi }}_i^E, \\ {\tilde{\phi }}_i^E(\mathscr {N}_E^C )&= 0, \end{aligned} \end{aligned}$$

where $\mathscr {N}_E^C=\{1,\dots ,n\}{\setminus }\mathscr {N}_E$. The full modifications are thereafter calculated according to

$$\begin{aligned} {\tilde{\phi }}_i = \sum _{E\in \mathscr {T}} {\tilde{\phi }}_i^E,\quad i\in \mathscr {M}, \end{aligned}$$

and used to assemble the modified basis matrix

$$\begin{aligned} {\tilde{B}}_M= [\{\lambda _i - {\tilde{\phi }}_i\}_{i\in \mathscr {M}}]. \end{aligned}$$

Finally, the following localized version of the global problem (3.5) is solved:

$$\begin{aligned} {\tilde{B}}_M^TK{\tilde{B}}_M {\tilde{U}}_\mathrm{ms} = {\tilde{B}}_M^TF, \end{aligned}$$

(3.6)

where the fine multiscale solution is calculated as $\tilde{u}_\mathrm{ms} = {\tilde{B}}_M{\tilde{U}}_\mathrm{ms}$.

Summarized, the algebraic formulation of the numerical multiscale method is:

3.4 Non-zero fixed boundary conditions

Consider a displacement problem where $F=0$ and the set of prescribed degrees of freedom, $\mathscr {N}_D$, includes degrees with non-zero displacement. Let

$$\begin{aligned} V_D = \{u\in V: u_i=g_i, \quad i\in \mathscr {N}_D\}, \end{aligned}$$

where $g_i$ are the prescribed displacement of degree of freedom $i\in \mathscr {N}_D$. Let still $V = \{u\in V: u_i=0, \quad i\in \mathscr {N}_D\}$. The displacement problem is

$$\begin{aligned} \text {Find } u\in V_D: \quad v^T Ku=0, \quad \forall v\in V. \end{aligned}$$

(3.8)

Assume $u=u_{0}+g_H$ where $u_{0}\in V$ and $g_H$ for simplicity is a linear combination of $\lambda _i, \,i=1,\dots , m$. The problem (3.8) is then equivalent to

$$\begin{aligned} \text {Find } u_{0}\in V: \quad v^TKu_{0} = - v^TKg_H, \quad \forall v\in V. \end{aligned}$$

(3.9)

The problem is in this way transformed back to the previous formulation and can be solved with the same method.

Proposition (3.4) states that the exact solution of the displacement problem is $u=u_\mathrm{ms} + u_f$. For the case with non-zero prescribed displacements and $g_H\in \text {span}(\{\lambda _i\}_{i=1}^m)$, which was described above, the correction vector $u_f$ turns out to be a linear combination of the vectors $\phi _i, \,i= 1,\dots , m$, as the following derivation shows. Consider the non-homogeneous displacement problem (3.9). The correction is attained from

$$\begin{aligned} \text {Find } u_f \in W: \quad w^TKu_f = -w^TKg_H, \quad \forall w\in W, \end{aligned}$$

(3.10)

where $g_H$ was assumed to be a linear combination of $\lambda _i$, i.e.:

$$\begin{aligned} g_H=\sum _{i=1}^m\alpha _i\lambda _i. \end{aligned}$$

(3.11)

The modification vectors $\phi _i, \,i=1,\dots , m$, is attained from the problems

$$\begin{aligned} \text {Find } \phi _i \in W: \quad w^TK\phi _i = w^TK\lambda _i, \quad \forall w\in W. \end{aligned}$$

(3.12)

Note that the problems are solved for all $i=1,\dots , m$, compared to before, when only $i\in \mathscr {M}$ were considered. This is necessary to construct the correction $u_f$ as will be seen next.

Inserting (3.11) into (3.10) and using (3.12) gives

$$\begin{aligned} w^TKu_f&= -w^TKg_H=-\sum _{i=1}^m \alpha _i w^TK\lambda _i=-\sum _{i=1}^m \alpha _i w^TK\phi _i=w^TK\left( -\sum _{i=1}^m \alpha _i \phi _i\right) , \end{aligned}$$

implying that the correction is the sum

$$\begin{aligned} u_f = -\sum _{i=1}^m\alpha _i\phi _i. \end{aligned}$$

For general fixed displacements $g_H$, not necessarily in $V_H$, see the work [6].

4 Error analysis

This paper concerns a quite general network model described by a connectivity matrix K and a right hand side load F. In this section, some error bounds are shown. Because of the generality of K, assumptions are needed. It is assumed that the coarse grid is a quasi-uniform finite element mesh with mesh parameter $H\approx m^{-1/2}$. The following two norms on the space V are introduced:

$$\begin{aligned} {\left| \left| \left| v \right| \right| \right| }&:=(v^TKv)^{1/2},\\ \Vert v\Vert _h&:=\left( \sum _i {\bar{h}}_i^2 v_i^2 \right) ^{1/2}, \end{aligned}$$

where ${\bar{h}}_i$ is the average length of all edges connected to the node corresponding to degree of freedom i. It is assumed that K is symmetric and positive definite on V, and that the smallest eigenvalue of K, denoted $\nu$, fulfils

$$\begin{aligned} \frac{{\left| \left| \left| v \right| \right| \right| }^2}{\Vert v\Vert _h^2}\ge \nu , \quad \forall v\in V, \end{aligned}$$

(4.1)

where $\nu$ is bounded from below by a constant independent of n. For the finite element method posed on a quasi-uniform mesh these definitions and assumptions correspond to the energy norm, a weighted $l^2$-norm and a Poincaré inequality. A bilinear weighted interpolant $\pi _H :V\rightarrow V_H$ is defined as $\pi _H v=\sum _{i\in \mathscr {M}}(\lambda _i^T v)\lambda _i$. i.e. $\pi _H v = B_HC_Hv$. Assume the following error bound:

$$\begin{aligned} \Vert v-\pi _H v\Vert _h\le C H{\left| \left| \left| v \right| \right| \right| }, \end{aligned}$$

(4.2)

which is expected to hold for $v\in V$.

Introducing the notation ${\bar{F}}$, with components ${\bar{F}}_i = F_i {\bar{h}}_i^{-2}$, the following variant of the Cauchy–Schwarz inequality can be derived using the original version,

$$\begin{aligned} |v^TF| = \left| \sum _i {\bar{h}}_i v_i {\bar{h}}_i^{-1} F_i\right| \le \left( \sum _i v_i^2 {\bar{h}}_i^2\right) ^{1/2} \left( \sum _i F_i^2 \bar{h}_i^{-2}\right) ^{1/2} = \Vert v\Vert _h \Vert {\bar{F}} \Vert _h. \end{aligned}$$

(4.3)

The weighted load vector ${\bar{F}}$ is O(1) in regard to n for distributed surface loads.

The main source of error in the proposed method is the localization of the multiscale basis functions to patches. To isolate this contribution the vector $\pi _H F \in V_H$ is introduced. Consider the modified model problem: find $\hat{u}\in V$ such that

$$\begin{aligned} v^T K \hat{u}=v^T \pi _H F,\,\,\,\forall v\in V. \end{aligned}$$

Given the assumptions, it is noted that

$$\begin{aligned} {\left| \left| \left| u-\hat{u} \right| \right| \right| }^2&=(u-{\hat{u}})^T (F-\pi _H F)\le \Vert u-{\hat{u}}\Vert _h \Vert {\bar{F}}-\pi _H {\bar{F}}\Vert _h\\&\le C{\left| \left| \left| u-\hat{u} \right| \right| \right| }\Vert {\bar{F}}-\pi _H {\bar{F}}\Vert _h, \end{aligned}$$

which together with (4.2) implies

$$\begin{aligned} {\left| \left| \left| u-\hat{u} \right| \right| \right| }\le C\Vert {\bar{F}}-\pi _H {\bar{F}}\Vert _h \le CH{\left| \left| \left| {\bar{F}} \right| \right| \right| }. \end{aligned}$$

(4.4)

This error is viewed as acceptable since the constant C is independent of data variation in the connectivity matrix, and without loss of generality it is assumed that $F\in V_H$.

First it is shown that the multiscale solution $u_\mathrm{ms}$ is equal to u if $F\in V_H$.

Proposition 4.1

Let u and $u_\mathrm{ms}$ be defined according to

$$\begin{aligned} u\in V: \, v^TKu&=v^TF, \,\,\,\forall v\in V,\\ u_\mathrm{ms}\in V_\mathrm{ms}: \, v^TKu_\mathrm{ms}&=v^TF, \,\,\,\forall v\in V_\mathrm{ms}. \end{aligned}$$

If $F\in V_H$, then it holds that $u=u_\mathrm{ms}$.

Proof

Since $F\in V_H$, there exists an $\bar{F}\in V$ such that $B_H C_H{\bar{F}}= F$. To show $v^TK u_\mathrm{ms}=v^T F, \, \forall v\in V$, it is noted that it holds for all $v\in V_{\text {ms}}$ from the definition of $u_\mathrm{ms}$. Due to the splitting $V=V_\mathrm{ms}\oplus W$ it is enough to show that it also is true for all test functions $v_f\in W$. For any $v_f\in W$ it holds that

$$\begin{aligned} v_f^T F-v_f^T K u_\mathrm{ms}=v_f^T B_H C_H {\bar{F}}=(\bar{F}^TC_H^T B_H^T v_f)^T=0, \end{aligned}$$

since $v_f^T K u_\mathrm{ms}=0$ by orthogonality and $B_H^T v_f=C_H v_f=0$ by the definition of the space W. Therefore $u_\mathrm{ms}$ solves the same equation as u, and due to uniqueness $u_\mathrm{ms}=u$. $\square$

The error committed by localization is difficult to study without stronger assumptions on the connectivity matrix K. It has however been analysed for several concrete cases, for instance when K is the stiffness matrix arising from discretizing the Poisson equation [15] and the elasticity equations [7] with the finite element method. For these cases it is true that

$$\begin{aligned} \max _{w\in V_\mathrm{ms}}\min _{v\in \tilde{V}_\mathrm{ms}}\frac{{\left| \left| \left| w-v \right| \right| \right| }}{{\left| \left| \left| w \right| \right| \right| }}\le Ce^{-c\rho }, \end{aligned}$$

(4.5)

where C and c are constants independent of H. In Sect. 6 it is shown, through numerical validation, that this relation is true for the more complicated network model considered there. The theoretical analysis of this result for the general network model proposed in Sect. 5 is complicated and postponed to future work. Assuming this result gives the following error bound.

Theorem 4.1

Assuming Eq. (4.5), the following error bound is true for any load vector $F\in V$:

$$\begin{aligned} {\left| \left| \left| u-{\tilde{u}}_\mathrm{ms} \right| \right| \right| }\le C\Vert {\bar{F}}-\pi _H \bar{F}\Vert _h+Ce^{-c\rho }\Vert {\bar{F}}\Vert _h. \end{aligned}$$

Proof

First consider the case $F\in V_H$. The vectors $u_\mathrm{}\in V_\mathrm{}$, $u_\mathrm{ms}\in V_\mathrm{ms}$ and ${\tilde{u}}_\mathrm{ms} \in {\tilde{V}}_\mathrm{ms}$ solve

$$\begin{aligned} \begin{aligned} v^TKu_\mathrm{}&=v^T F,\quad \forall v\in V_\mathrm{},\\ v^TKu_\mathrm{ms}&=v^T F,\quad \forall v\in V_\mathrm{ms},\\ v^TK{\tilde{u}}_\mathrm{ms}&=v^T F,\quad \forall v\in \tilde{V}_\mathrm{ms}. \end{aligned} \end{aligned}$$

(4.6)

Since ${\tilde{V}}_\mathrm{ms}\subset V$, Galerkin orthogonality gives

$$\begin{aligned} {\left| \left| \left| u-{\tilde{u}}_\mathrm{ms} \right| \right| \right| }\le \min _{v\in \tilde{V}_\mathrm{ms}} {\left| \left| \left| u-v \right| \right| \right| }. \end{aligned}$$

Moreover, since $F\in V_H$ it holds that $u=u_\mathrm{ms}$. Using this together with the assumption (4.5) with $w=u=u_\mathrm{ms}$ leads to

$$\begin{aligned} {\left| \left| \left| u_\mathrm{ms}-{\tilde{u}}_\mathrm{ms} \right| \right| \right| }\le \min _{v\in \tilde{V}_{\text {ms}}}{\left| \left| \left| u_\mathrm{ms}-v \right| \right| \right| }\le Ce^{-c\rho }{\left| \left| \left| u_\mathrm{ms} \right| \right| \right| }. \end{aligned}$$

This together with the fact that

$$\begin{aligned} {\left| \left| \left| u_\mathrm{ms} \right| \right| \right| }^2 = |u_\mathrm{ms}^TK u_\mathrm{ms}| = |u_\mathrm{ms}^TF|\le \Vert u_\mathrm{ms}\Vert \Vert {\bar{F}}\Vert _h \le C{\left| \left| \left| u_\mathrm{ms} \right| \right| \right| }\Vert {\bar{F}}\Vert _h, \end{aligned}$$

gives

$$\begin{aligned} {\left| \left| \left| u_\mathrm{ms}-\tilde{u}_\mathrm{ms} \right| \right| \right| }\le Ce^{-c\rho }\Vert {\bar{F}}\Vert _h. \end{aligned}$$

(4.7)

For the general case $F\in V$, let ${\hat{u}}, {\hat{u}}_\mathrm{ms}$ and $\hat{{\tilde{u}}}_\mathrm{ms}$ denote the different solutions to the problems in (4.6), but with load vectors $\pi _H F$. Using ${\hat{u}} = {\hat{u}}_\mathrm{ms}$ and the inequalities (4.4) and (4.7) finally result in

$$\begin{aligned} \begin{aligned} {\left| \left| \left| u-{\tilde{u}}_\mathrm{ms} \right| \right| \right| }&= {\left| \left| \left| u-{\hat{u}} + {\hat{u}}_\mathrm{ms} - \hat{{\tilde{u}}}_\mathrm{ms} + \hat{{\tilde{u}}}_\mathrm{ms} - {\tilde{u}}_\mathrm{ms} \right| \right| \right| } \\&\le {\left| \left| \left| u-{\hat{u}} \right| \right| \right| } + {\left| \left| \left| {\hat{u}}_\mathrm{ms} - \hat{{\tilde{u}}}_\mathrm{ms} \right| \right| \right| } + {\left| \left| \left| \hat{{\tilde{u}}}_\mathrm{ms} - {\tilde{u}}_\mathrm{ms} \right| \right| \right| }\\&\le C\Vert {\bar{F}}-\pi _H {\bar{F}}\Vert _h+Ce^{-c\rho }\Vert \pi _H {\bar{F}}\Vert _h + C\Vert {\bar{F}}-\pi _H {\bar{F}}\Vert _h \\&\le C\Vert {\bar{F}}-\pi _H {\bar{F}}\Vert _h+ Ce^{-c\rho }\left( \Vert {\bar{F}}-\pi _H {\bar{F}}\Vert _h + \Vert {\bar{F}}\Vert _h\right) \\&\le C\Vert {\bar{F}}-\pi _H {\bar{F}}\Vert _h+Ce^{-c\rho }\Vert {\bar{F}}\Vert _h. \end{aligned} \end{aligned}$$

$\square$

5 Network model for paper-based materials

In this section, a two-dimensional elasticity network model is presented, which can be used to model fiber networks in paper-based materials. An elasticity network consists of nodes and edges. When the nodes are displaced, internal forces act to restore the displacements. These forces act at two types of elements, either on edges or on edge pairs (two edges connected at a joint node). Three types of internal forces are included in this model, one type is related to edges, and two types are related to edge pairs. The model is two-dimensional, static and assumes small deformations.

The network mechanics are governed by force equilibrium equations assembled at each node. The general form of the equation for each node i reads

$$\begin{aligned} F^\text {Internal}_i + F^\text {External}_i = 0, \end{aligned}$$

where $F^\text {External}_i\in \mathbb {R}^2$ are all externally applied forces. The internal force $F^\text {Internal}_i\in \mathbb {R}^2$ is, as mentioned, a sum of three contributions:

$$\begin{aligned} F^\text {Internal}_i = F^\text {I}_i + F^\text {II}_i + F^\text {III}_i . \end{aligned}$$

The first force contribution is related to the edges of the network and acts to compensate for changes in length of the edges. The second force contribution acts on edge pairs to compensate changes in the angle between the edges of each edge pair. The third force contribution is included to model the Poisson effect. It acts on edge pairs by introducing a resistance to changes in the total length of the two edges of the pair. In the following sections, the three forces are described. Preparatory some nomenclature is introduced.

Let the network consist of N nodes. Let (i, j) denote the edge connecting node i to j and let $\mathscr {E}$ denote the set of all edges. Note that $(i, j) = (j, i)$. Edge pairs are denoted by (i, j, l) where j is the central node. Denote by $\mathscr {P}$ the set of all edge pairs. Note that $(i,j,l)=(l,j,i)$. Each node i has two degrees of freedom, the x-directed displacement and the y-directed displacement, contained in the vector $\delta _i\in \mathbb {R}^2$.

All force equilibrium equations can be assembled into a system of the form $-Ku + F = 0$ where $K\in \mathbb {R}^{n\times n}$ is called the elasticity matrix, $u\in \mathbb {R}^n$ is the node displacements, and $F\in \mathbb {R}^n$ is the external forces. The elasticity matrix K is attained by summation of matrices assembled at edges and edge pairs. The node displacement vector u is arranged according to

$$\begin{aligned} \begin{bmatrix} u(2i-1) \\ u(2i) \end{bmatrix} = \delta _i , \quad \quad 1 \le i \le N, \end{aligned}$$

and the elasticity matrix is assembled such that the force equilibrium equation for node i is at row $2i-1$ for the x-component, and at row 2i for the y-component.

Let the length of edge (i, j) be denoted $L_{ij}$ and assume that the edge has a width $w_{ij}$. All edges is assumed to have a uniform thickness z in the direction into the plane. The direction vector, $d_{ij}^a$, of an edge (i, j), with respect to node $a\in \{i,j\}$, is defined as

$$\begin{aligned} d_{ij}^a = \frac{p_b - p_a}{|p_b-p_a|},\quad b \in \{i,j\}, \,\, b\ne a. \end{aligned}$$

The length change of edge (i, j) is denoted $\varDelta L_{ij}$ and given by

$$\begin{aligned} \varDelta L_{ij} = \left( \delta _j-\delta _i\right) \cdot d_{ij}^i. \end{aligned}$$

In Fig. 7 some of the introduced notation is depicted. The length change $\varDelta L_{ij}$ of an edge, which is not exact but approximated by taking the dot product of the displacement difference onto the direction vector, is illustrated.

Fig. 7

Sketches showing the notation used for the network nodes, edges and edge pairs. To the right it is illustrated how the approximate length change $\varDelta L_{ij}$ of an edge (i, j) is calculated

5.1 Extension of edges

The first force contribution acts at edges due to their internal resistance to length change. When the nodes of an edge are displaced so that the projection of the difference of the node displacements onto the initial edge direction is nonzero, anti-parallel forces arise at the nodes of the edge to restore the length. The tendency of an edge to restore its length is described by the elastic modulus $k_{ij}$.

Consider an edge (i, j) as shown in Fig. 7b. When the nodes are displaced, the edge (i, j) will give rise to two forces $F_a^\text {I}(i,j)$, $a\in \{i,j\}$, acting on node a according to

$$\begin{aligned} F_a^\text {I}(i,j) = k_{ij} \frac{w_{ij}z}{L_{ij}}\varDelta L_{ij}d_{ij}^a, \quad a\in \{i,j\}. \end{aligned}$$

5.2 Angular deviations of edge pairs

The second force contribution acts at edge pairs from their internal tendency to resist change of the angle between their two edges. When a change in angle occurs, two torques arise at the connecting node acting on one edge each to restore the change. By transforming these torques to force couples the effect can be converted into the force equilibrium equations.

Consider an edge pair (i, j, l) as depicted in Fig. 7a. When the nodes are displaced, an angular change ${\varDelta \theta }_{ijl}$ occurs, giving rise to two torques

$$\begin{aligned} \begin{aligned} \tau _i&= \kappa _{ijl} V_{ijl} {\varDelta \theta }_{ijl} {\hat{z}},\\ \tau _l&= - \tau _i, \end{aligned} \end{aligned}$$

acting on edge (i, j) and (j, l) respectively, at the position of node j. The angular change is a sum of two contributions according to

$$\begin{aligned} {\varDelta \theta }_{ijl} = \delta \theta _{ji} + \delta \theta _{jl}. \end{aligned}$$

Each term, $\delta _{ji}$ and $\delta _{jl}$, is the angle deviation of respective edge from its initial orientation, as can be seen in Fig. 8. By using the assumption $\alpha \approx \tan \alpha$, the angles $\delta \theta _{ja}, a\in \{i,l\}$ can be calculated according to

$$\begin{aligned} \delta \theta _{ja} \approx \tan \theta _{ja} = \frac{(\delta _{a}-\delta _j)\cdot n^j_{ja}}{L_{ja}}, \quad a\in \{i,l\}, \end{aligned}$$

where the edge normals are calculated according to

$$\begin{aligned} \begin{aligned} n_{ji}^j&= d^j_{ji} \times {\hat{z}}, \\ n_{jl}^j&= -d^j_{jl} \times {\hat{z}}. \end{aligned} \end{aligned}$$

Transforming the torques to force couples gives the resulting three forces $F^\text {II}_a(i, j, l), a\in \{i,j,l\}$ from edge pair (i, j, l), acting on node a, according to

$$\begin{aligned} \begin{aligned} F^\text {II}_a(i,j, l)&= -\frac{\kappa _{ijl} V_{ijl} {\varDelta \theta }_{ijl}}{L_{aj}}n^j_{ja}, \quad a\in \{i,l\},\\ F_j^\text {II}(i,j,l)&= - F^\text {II}_i(i,j, l) - F^\text {II}_l (i,j, l). \end{aligned} \end{aligned}$$

These forces are illustrated in Fig. 8.

Fig. 8

The initial position of the edge pair (i, j, l) is shown with dashed lines. After displacement, forces act at the nodes of the edge pair to restore the angular change between the edges. The angle deviation of each edge is denoted $\delta _{ji}$ and $\delta _{jl}$, respectively

5.3 Poisson effect of edge pairs

The third force contribution results in an effect similar to the Poisson effect and acts at edge pairs. The idea is to add forces that work to keep the total length of the two edges of the pair constant. Hence, when one edge changes length, two kind of forces occur, on one hand forces acting to restore the length of the specific edge, on the other hand forces acting to change the length of the other edge in the pair.

Consider an edge pair (i, j, l), as shown in Fig. 7a. The forces acting at the outer nodes $a\in \{i,l\}$ will be

$$\begin{aligned} F^\text {III}_a(i,j,l)&= -\eta _{ijl} \frac{w_{aj}z}{L_{aj}}\left( \varDelta L_{aj}+ \gamma _{ijl}\frac{w_{bj}}{2}\frac{\varDelta L_{bj}}{L_{bj}}|n^j_{aj}\cdot d^j_{bj}|\right) d^j_{aj}, \quad \\&\quad a,b\in \{i,l\}, \,\, b\ne a, \end{aligned}$$

and at the central node

$$\begin{aligned} F_j^\text {III}(i,j,l) = - F^\text {III}_i(i,j,l) - F^\text {III}_l(i,j,l). \end{aligned}$$

5.4 Assembly of elasticity matrix

The governing equation $-{\bar{K}}u+F=0$ contains all node equilibrium equations as a matrix system. Here ${\bar{K}}$ is the modification of K by explicitly including boundary conditions. Since each edge and each edge pair leads to separate force contributions, the total elasticity matrix K can be assembled from separate element matrices for each of the three force contributions. Let $K^\text {I}_{ij}, K^\text {II}_{ijl}, K^\text {III}_{ijl}\in \mathbb {R}^{n\times n}$ denote the matrices assembled from the first, second and third force contribution respectively, at different elements [edges (i, j) or edge pairs (i, j, l)]. These matrices are sparse and the only nonzero elements are defined by the relations

$$\begin{aligned} \begin{aligned} K^\text {I}_{ij}(\{2a-1, 2a\}, \{1, \dots , n\}) u&= F_a^\text {I}(i,j), \quad \quad \, a\in \{i,j\}, \\ K^\text {II}_{ijl}(\{2a-1, 2a\}, \{1, \dots , n\}) u&= F_a^\text {II}(i,j, l), \quad \, a\in \{i,j, l\},\\ K^\text {III}_{ijl}(\{2a-1, 2a\}, \{1, \dots , n\}) u&= F_a^\text {III}(i,j, l), \quad a\in \{i,j, l\}. \end{aligned} \end{aligned}$$

The elasticity matrix K is assembled according to

$$\begin{aligned} K =&-\sum _{(i,j)\in \mathscr {E}} K^\text {I}_{ij} - \sum _{(i,j,l)\in \mathscr {P}}\left( K^\text {II}_{ijl} + K^\text {III}_{ijl}\right) . \end{aligned}$$

The matrix K is symmetric and semi-positive definite. With proper fixation of nodes, the matrix will be positive definite on the restricted solution space. Moreover, for a regular network with uniform coefficients k, $\kappa$, $\eta$ and $\gamma$, the presented model is equivalent to the finite difference discretization of the two-dimensional linear elasticity equations.

6 Numerical results

In this section, two network problems are solved using the proposed multiscale method and the fiber network model. Both problems are similar, with different boundary conditions and load vectors. Consider a unit square network with nodes and edges in a regular grid pattern, as shown in Fig. 9a. Let the number of nodes be $(r+1)^2$. In the examples $r=2^7=128$ will be used. Let $l_i$ and $\mu _i$ represent the standard Lamé parameters, with one value for each node i. Set the parameters of the network model to

$$\begin{aligned} \begin{aligned} k_{ij} = \frac{2{\bar{\mu }}_{ij} + {\bar{l}}_{ij}}{5c}, \quad \kappa _{ij} = \frac{{\bar{\mu }}_{ij} }{4c^2}, \quad \eta _{ijl} = \frac{2\mu _{j} + l_{j}}{5c},\quad \gamma _{ijl} = \frac{2 l_j}{4\eta _{ijl}c^2},\quad c = \frac{1}{2}, \end{aligned} \end{aligned}$$

where ${\bar{\mu }}_{ij} = \frac{\mu _i + \mu _j}{2}$ and ${\bar{l}}_{ij} = \frac{l_i + l_j}{2}$ are the mean values of the two nodes of edge (i, j). Using the above parameters, the elasticity matrix K is assembled. Let $h=1/r$ denote the length of each network edge.

Fig. 9

Example of a regular square network with $r=16$ and $R=4$

The first network problem, called the fixed boundary problem, is

$$\begin{aligned} {\bar{K}}_1 u = {\bar{F}}, \end{aligned}$$

(6.1)

where $u(i)=0$ for all network nodes on the unit square boundary, and $F(i) = \frac{1}{h^2}$. As mentioned earlier, ${\bar{K}}$ and ${\bar{F}}$ correspond to the modification of $K_1$ and F by explicitly including the boundary conditions.

The second problem, called the displaced boundary problem, is

$$\begin{aligned} {\bar{K}}_2 u = 0 \end{aligned}$$

(6.2)

where $u(i)=0$ for all network nodes with $x=0$, and $u(i) = 0.1$ for all x-directed degrees of freedom with $x=1$.

To solve the two problems using the proposed numerical multiscale method, a coarse FEM grid is introduced. The grid is similar to the network but with $(R+1)^2$ nodes. The basis functions $\varLambda _i$ are chosen as classic bilinear. See Fig. 9b for an illustration of the network and coarse FEM grid. Let $H= 1/R$ be the width of the coarse elements. With the described network geometry, the fixed boundary conditions correspond to fixation of coarse nodes at the boundary.

Each problem is solved with three different setups, first a basic setup with $l_i=\mu _i=1$. Secondly, by using a realization of random coefficients $l_i$ and $\mu _i$ sampled in [0.1, 10] with uniform distribution. Thirdly, each node i is displaced $[\delta x_i,\, \delta y_i]$ where $\delta x_i$ and $\delta y_i$ is randomly sampled in $[-0.4h, 0.4h]$ with uniform distribution. For nodes with initial coordinate $x=0$ or $x=1$, it is enforced that $\delta x_i = 0$, and similarly for nodes with $y=0$ or $y=1$, $\delta y_i=0$. In Fig. 10, a network with such random structure is shown for $r=32$.

Fig. 10

Example of a unit square network with $r=32$ and randomly perturbed nodes

The benefits of the multiscale method will be demonstrated by utilizing localization as described in Sect. 3.2, by introducing patches that the local modification problems are solved over. Which degrees of freedom that are included in each patch set $\mathscr {N}_E$, are chosen based on the radius $\rho H$. A degree of freedom $i\in \mathscr {N}$ will be in $\mathscr {N}_E$, if $|p_i-c_E|\le \rho H$, where $c_E\in \mathbb {R}^2$ is the center of element E. Hence patches will have the circular form as was illustrated in Fig. 6a. The rapid decay of the modified basis $\{\lambda _i-\phi _i\}_{i\in \mathscr {M}}$ is demonstrated by solving the modification ${\tilde{\phi }}_i$ with different $\rho$ and computing the relative error $\vert \vert \vert \phi _i-{\tilde{\phi }}_i\vert \vert \vert /\vert \vert \vert \phi _i\vert \vert \vert$. The degree of freedom i is chosen as one of the central nodes but the trend is similar for all nodes. The resulting errors for the first problem (6.1) are seen in Fig. 11.

Fig. 11

Relative error $\vert \vert \vert \phi _i-{\tilde{\phi }}_i\vert \vert \vert /\vert \vert \vert \phi _i\vert \vert \vert$ for a localized modified basis function ${\tilde{\phi }}_i$ for a central node for the first problem (6.1)

Moreover, the validity of the assumed Poincaré inequality (4.1) is investigated by solving the generalized eigenvalue problem

$$\begin{aligned} Kx = \nu Dx, \quad x\in V, \end{aligned}$$

(6.3)

and calculating the smallest eigenvalue. Here D denotes the diagonal matrix with values ${\bar{h}}_i^2, \,\, i=1, \dots , n$ on the diagonal. This is done for the three different setups and $r \in \{2^2, 2^3, 2^4, 2^5, 2^6\}$. In Fig. 12, it is seen that the resulting eigenvalues are independent of n. The values are normalized with the value for $r=2^6$.

Fig. 12

The smallest eigenvalue to the generalized eigenvalue problem (6.3) for the three different network setups and $r \in \{2^2, 2^3, 2^4, 2^5, 2^6\}$. For each case the smallest eigenvalue is normalized with the smallest eigenvalue for $r= 2^6$

Fig. 13

Relative errors for the first problem (6.1). The coarse mesh size H and its square $H^2$ are included in the plots to clarify convergence rates. For the basic setup and the setup with random coefficients, also the so-called FEM-error is included

Next, the two problems (6.1) and (6.2) are solved using $r= 2^7 = 128$ for different coarse grids with $R = 2, 4, 8, 16, 32$. The problems are solved using localization with patch radius $\rho H =C H \log _2(H^{-1})=k/2^k$ where $k=1,2,3,4,5$. For the first problem (6.1), the constant is chosen to $C=1$, and for the second problem (6.2) it is $C=1.5$. The second problem, with non-zero fixed displacement, is solved with correction as described in Sect. 3.4. The resulting relative errors for the two problem types and their three different setups are shown in Figs. 13 and 14. In some setups the errors can be compared with the so-called FEM-error, corresponding to solving the multiscale problem (3.6) with non-modified basis $B_H$ instead of the modified multiscale basis ${\tilde{B}}_\mathrm{ms}$. It can be seen that this FEM-solution behaves poorly for the setup with random coefficients.

Fig. 14

Realtive errors for the second problem (6.2) with correction as in Sect. 3.4. The coarse mesh size H and its square $H^2$ are included in the plots to clarify convergence rates

7 Conclusion and future work

In this paper a numerical multiscale method for discrete networks is proposed. For a set of different numerical examples, the convergence rates of the proposed method are examined. For regular networks with low connectivity variation the method resembles the convergence rates of the ordinary FEM. For networks with randomly varying connectivity it is shown that the multiscale method performs better than ordinary FEM. The method is moreover used to solve network problems with random structure, indicating error convergence rates at least linear in energy norm and quadratic in $l^2$-norm.

A challenging theoretical problem for the future is to extend the result in (4.5) beyond finite element based discretization to more general networks, including the presented network model. Further on, other fiber network models, for example beam based, as well as three-dimensional, should be considered. Thereby the presented numerical upscaling method will be applied to realistic macroscale paper networks, investigating the problem size which can be studied and the computational efficiency of the proposed method. Moreover, the proposed method will be used together with the paper forming simulation framework presented in [12, 14, 19] to study virtual paper sheets and macroscale mechanical properties such as tensile strength, tensile stiffness, bending resistance, z-strength and fracture propagation.