Multiscale characterization of the mechanics of curved fibered structures with application to biological and engineered materials

Computers and Structures 310 (2025) 107690 Available online 21 February 2025 0045-7949/© 2025 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). Contents lists available at ScienceDirect Computers and Structures journal homepage: www.elsevier.com/locate/compstruc Multiscale characterization of the mechanics of curved fibered structures with application to biological and engineered materials J.A. Sanz-Herrera ,∗, A. Apolinar-Fernandez, A. Jimenez-Aires, P. Perez-Alcantara, J. Dominguez, E. Reina-Romo Escuela Técnica Superior de Ingeniería, Universidad de Sevilla, Spain A R T I C L E I N F O A B S T R A C T Keywords: Multiscale formulation Finite element method Virtual tests Fibered microstructures Hydrogels Mechanics of biological tissues Curved fibered structures are ubiquitous in nature and the mechanical behavior of these materials is of pivotal importance in the biomechanics and mechanobiology fields. We develop a multiscale formulation to characterize the macroscopic mechanical nonlinear behavior from the microstructure of fibered matrices. From the analysis of the mechanics of a randomly curved single fiber, a fibered matrix model is built to determine the macroscopic behavior following a homogenization approach. The model is tested for tensile, compression and shear loads in different applications. The presented approach naturally recovers instabilities at compression as well as the strain stiffening regime, which are observed experimentally in the mechanical behavior of collagen matrices. Indeed, it was found that the bending energy associated to fiber unrolling, is the most important source of energy developed by fibers for the analyzed cases in tensile and shear in all deformation regions (except the strain stiffening region), whereas bending energy dominates at compression too during buckling. The proposed computational framework can also be used to perform multiscale simulations in engineered fibered materials. Therefore, the developed methodology may be an interesting and complementary tool to characterize the nonlinear behavior and evolution of curved fibered structures present in biology and engineering. 1. Introduction Fibrous networks play a fundamental role in many biological materials. Fibered tissues like the extracellular matrix (ECM) of biological tissues are composed by fibrils (e.g. collagen) connected to each other and surrounding nonfibrillar matrix. This hierarchical structure makes fibered materials macroscopic behavior to be largely dependent on its microstructural organization. Therefore, it is important to understand and characterize the behavior of fibrous tissues to provide insights into the pathophysiology of different diseases such as cancer or arteriosclerosis [1]. In many solid tumors there is a significant collagen fiber rearrangement and an increase of collagen fibers deposition and crosslinking [2,3]. These observations suggest that collagen fiber reorganization favors cancer progression. Understanding the mechanical behavior of fibered structures has also a huge biological interest to produce biomimicking materials within tissue engineering. Reproducing structures with similar properties to the ECM allows in vitro studies to resemble in vivo material behavior [4,5] that serves to organize and *Corresponding author at: Camino de los descubrimientos s/n, 41092 Seville, Spain. E-mail address: [email protected] (J.A. Sanz-Herrera). regulate cells and that may be relevant in developing new therapies for pathologies [6]. Due to the variety of fibered structures and their relevance in many biological processes, there is a widespread interest in the scientific community to understand and characterize the highly non linear mechanical behavior of these networks [7]. Based on high resolution fiber images available in the literature [4,8–10], fibrils can be found in different curved shapes such as twisted, crimped or non-regular curved fibrils. When fibrils present a curved shape, applying uniaxial deformations will unroll the fibril under a low stress state. Eventually, maintaining the strain rate will lead to a straighter shape of the fibril. As the distance between fibril ends becomes similar to the fibril length, fibril mechanical behavior shows rigidization [10,11]. Fibrils suffer higher stress in order to maintain strain rate, showing a linear behavior once fibril is completely straight [12]. This rigidization occurs not only in isolated fibrils. As fibril matrices are stretched, fibrils axial direction is aligned with the stretch direction, leading to the stiffening of the matrix mechanical response [13]. Therefore, this randomly-curved shape in the fiber https://doi.org/10.1016/j.compstruc.2025.107690 Received 17 December 2024; Accepted 15 February 2025 Computers and Structures 310 (2025) 107690 2 J.A. Sanz-Herrera, A. Apolinar-Fernandez, A. Jimenez-Aires et al. Fig. 1. Schematics of the fiber generator. (a) Definition of the discrete points 𝒙𝑖(for 𝑛=3). (b) Different individual randomly-curved fibers in 3D representation, for input parameters 𝑛= 3, 5 and 10; and 𝐿𝑓∕𝐿= 1.11, 1.21 and 1.4 for blue, green and magenta fibers, respectively. level (micro level) introduces a non-linear behavior in the stress-strain response at the tissue level (macro level) due to the stiffening associated to fibril unrolling and alignment in the loading direction. Modeling the mechanical behavior of fibered tissues has been the focus of the development of a wide variety of models (see [14] for a review). The mechanics of the fibril network is dependent of the spatial distribution of the fibers and the constitutive behavior of fibers and crosslinks [14]. An essential input for any fiber network model is the mechanical behavior of a single fiber. Its non-linear fibrous mechanical response can be modeled through either phenomenological or structural models. Phenomenological models are ad hoc models that fit fibril rigidization stress-strain curves considering exponential, polynomial or logarithmic functions [15–18]. In this modeling strategy, a strain energy density function is deployed to express the mechanical property of individual fibers and correlations with experimental results are needed to fit the model parameters [19–23]. Although phenomenological models predict the behavior of the fibers accurately, they do not account for any structural information. On the other hand, fibrils with different shape on their twisted state show different mechanical response in the early stretching phase before the completely straight fibril shape with the same physiological parameters. As fibril can present randomly-curved shapes, it is interesting to evaluate their structural evolution during the stretching process. Hence, it is possible to use structural models in order to describe explicitly fibril mechanical response, where physiological parameters will give the accurate stress-strain results for a specific fibril shape. This way, changes in fibril mechanical response due to only fibril shape differences are reproducible. Several structural models have been developed for fibrous structures by incorporating the network geometry in the model, either through image processing [24–30] or numerically generated [13,31–40]. Voronoi tessellations, Delaunay triangulations and other networks are used to generate the random fiber systems and use straight fibrils to replicate the fiber microstructure. Therefore, these models are limited as their non linear constitutive behavior comes from the gradual alignment of the fibers in the direction of the load applied but not from the physiological unrolling of the fibers. Others include waviness with collagen fibers in a load free state through helical springs [10] and predefined curved fibers [41–43]. In this work, a multiscale approach is used to mechanically characterize the non linear behavior of randomly curved fibered structures from the microstructure. The formulation is stated at finite strains and considers a linear elastic response of the fibers. Contrary to previous existing structural models, a fibered matrix model is built from the analysis of the mechanics of a randomly curved single fiber. The isotropic / anisotropic fiber distribution, their evolution during straining and the fluctuations of density are modeled at the micro level. Upon homogenizing these fibered matrices, the macroscopic behavior is captured at the macroscale. Therefore, the key aspect of this study is to recognize the importance of the curved shape of fibered structures on their three dimensional mechanical behavior under tensile, compression and shear loads. 2. Single fiber This section describes the mechanical analysis of a single fiber, from the generation of randomly curved single fibers, to the mathematical analysis of their mechanical behavior at tensile and compression tests. 2.1. Fiber generator Here we present the methodology to build different 3D randomly curved collagen fibers (see Fig. 1). We define the input parameter Las the distance between both ends of the fiber. This distance is discretized into nsegments, in which wavy portions of fibers are generated. An auxiliary polar system (𝑟, 𝜃)is defined along these segments (see Fig. 1a). The polar angle is randomly selected from a uniform distribution such that 𝜃∈[0,2𝜋). Therefore, the coordinates of the different discretized points along the fiber are given as: x𝑖=(𝑟cos 𝜃, 𝑟sin𝜃, 𝐿𝑖 𝑛 ), 𝑖=1...𝑛. (1) The parameter nrepresents in Eq. (1) the number of different portions that compose the fiber, along which the curved fiber is twisted throughout its length. On the other hand, the fiber eccentricity from the main axis, 𝑟in Eq. (1), is selected from a uniform distribution 𝑟∈[0,2𝑅]. 𝑅is defined as follows (see also Fig. 1a): 𝑅=√[(1+𝑃)𝐿 𝑛 ]2 −(𝐿 𝑛 )2 ,(2) with 𝑃(𝑃≥0) being a relative length parameter of the fiber length versus the initial distance between the end points. Once the 3-D spatial points 𝒙𝑖have been created (coarse discretization of the fiber), the curved fiber is smoothed by defining a spline through points 𝒙𝑖. This curve is then discretized (fine discretization of the fiber) into 𝑁𝐸 elements. The final length of the fiber 𝐿𝑓is then defined as ∑𝑁𝐸 𝑖=1 𝑙𝑖, being 𝑙𝑖the element length. Different generated individual fibers can be seen in Fig. 1b. 2.2. Mathematical formulation The mechanics of the 3D random fibers exposed above does not hold an analytical solution. Then, the generated curved fiber is discretized into 𝑁𝐸 Euler-Bernouilli beam finite elements. Each element contains 2 nodes and 6 degrees of freedom (3 translations and 3 rotations) per node. A linear elastic behavior of the fiber is assumed. We follow an updated lagrangian formulation such that the fiber geometry is updated at each current time 𝑡of analysis: x𝑖,𝑡+Δ𝑡=x𝑖,𝑡 +Δu𝑖(3) Computers and Structures 310 (2025) 107690 3 J.A. Sanz-Herrera, A. Apolinar-Fernandez, A. Jimenez-Aires et al. Fig. 2. Exact solution of a nonlinear cantilever rod subjected to a point dimensionless bending moment  𝑀=𝑀𝐿 𝐸𝐼 versus its finite element implementation. The analytical solution is an arc with dimensionless radius  𝑅=1∕  𝑀[45]. The rod was discretized into 150 elements, and the moment step was Δ 𝑀=𝜋∕500. with x𝑖,𝑡+Δ𝑡and x𝑖,𝑡 the position vectors of node 𝑖at configurations 𝑡+Δ𝑡 and 𝑡, respectively, and Δu𝑖the vector of displacements of node 𝑖from configuration 𝑡to 𝑡+Δ𝑡. This vector Δu𝑖is obtained after the discretization and assembly of global vector and matrix following a structural matrix finite element analysis [44]: 𝔸𝑁𝐸 𝑒=1 ΔF𝑒=𝔸𝑁𝐸 𝑒=1 {K𝑒(x𝑒 𝑖,𝑡)⋅Δu𝑒}(4) being 𝔸the assembly operator; and ΔF𝑒and Δu𝑒the incremental nodal vector of structural forces and displacements/rotations at element 𝑒, respectively. K𝑒(x𝑒 𝑖,𝑡)is the matrix of element 𝑒computed in the configuration 𝑡with element nodal coordinates x𝑒 𝑖,𝑡. After assembly, the (global) system in Eq. (4) yields, ΔF=K𝑡⋅Δu(5) The total nodal force vector is then updated as: F𝑡+Δ𝑡=F𝑡+ΔF(6) The finite element implementation of Eq. (5) was developed in Matlab®. The computer implementation was validated in Fig. 2versus the exact solution of a nonlinear cantilever rod subjected to a point bending moment. 2.3. Results In this section, we analyze the tensile/compression behavior of a single randomly curved fiber for a range of model parameters, namely, fiber end points length 𝐿, fiber diameter 𝑑, fiber elastic modulus 𝐸; and fiber length to fiber end points length ratio 𝐿𝑓∕𝐿(being this parameter an indication of the wavy/curly shape of the fiber). A baseline case is selected, with parameters 𝐿=20μm, 𝑑=0.2μm, 𝐸= 100 MPa and 𝐿𝑓∕𝐿=1.1. These parameters are taken as an estimation of the order of magnitude found in the literature for collagen (see Table 1). The displacements of the fiber are prescribed at its end (boundary) nodes, except in the point in which a uniaxial stretch is applied. On the other hand, the fiber can freely rotate at the end points. Fig. 3shows the mechanical behavior of a single fiber in uniaxial stretch-stress tests for the range of analyzed model parameters. The tensile region (𝜆>1) of the curves in Fig. 3shows a nonlinear behavior divided in two different regimes: first a nonlinear part followed by a strain stiffening due to fiber straightening. This behavior can be seen in Table 1 Range of parameters of natural collagen structures available in the literature. Dispersion of parameters might be due to different animal origin. Property Value References Linear modulus (MPa) 110-1470 [46] 0-32 [47] 2000-7000 [12] Fiber length (μm) 100-200 [12] 6-18 [46] 10-100 [10] Fiber diameter (μm) 0.1-0.5 [12] 0.418-0.446 [4] 0.01-0.25 [48] 0.205-0.905 [46] video S1 of the supplementary material, for the baseline fiber. Interestingly, it can be seen in the compression region (𝜆<1) in Fig. 3, that our formulation naturally recovers the instability (buckling) phenomenon developed by the fiber at compression. On the other hand, the incremental elastic energy of the fiber can be defined as follows: Δ𝑈= 𝑁𝐸 ∑ 𝑒=1 {F𝑒 𝑡 ⋅Δu𝑒}(7) The total accumulated elastic energy of the fiber is then obtained as: 𝑈𝑡+Δ𝑡=𝑈𝑡+Δ𝑈(8) The total accumulated elastic energy of the fiber is split from Eq. (8) into axial, bending and torsional energies as follows: •Axial energy is obtained by selecting the axial forces and axial displacements product components in Eq. (8). •Bending energy is obtained by selecting the shear forces/bending moments and shear displacements/bending rotations product components in Eq. (8). •Torsional energy is obtained by selecting the torque moment and torque rotation product components in Eq. (8). Fig. 4shows the different fiber accumulated energies for the uniaxial tensile/compression stretch-stress tests of the baseline single fiber. Fig. 4additionally plots the different fiber accumulated energies rela- Computers and Structures 310 (2025) 107690 4 J.A. Sanz-Herrera, A. Apolinar-Fernandez, A. Jimenez-Aires et al. Fig. 3. Parametric analysis of the mechanical behavior of elastic single randomly curved fiber. Model parameters are varied within the range found in the literature (see Table 1), from a baseline case that with parameters 𝐿=20μm, 𝑑=0.2μm, 𝐸= 100 MPa and 𝐿𝑓∕𝐿=1.1. Fig. 4. Total, axial, bending and torsional components of the accumulated energy (left) and relative energy to the total energy (right), for a single randomly curved fiber with parameters 𝐿=20μm, 𝑑=0.2μm, 𝐸= 100 MPa and 𝐿𝑓∕𝐿=1.1. Inset fibers in the figure on the right represent the deformation state of the fiber for the selected stretch points. tive to the total energy. According to Fig. 4, the first tensile region of the curve (represented by point 1) is governed mostly by fiber bending (with minor contribution of the axial and torsional energies), up to a certain stretch (point 2) in which the fiber is almost straight except a small portion of the fiber. Indeed, the deformation point 2 represents the onset of straightening of the fiber, according to the overall length of the fiber given by parameter 𝐿𝑓∕𝐿(equal to 1.1 in this case). Once the fiber is completely straight (point 3) the most important contribution of energy is due to axial deformation. In the compression regime (points 4 and 5) the most important contribution is the bending energy as the fiber is getting wavy (curly). The torsional energy is almost negligible in the different regimes of deformation of the analyzed fiber. It was checked (data not shown) that the total accumulated elastic internal energy of the fiber, obtained from Eqs. (7)-(8), is equal to the external energy developed by the applied force in the fiber. Note that the point 𝜆=1, at which all the energies are null, was not computed in Fig. 4to avoid this singularity. Keeping in mind that the first tensile/compression regions of the stretch-stress curve of the fiber behavior are governed by bending, and the strain stiffening region by the axial energy, the impact of model parameters shown in Fig. 3can be easily understood. Consistently, we can see a stiffer behavior (both at tensile and compression) for increasing fiber stiffness (𝐸) and diameter (𝑑); and decreasing length (𝐿). Indeed, in qualitative terms, the fiber behavior depends as ∝𝐸𝑑4∕𝐿3in the Computers and Structures 310 (2025) 107690 5 J.A. Sanz-Herrera, A. Apolinar-Fernandez, A. Jimenez-Aires et al. bending region; and as ∝𝐸𝑑2∕𝐿in the axial (strain stiffening) region. The effect of fiber wavy/curly (parameter 𝐿𝑓∕𝐿) is seen in Fig. 3as a stiffer behavior for decreasing 𝐿𝑓∕𝐿(less wavy fiber) with shorter bending region and, consequently, with faster development of the strain stiffening region. Indeed, the onset of the development of the strain stiffening region is correlated with parameter 𝐿𝑓∕𝐿. 3. Fibered matrix This section describes the multiscale mechanical analysis of a fibered matrix, composed of randomly curved fibers, from the mechanical interaction of fibers in a representative volume element (RVE). In particular, the impact of model parameters on the overall material behavior, including anisotropy, is analyzed at the matrix macroscopic level. It is assumed that the RVE is composed solely by fibers, neglecting the effect of surrounding aqueous substances in the matrix. 3.1. Matrix generator The generation of the RVE is established under the basis of the definition of RVE, that is a selected volume of the fiber microstructure that statistically represents the heterogeneity of the matrix [49]. Our RVE is composed of isolated fibers 𝑁𝑓𝑖and crosslink fibers 𝑁𝑓𝑥, such that the total number of fibers of the RVE is 𝑁𝑓 =𝑁𝑓𝑖+𝑁𝑓𝑥. Therefore, the volume concentration of fibers is defined as 𝑉𝑐=𝑁𝑓∕𝑉𝑅𝑉 𝐸 , with 𝑉𝑅𝑉 𝐸 being the volume of the RVE. The algorithm to generate fibered matrices proceeds as follows: Box 1: Algorithm to generate fibered matrices. 0. Set the volume concentration of fibers 𝑉𝑐and the volume of the RVE, and consequently, 𝑁𝑓, 𝑁𝑓𝑖and 𝑁𝑓𝑥. 1. Generate 𝑁𝑓𝑖isolated fibers (see section 2.1). 2. FOR 𝑚=1..𝑁𝑓𝑖 2.1 Randomly select 𝑁𝑚𝑥 crosslinking nodes of fiber 𝑚. 2.2 FOR 𝑗=1..𝑁𝑚𝑥 i. Create a set with potential connecting candidate crosslinking nodes 𝑘, for fibers 𝑛=1..𝑁𝑓𝑖(𝑛≠𝑚) such that 𝐿⋅(1 − 𝜖)≤𝐿𝑗−𝑘≤𝐿⋅(1 + 𝜖). ii. Randomly select node 𝑗from the set. iii. Generate a crosslinking fiber from nodes 𝑗−𝑘(with length 𝐿) as in section 2.1. END FOR. END FOR. We set 𝑁𝑚𝑥 =1 in Box 1, meaning that a crosslinking fiber is generated per isolated fiber. Therefore, 𝑁𝑚𝑥 represents the degree of crosslinking of the matrix. On the other hand, the endpoints length of the crosslinking fiber is set to the same length 𝐿of the isolated fiber, according to item iin Box 1, up to a tolerance 𝜖which is set to 0.01 in our code. 3.2. Multiscale formulation The overall macroscopic mechanical behavior of a fibered matrix and its evolution, is obtained from the micromechanical interaction of fibers within the RVE in a multiscale fashion (see Fig. 5). In the macroscale, we define the (Cauchy) stress tensor and the (logarithmic) strain tensor as 𝝈𝑀 𝐼,𝑡(𝐗𝐼,𝑡),E𝑀 𝐼,𝑡(𝐗𝐼,𝑡), respectively, in a macroscopic (Gauss) material point 𝐼of the macroscale (with coordinates 𝐗𝐼,𝑡), for load time 𝑡(see Fig. 5). The (true) Cauchy’s stress components are used in the formulation for convenience and subsequent validation with experimental outcomes. On the other hand, in the microscale, we define the variables F𝑚 𝑡(𝐱𝑖,𝑡),u𝑚 𝑡(𝐱𝑖,𝑡)associated to an Euler-Bernoulli beam in a material point (node) 𝑖of the microscale (with coordinates 𝐱𝑖,𝑡) for load time Fig. 5. Schematics of the proposed multiscale formulation of fibered matrices. 𝑡. F𝑚 𝑡(𝐱𝑖,𝑡)represents the vector of nodal structural forces (forces and moments), whilst u𝑚 𝑡(𝐱𝑖,𝑡)is the vector that contains nodal displacements and rotations (see Fig. 5). As stated in the analysis of the single fiber, we follow an updated lagrangian approach such that the configuration (finite element mesh) is updated at each load time 𝑡. Consequently, the macroscopic behavior of the matrix is geometrically nonlinear and established at finite strains through the logarithmic strain tensor. The logarithmic strain tensor is obtained by integrating the strain rate, numerically, in a material frame of reference [50]: E𝑀 𝐼,𝑡+Δ𝑡=ΔR⋅E𝑀 𝐼,𝑡 ⋅ΔR𝑇+ΔE𝑀 𝐼(9) where E𝑀 𝐼,𝑡+Δ𝑡and E𝑀 𝐼,𝑡 are the total strains at increments 𝑡+Δ𝑡and 𝑡, respectively; ΔRis the incremental rotation tensor; and ΔE𝑀 𝐼is the total strain increment from increment 𝑡to 𝑡+Δ𝑡. Furthermore, the macroscopic strain energy density from increment 𝑡to 𝑡+Δ𝑡is defined as: Δ𝑈𝑀 𝐼=𝝈𝑀 𝑡+Δ𝑡,𝐼 ⋅ΔE𝑀 𝐼(10) Analogously, the kinematics at the microscale is defined as: u𝑚 𝑡+Δ𝑡=u𝑚 𝑡+Δu𝑚(11) where u𝑚 𝑡+Δ𝑡and u𝑚 𝑡are the total displacements/rotations vectors at increments 𝑡+Δ𝑡and 𝑡, respectively; and Δu𝑚is the total displacement increment vector from increment 𝑡to 𝑡+Δ𝑡. The displacements increment vector is defined at nodes 𝑘of the microscale through the strain increment at the macroscale as: Δu𝑚 𝑘=ΔE𝑀 𝐼 ⋅x𝑘,𝑡 (12) where x𝑘,𝑡 are the (updated) microscopic coordinates, at load time 𝑡, of the points 𝑘in which the macroscopic strain is affinely transmitted to the microscale. We assume that points 𝑘are the endpoints (knots) of the fibers. After discretization and FE assembly, the global nodal forces vector is updated as follows, F𝑚 𝑡+Δ𝑡=F𝑚 𝑡+ΔF𝑚(13) Using Eq. (5) in (13) yields: F𝑚 𝑡+Δ𝑡=F𝑚 𝑡+𝐊𝑡⋅Δu𝑚(14) Finally, the total microscopic strain energy density from increment 𝑡to 𝑡+Δ𝑡is defined as: Δ𝑈𝑚=1 𝑉𝑅𝑉 𝐸 𝑁 ∑ 𝑖=1 𝐹𝑚 𝑖,𝑡+Δ𝑡 ⋅Δ𝑢𝑚 𝑖(15) Computers and Structures 310 (2025) 107690 6 J.A. Sanz-Herrera, A. Apolinar-Fernandez, A. Jimenez-Aires et al. where 𝑁is the number of nodes of the microstructure. 𝐹𝑚 𝑖,𝑡+Δ𝑡and Δ𝑢𝑚 𝑖 represent the 𝑖components of vectors F𝑚 𝑡+Δ𝑡and Δu𝑚, respectively. On the other hand, we use the Hill’s equality of energies from the transition from the macro to the microscale [49,51] Δ𝑈𝑀 𝐼=Δ𝑈𝑚(16) Thus, using (10) and (15) in (16) yields, 𝝈𝑀 𝑡+Δ𝑡,𝐼 ⋅ΔE𝑀 𝐼=1 𝑉𝑅𝑉 𝐸 𝑁 ∑ 𝑖=1 𝐹𝑚 𝑖,𝑡+Δ𝑡 ⋅Δ𝑢𝑚 𝑖(ΔE𝑀 𝐼)(17) Note that in Eq. (15), Δ𝑢𝑚 𝑖(ΔE𝑀 𝐼)represents the incremental nodal displacements/rotations of the microstructure, as a result of the solution of the discrete finite element system ΔF𝑚=K𝑡⋅Δu𝑚with prescribed displacements at discretized nodes of knots 𝑘according to Eq. (12). Therefore, the different 𝑘𝑙 components of the stress tensor 𝜎𝑀 𝑡+Δ𝑡,𝐼,𝑘𝑙 at macroscopic point 𝐼(see Fig. 5) are obtained using Eq. (17) as follows: 𝜎𝑀 𝑡+Δ𝑡,𝐼,𝑘𝑙 =1 𝑉𝑅𝑉 𝐸 𝑁 ∑ 𝑖=1 𝐹𝑚 𝑖,𝑡+Δ𝑡 ⋅Δ𝑢𝑚 𝑖(Δ𝐸𝑀 𝐼,𝑘𝑙)⋅ 1 Δ𝐸𝑀 𝐼,𝑘𝑙 (18) Eq. (18) is equivalent to, 𝝈𝑀 𝑡+Δ𝑡,𝐼 =1 𝑉𝑅𝑉 𝐸 𝑁 ∑ 𝑖=1 𝐹𝑚 𝑖,𝑡+Δ𝑡 ⋅Δ𝑢𝑚 𝑖(𝕀)(19) with 𝕀being the unit tensor: 𝕀=1 2(𝛿𝑖𝑘𝛿𝑗𝑙 +𝛿𝑖𝑙 +𝛿𝑗𝑘)(20) Δ𝑢𝑚 𝑖(𝕀)represents the incremental nodal displacements/rotations of the microstructure, as a result of the solution of the discrete finite element system ΔF𝑚=K𝑡⋅Δu𝑚with prescribed displacements at discretized nodes of knots 𝑘for unit strain tensors 𝕀, following Eq. (12). Using the decomposition of the nodal forces vector (13) in (18), the incremental stress tensor can be obtained as: Δ𝝈𝑀 𝐼=1 𝑉𝑅𝑉 𝐸 𝑁 ∑ 𝑖=1 (𝐾𝑡,𝑖𝑗 ⋅Δ𝑢𝑚 𝑗(ΔE𝑀 𝐼))⋅Δ𝑢𝑚 𝑖(𝕀)(21) The macroscopic tangent stiffness tensor can be then obtained as: ℂ𝑡=𝜕Δ𝝈𝑀 𝐼 𝜕ΔE𝑀 𝐼 =1 𝑉𝑅𝑉 𝐸 𝑁 ∑ 𝑖=1 (𝐾𝑡,𝑖𝑗 ⋅Δ𝑢𝑚 𝑗(𝕀))⋅Δ𝑢𝑚 𝑖(𝕀)(22) Box 2 summarizes the main steps of the developed multiscale algorithm. Box 2: Algorithm to compute macroscopic stress and strain quantities from microscopic fibered matrices. 1. Macroscale: E𝑀 𝐼,𝑡, 𝝈𝑀 𝐼,𝑡 available at load increment 𝑡. 2. Macroscale: Set next strain increment ΔE𝑀 𝐼. 3. Microscale: Obtain Δ𝑢𝑚 𝑗(ΔE𝑀 𝐼)and Δ𝑢𝑚 𝑖(𝕀)from result of the solution of the system ΔF𝑚=K𝑡⋅Δu𝑚with prescribed displacements at nodes 𝑘for strain tensor ΔE𝑀 𝐼and unit strain tensors 𝕀, respectively, following Eq. (12). 4. Microscale: Update microscopic coordinates as x𝑖,𝑡+Δ𝑡=x𝑖,𝑡 + Δu𝑚 𝑖(ΔE𝑀 𝐼). 5. Macroscale: Compute Δ𝝈𝑀 𝐼according to Eq. (21). Compute the macroscopic tangent stiffness tensor ℂaccording to Eq. (22). 6. Macroscale: Update macroscopic quantities as 𝝈𝑀 𝐼,𝑡+Δ𝑡= 𝝈𝑀 𝐼,𝑡 +Δ𝝈𝑀 𝐼, E𝑀 𝐼,𝑡+Δ𝑡=E𝑀 𝐼,𝑡 +ΔE𝑀 𝐼 7. 𝑡←𝑡+Δ𝑡. GO TO 1. Finally, the incremental elastic internal microscopic energy of a fiber 𝑓of the microstructure of the matrix, Δ𝑈𝑚,𝑓 , is computed as, Δ𝑈𝑚,𝑓 = 𝑁𝐸𝑓 ∑ 𝑒=1 {F𝑚,𝑒 𝑡 ⋅Δu𝑚,𝑒}(23) where F𝑚,𝑒 𝑡and Δu𝑚,𝑒 are the total structural forces and displacements/rotations nodal vectors at element 𝑒, respectively. 𝑁𝐸𝑓is the number of elements of the fiber 𝑓. The total accumulated elastic energy of each fiber 𝑓of the microstructure of the matrix is obtained as: 𝑈𝑚,𝑓 𝑡+Δ𝑡=𝑈𝑚,𝑓 𝑡+Δ𝑈𝑚,𝑓 .(24) On the other hand, the total accumulated microscopic strain energy density is computed as, 𝑈𝑚 𝑡+Δ𝑡=𝑈𝑚 𝑡+Δ𝑈𝑚(25) with Δ𝑈𝑚obtained from Eq. (15). 3.3. Results: tensile/compression In this section, different fibered microstructures—with different microstructural characteristics—are subjected to a macroscopic tensile and compression stress state. The stress-stretch curves of the virtual tests are then obtained following the multiscale approach exposed in section 3.2. For this purpose, a vertical stretch velocity  𝜆(𝑌-direction) is prescribed. The macroscopic incremental strain tensor (point 2 of Box 2 of the multiscale formulation) is obtained from the solution of: Δ𝝈𝑀=ℂ𝑡⋅ΔE𝑀(26) ΔE𝑀is obtained from the solution of the system in Eq. (26) with uniaxial (mixed boundary conditions) stress state. In this sense, Δ𝐸𝑀 𝑌= Δ𝑡⋅ 𝜆for the macroscopic incremental strain tensor, and Δ𝜎𝑀 𝑋=Δ𝜎𝑀 𝑍= Δ𝜎𝑀 𝑋𝑌 =Δ𝜎𝑀 𝑋𝑍 =Δ𝜎𝑀 𝑌𝑍 =0. Fig. 6shows the mechanical behavior of the considered matrices in uniaxial stretch-stress tests. In this figure, different randomly curved fibered matrices, at different fiber versus volume concentration (parameter 𝑉𝑐), were generated for the range of analyzed model parameters. These parameters were varied within the order of magnitude of real biological fibered matrices (see Table 1), from a baseline case with parameters 𝐿=30μm, 𝑑=0.4μm, 𝐸= 150 MPa, 𝐿𝑓∕𝐿=1.15 and 𝑉𝑐=2.05 ⋅106fibers∕𝑚𝑚3. The microstructures were computed using the matrix generator considering that fibers can be distributed and oriented freely along the RVE without any preferential direction, such that the microstructures are assumed isotropic. Similarly to the single fiber, the tensile region (𝜆>1) of the curves in Fig. 6shows a linear behavior at small strains followed by a nonlinear behavior. In this case, the nonlinear regime is the consequence of both fiber straightening and fiber alignment with the direction of load. Besides, the loss of stiffness due to the instability (buckling) phenomenon developed by the fibers at compression, can be seen in the compression region (𝜆<1) in Fig. 6. According to Fig. 6, and analogously to the single fiber behavior, the overall tensile/compression behavior of the matrices is stiffer for increasing fiber stiffness (𝐸), fiber diameter (𝑑) and concentration of fibers. On the other hand, the overall tensile/compression behavior of the matrices is softer for increasing fiber length (𝐿) and fiber curvature (𝐿∕𝐿𝑓). The total accumulated strain energy density developed by the microstructure, computed from Eq. (25), can be seen in Fig. 7for the uniaxial tensile/compression stretch-stress test. Results are therefore shown as the (accumulated) microscopic contribution of all the fibers of the matrix along the deformation path, for the considered baseline matrix defined above. Note that this quantity is equal to the total accumulated strain energy density developed at the macrostructure (Eq. (10)), by virtue of Eq. (16). In Fig. 7, the axial, bending and torsional components of the total accumulated strain energy density are also shown. Computers and Structures 310 (2025) 107690 7 J.A. Sanz-Herrera, A. Apolinar-Fernandez, A. Jimenez-Aires et al. Fig. 6. Tensile/compression tests. Parametric analysis of the mechanical behavior of elastic randomly curved fibered matrices. Model parameters are varied within the range found in the literature (see Table 1), from a baseline case with parameters 𝐿=30μm, 𝑑=0.4μm, 𝐸= 150 MPa, 𝐿𝑓∕𝐿=1.15 and 𝑉𝑐=2.05 ⋅106fibers∕𝑚𝑚3. These quantities were obtained from Eq. (25), after selecting the axial, bending and torsional components similarly to the single fiber case (section 2.3). As in the single fiber case, it can be seen in this figure the asymmetry of the total strain energy density at tensile (𝜆>1) versus the compression regime (𝜆<1). Indeed, matrix compression induces fiber instabilities (buckling) whereas matrix tensile produces fiber unrolling and fiber alignment with load (strain stiffening). Fig. 7additionally plots the different accumulated strain energy densities relative to the total strain energy density. On the other hand, Fig. 8shows the accumulated multiscale energies in the fibers of the microstructure of the matrix for different tensile/compression deformation points, as defined in Fig. 7. These quantities were computed from Eq. (24) after selecting the axial, bending and torsional components. The tensile/compression behavior of the baseline matrix, the evolution of its microstructure, and both the axial, bending and torsional accumulated energies in the fibers; can be seen in videos S2, S3 and S4 of the supplementary material. It can be seen in Fig. 7that bending energy dominates in all the analyzed deformation regime, except the strain stiffening region. At the tensile region (𝜆>1) this energy is employed in fiber unrolling and fiber alignment with load, as corroborated in the microstructure evolution by Fig. 8and video S3. At the compression region (𝜆<1) the bending energy is employed in fiber buckling, as corroborated in the microstructure evolution (see also Fig. 8and video S3). Moreover, the torsional energy is almost negligible in all the analyzed deformation range, according to Figs. 7and 8. 3.4. Results: simple shear Analogously to the previous section, different fibered microstructures (assumed isotropic) are subjected to a macroscopic simple shear stress state. In this case, a shear strain velocity (𝛾) is prescribed, such that the macroscopic incremental strain tensor (point 2 of Box 2 of the Computers and Structures 310 (2025) 107690 8 J.A. Sanz-Herrera, A. Apolinar-Fernandez, A. Jimenez-Aires et al. Fig. 7. Multiscale energies: Macroscopic. Total, axial, bending and torsional accumulated strain energy densities (left), and relative to the total strain energy density (right), developed by the microstructure for tensile/compression tests of the baseline matrix. Fig. 8. Multiscale energies: Microscopic (fiber scale). Plot in the microstructure of the axial (first row), bending (second row) and torsional (third row) of the accumulated energies in the fibers of the baseline matrix for tensile/compression tests. Results are shown at different deformation points (1,2,3) defined in Fig. 7. Computers and Structures 310 (2025) 107690 9 J.A. Sanz-Herrera, A. Apolinar-Fernandez, A. Jimenez-Aires et al. Fig. 9. Simple shear tests. Parametric analysis of the mechanical behavior of elastic randomly curved fibered matrices. Model parameters are varied within the range found in the literature (see Table 1), from a baseline case with parameters 𝐿=30μm, 𝑑=0.4μm, 𝐸= 150 MPa, 𝐿𝑓∕𝐿=1.15 and 𝑉𝑐=2.05 ⋅106fibers∕𝑚𝑚3. multiscale formulation) is fully prescribed as a null value for all components except for Δ𝐸𝑀 𝑋𝑌 =Δ𝑡⋅𝛾∕2. Fig. 9shows the mechanical behavior of the considered matrices in simple shear stress tests. In this figure, we assumed the same range of variation of the parameters, as well as the same baseline case, than the tensile/compression tests. Therefore, the same microstructures (RVEs) defined for the tensile/compression tests in the previous section, are used here to simulate simple shear tests. The shear behavior of the (virtually) tested matrices, Fig. 9, shows also a nonlinear profile. First, a linear region at small strains followed by a nonlinear region, also as a consequence of both fiber straightening (unrolling) and fiber alignment with the direction of load. Analogously to the tensile tests, the shear behavior of the matrices is stiffer for increasing fiber stiffness (𝐸), fiber diameter (𝑑) and concentration of fibers. On the other hand, the overall tensile/compression behavior of the matrices is softer for increasing fiber length (𝐿) and fiber curvature (𝐿∕𝐿𝑓). As in the previous section, the accumulated strain energy densities developed by the microstructure for the shear test can be seen in Fig. 10. In this figure, results are shown as the (accumulated) microscopic contribution of all the fibers of the matrix along the deformation path, for the considered baseline matrix. The axial, bending and torsional components of the strain energy densities were computed analogously to the previous section. Fig. 10 additionally plots the different accumulated strain energy densities relative to the total strain energy density. On the other hand, Fig. 11 shows the accumulated multiscale energies in the fibers of the microstructure of the matrix for different shear deformation points, as defined in Fig. 10. The shear behavior of the baseline matrix, the evolution of its microstructure, and both the axial, bending and torsional accumulated energies in the fibers; can be seen