Compu e s and S uc u es 310 (2025) 107690
A ailable online 21 Feb ua y 2025
0045-7949/© 2025 The Au ho s. Published by Else ie L d. This is an open access a icle unde he CC BY-NC-ND license (h p://c ea i ecommons.o g/licenses/by-
nc-nd/4.0/).
Con en s lis s a ailable a ScienceDi ec
Compu e s and S uc u es
jou nal homepage: www.else ie .com/loca e/comps uc
Mul iscale cha ac e iza ion o he mechanics o cu ed fibe ed s uc u es
wi h applica ion o biological and enginee ed ma e ials
J.A. Sanz-He e a ,∗, A. Apolina -Fe nandez, A. Jimenez-Ai es, P. Pe ez-Alcan a a,
J. Dominguez, E. Reina-Romo
Escuela Técnica Supe io de Ingenie ía, Uni e sidad de Se illa, Spain
A R T I C L E I N F O A B S T R A C T
Keywo ds:
Mul iscale o mula ion
Fini e elemen me hod
Vi ual es s
Fibe ed mic os uc u es
Hyd ogels
Mechanics o biological issues
Cu ed fibe ed s uc u es a e ubiqui ous in na u e and he mechanical beha io o hese ma e ials is o pi o al
impo ance in he biomechanics and mechanobiology fields. We de elop a mul iscale o mula ion o cha ac e ize
he mac oscopic mechanical nonlinea beha io om he mic os uc u e o fibe ed ma ices. F om he analysis
o he mechanics o a andomly cu ed single fibe , a fibe ed ma ix model is buil o de e mine he mac oscopic
beha io ollowing a homogeniza ion app oach. The model is es ed o ensile, comp ession and shea loads in
diffe en applica ions. The p esen ed app oach na u ally eco e s ins abili ies a comp ession as well as he s ain
s iffening egime, which a e obse ed expe imen ally in he mechanical beha io o collagen ma ices. Indeed, i
was ound ha he bending ene gy associa ed o fibe un olling, is he mos impo an sou ce o ene gy de eloped
by fibe s o he analyzed cases in ensile and shea in all de o ma ion egions (excep he s ain s iffening egion),
whe eas bending ene gy domina es a comp ession oo du ing buckling. The p oposed compu a ional amewo k
can also be used o pe o m mul iscale simula ions in enginee ed fibe ed ma e ials. The e o e, he de eloped
me hodology may be an in e es ing and complemen a y ool o cha ac e ize he nonlinea beha io and e olu ion
o cu ed fibe ed s uc u es p esen in biology and enginee ing.
1. In oduc ion
Fib ous ne wo ks play a undamen al ole in many biological ma-
e ials. Fibe ed issues like he ex acellula ma ix (ECM) o biologi-
cal issues a e composed by fib ils (e.g. collagen) connec ed o each
o he and su ounding nonfib illa ma ix. This hie a chical s uc u e
makes fibe ed ma e ials mac oscopic beha io o be la gely dependen
on i s mic os uc u al o ganiza ion. The e o e, i is impo an o un-
de s and and cha ac e ize he beha io o fib ous issues o p o ide
insigh s in o he pa hophysiology o diffe en diseases such as cance
o a e ioscle osis [1]. In many solid umo s he e is a significan col-
lagen fibe ea angemen and an inc ease o collagen fibe s deposi ion
and c osslinking [2,3]. These obse a ions sugges ha collagen fibe
eo ganiza ion a o s cance p og ession. Unde s anding he mechani-
cal beha io o fibe ed s uc u es has also a huge biological in e es o
p oduce biomimicking ma e ials wi hin issue enginee ing. Rep oduc-
ing s uc u es wi h simila p ope ies o he ECM allows in i o s udies
o esemble in i o ma e ial beha io [4,5] ha se es o o ganize and
*Co esponding au ho a : Camino de los descub imien os s/n, 41092 Se ille, Spain.
E-mail add ess: [email p o ec ed] (J.A. Sanz-He e a).
egula e cells and ha may be ele an in de eloping new he apies o
pa hologies [6].
Due o he a ie y o fibe ed s uc u es and hei ele ance in many
biological p ocesses, he e is a widesp ead in e es in he scien ific com-
muni y o unde s and and cha ac e ize he highly non linea mechanical
beha io o hese ne wo ks [7]. Based on high esolu ion fibe images
a ailable in he li e a u e [4,8–10], fib ils can be ound in diffe en
cu ed shapes such as wis ed, c imped o non- egula cu ed fib ils.
When fib ils p esen a cu ed shape, applying uniaxial de o ma ions
will un oll he fib il unde a low s ess s a e. E en ually, main aining
he s ain a e will lead o a s aigh e shape o he fib il. As he dis ance
be ween fib il ends becomes simila o he fib il leng h, fib il mechan-
ical beha io shows igidiza ion [10,11]. Fib ils suffe highe s ess in
o de o main ain s ain a e, showing a linea beha io once fib il is
comple ely s aigh [12]. This igidiza ion occu s no only in isola ed
fib ils. As fib il ma ices a e s e ched, fib ils axial di ec ion is aligned
wi h he s e ch di ec ion, leading o he s iffening o he ma ix mechan-
ical esponse [13]. The e o e, his andomly-cu ed shape in he fibe
h ps://doi.o g/10.1016/j.comps uc.2025.107690
Recei ed 17 Decembe 2024; Accep ed 15 Feb ua y 2025
Compu e s and S uc u es 310 (2025) 107690
2
J.A. Sanz-He e a, A. Apolina -Fe nandez, A. Jimenez-Ai es e al.
Fig. 1. Schema ics o he fibe gene a o . (a) Defini ion o he disc e e poin s 𝒙𝑖( o 𝑛=3). (b) Diffe en indi idual andomly-cu ed fibe s in 3D ep esen a ion, o
inpu pa ame e s 𝑛= 3, 5 and 10; and 𝐿𝑓∕𝐿= 1.11, 1.21 and 1.4 o blue, g een and magen a fibe s, espec i ely.
le el (mic o le el) in oduces a non-linea beha io in he s ess-s ain
esponse a he issue le el (mac o le el) due o he s iffening associa ed
o fib il un olling and alignmen in he loading di ec ion.
Modeling he mechanical beha io o fibe ed issues has been he
ocus o he de elopmen o a wide a ie y o models (see [14] o a
e iew). The mechanics o he fib il ne wo k is dependen o he spa-
ial dis ibu ion o he fibe s and he cons i u i e beha io o fibe s and
c osslinks [14]. An essen ial inpu o any fibe ne wo k model is he
mechanical beha io o a single fibe . I s non-linea fib ous mechanical
esponse can be modeled h ough ei he phenomenological o s uc-
u al models. Phenomenological models a e ad hoc models ha fi fib il
igidiza ion s ess-s ain cu es conside ing exponen ial, polynomial o
loga i hmic unc ions [15–18]. In his modeling s a egy, a s ain ene gy
densi y unc ion is deployed o exp ess he mechanical p ope y o in-
di idual fibe s and co ela ions wi h expe imen al esul s a e needed o
fi he model pa ame e s [19–23]. Al hough phenomenological models
p edic he beha io o he fibe s accu a ely, hey do no accoun o any
s uc u al in o ma ion. On he o he hand, fib ils wi h diffe en shape
on hei wis ed s a e show diffe en mechanical esponse in he ea ly
s e ching phase be o e he comple ely s aigh fib il shape wi h he
same physiological pa ame e s. As fib il can p esen andomly-cu ed
shapes, i is in e es ing o e alua e hei s uc u al e olu ion du ing he
s e ching p ocess. Hence, i is possible o use s uc u al models in o de
o desc ibe explici ly fib il mechanical esponse, whe e physiological pa-
ame e s will gi e he accu a e s ess-s ain esul s o a specific fib il
shape. This way, changes in fib il mechanical esponse due o only fib il
shape diffe ences a e ep oducible. Se e al s uc u al models ha e been
de eloped o fib ous s uc u es by inco po a ing he ne wo k geome y
in he model, ei he h ough image p ocessing [24–30] o nume ically
gene a ed [13,31–40]. Vo onoi essella ions, Delaunay iangula ions
and o he ne wo ks a e used o gene a e he andom fibe sys ems and
use s aigh fib ils o eplica e he fibe mic os uc u e. The e o e, hese
models a e limi ed as hei non linea cons i u i e beha io comes om
he g adual alignmen o he fibe s in he di ec ion o he load applied
bu no om he physiological un olling o he fibe s. O he s include
wa iness wi h collagen fibe s in a load ee s a e h ough helical sp ings
[10] and p edefined cu ed fibe s [41–43].
In his wo k, a mul iscale app oach is used o mechanically cha -
ac e ize he non linea beha io o andomly cu ed fibe ed s uc u es
om he mic os uc u e. The o mula ion is s a ed a fini e s ains and
conside s a linea elas ic esponse o he fibe s. Con a y o p e ious
exis ing s uc u al models, a fibe ed ma ix model is buil om he anal-
ysis o he mechanics o a andomly cu ed single fibe . The iso opic
/ aniso opic fibe dis ibu ion, hei e olu ion du ing s aining and he
fluc ua ions o densi y a e modeled a he mic o le el. Upon homoge-
nizing hese fibe ed ma ices, he mac oscopic beha io is cap u ed a
he mac oscale. The e o e, he key aspec o his s udy is o ecognize
he impo ance o he cu ed shape o fibe ed s uc u es on hei h ee
dimensional mechanical beha io unde ensile, comp ession and shea
loads.
2. Single fibe
This sec ion desc ibes he mechanical analysis o a single fibe , om
he gene a ion o andomly cu ed single fibe s, o he ma hema ical
analysis o hei mechanical beha io a ensile and comp ession es s.
2.1. Fibe gene a o
He e we p esen he me hodology o build diffe en 3D andomly
cu ed collagen fibe s (see Fig. 1). We define he inpu pa ame e Las
he dis ance be ween bo h ends o he fibe . This dis ance is disc e ized
in o nsegmen s, in which wa y po ions o fibe s a e gene a ed. An
auxilia y pola sys em (𝑟, 𝜃)is defined along hese segmen s (see Fig. 1a).
The pola angle is andomly selec ed om a uni o m dis ibu ion such
ha 𝜃∈[0,2𝜋). The e o e, he coo dina es o he diffe en disc e ized
poin s along he fibe a e gi en as:
x𝑖=(𝑟cos 𝜃, 𝑟sin𝜃,
𝐿𝑖
𝑛 ), 𝑖=1...𝑛. (1)
The pa ame e n ep esen s in Eq. (1) he numbe o diffe en po -
ions ha compose he fibe , along which he cu ed fibe is wis ed
h oughou i s leng h. On he o he hand, he fibe eccen ici y om he
main axis, 𝑟in Eq. (1), is selec ed om a uni o m dis ibu ion 𝑟∈[0,2𝑅].
𝑅is defined as ollows (see also Fig. 1a):
𝑅=√[(1+𝑃)𝐿
𝑛
]2
−(𝐿
𝑛
)2
,(2)
wi h 𝑃(𝑃≥0) being a ela i e leng h pa ame e o he fibe leng h
e sus he ini ial dis ance be ween he end poin s.
Once he 3-D spa ial poin s 𝒙𝑖ha e been c ea ed (coa se disc e iza-
ion o he fibe ), he cu ed fibe is smoo hed by defining a spline
h ough poin s 𝒙𝑖. This cu e is hen disc e ized (fine disc e iza ion o
he fibe ) in o 𝑁𝐸 elemen s. The final leng h o he fibe 𝐿𝑓is hen
defined as
∑𝑁𝐸
𝑖=1 𝑙𝑖, being 𝑙𝑖 he elemen leng h. Diffe en gene a ed in-
di idual fibe s can be seen in Fig. 1b.
2.2. Ma hema ical o mula ion
The mechanics o he 3D andom fibe s exposed abo e does no hold
an analy ical solu ion. Then, he gene a ed cu ed fibe is disc e ized
in o 𝑁𝐸 Eule -Be nouilli beam fini e elemen s. Each elemen con ains
2 nodes and 6 deg ees o eedom (3 ansla ions and 3 o a ions) pe
node. A linea elas ic beha io o he fibe is assumed. We ollow an
upda ed lag angian o mula ion such ha he fibe geome y is upda ed
a each cu en ime 𝑡o analysis:
x𝑖,𝑡+Δ𝑡=x𝑖,𝑡 +Δu𝑖(3)
Compu e s and S uc u es 310 (2025) 107690
3
J.A. Sanz-He e a, A. Apolina -Fe nandez, A. Jimenez-Ai es e al.
Fig. 2. Exac solu ion o a nonlinea can ile e od subjec ed o a poin dimensionless bending momen
𝑀=𝑀𝐿
𝐸𝐼 e sus i s fini e elemen implemen a ion. The
analy ical solu ion is an a c wi h dimensionless adius
𝑅=1∕
𝑀[45]. The od was disc e ized in o 150 elemen s, and he momen s ep was Δ
𝑀=𝜋∕500.
wi h x𝑖,𝑡+Δ𝑡and x𝑖,𝑡 he posi ion ec o s o node 𝑖a configu a ions 𝑡+Δ𝑡
and 𝑡, espec i ely, and Δu𝑖 he ec o o displacemen s o node 𝑖 om
configu a ion 𝑡 o 𝑡+Δ𝑡. This ec o Δu𝑖is ob ained a e he disc e iza-
ion and assembly o global ec o and ma ix ollowing a s uc u al
ma ix fini e elemen analysis [44]:
𝔸𝑁𝐸
𝑒=1 ΔF𝑒=𝔸𝑁𝐸
𝑒=1 {K𝑒(x𝑒
𝑖,𝑡)⋅Δu𝑒}(4)
being 𝔸 he assembly ope a o ; and ΔF𝑒and Δu𝑒 he inc emen al nodal
ec o o s uc u al o ces and displacemen s/ o a ions a elemen 𝑒,
espec i ely. K𝑒(x𝑒
𝑖,𝑡)is he ma ix o elemen 𝑒compu ed in he configu-
a ion 𝑡wi h elemen nodal coo dina es x𝑒
𝑖,𝑡. A e assembly, he (global)
sys em in Eq. (4) yields,
ΔF=K𝑡⋅Δu(5)
The o al nodal o ce ec o is hen upda ed as:
F𝑡+Δ𝑡=F𝑡+ΔF(6)
The fini e elemen implemen a ion o Eq. (5) was de eloped in Ma -
lab®. The compu e implemen a ion was alida ed in Fig. 2 e sus he
exac solu ion o a nonlinea can ile e od subjec ed o a poin bending
momen .
2.3. Resul s
In his sec ion, we analyze he ensile/comp ession beha io o a
single andomly cu ed fibe o a ange o model pa ame e s, namely,
fibe end poin s leng h 𝐿, fibe diame e 𝑑, fibe elas ic modulus 𝐸;
and fibe leng h o fibe end poin s leng h a io 𝐿𝑓∕𝐿(being his pa-
ame e an indica ion o he wa y/cu ly shape o he fibe ). A baseline
case is selec ed, wi h pa ame e s 𝐿=20μm, 𝑑=0.2μm, 𝐸= 100 MPa
and 𝐿𝑓∕𝐿=1.1. These pa ame e s a e aken as an es ima ion o he o -
de o magni ude ound in he li e a u e o collagen (see Table 1). The
displacemen s o he fibe a e p esc ibed a i s end (bounda y) nodes,
excep in he poin in which a uniaxial s e ch is applied. On he o he
hand, he fibe can eely o a e a he end poin s.
Fig. 3shows he mechanical beha io o a single fibe in uniaxial
s e ch-s ess es s o he ange o analyzed model pa ame e s. The en-
sile egion (𝜆>1) o he cu es in Fig. 3shows a nonlinea beha io
di ided in wo diffe en egimes: fi s a nonlinea pa ollowed by a
s ain s iffening due o fibe s aigh ening. This beha io can be seen in
Table 1
Range o pa ame e s o na u al collagen s uc u es
a ailable in he li e a u e. Dispe sion o pa ame e s
migh be due o diffe en animal o igin.
P ope y Value Re e ences
Linea modulus (MPa) 110-1470 [46]
0-32 [47]
2000-7000 [12]
Fibe leng h (μm) 100-200 [12]
6-18 [46]
10-100 [10]
Fibe diame e (μm) 0.1-0.5 [12]
0.418-0.446 [4]
0.01-0.25 [48]
0.205-0.905 [46]
ideo S1 o he supplemen a y ma e ial, o he baseline fibe . In e es -
ingly, i can be seen in he comp ession egion (𝜆<1) in Fig. 3, ha ou
o mula ion na u ally eco e s he ins abili y (buckling) phenomenon
de eloped by he fibe a comp ession.
On he o he hand, he inc emen al elas ic ene gy o he fibe can be
defined as ollows:
Δ𝑈=
𝑁𝐸
∑
𝑒=1
{F𝑒
𝑡
⋅Δu𝑒}(7)
The o al accumula ed elas ic ene gy o he fibe is hen ob ained as:
𝑈𝑡+Δ𝑡=𝑈𝑡+Δ𝑈(8)
The o al accumula ed elas ic ene gy o he fibe is spli om Eq. (8)
in o axial, bending and o sional ene gies as ollows:
•Axial ene gy is ob ained by selec ing he axial o ces and axial dis-
placemen s p oduc componen s in Eq. (8).
•Bending ene gy is ob ained by selec ing he shea o ces/bending
momen s and shea displacemen s/bending o a ions p oduc com-
ponen s in Eq. (8).
•To sional ene gy is ob ained by selec ing he o que momen and
o que o a ion p oduc componen s in Eq. (8).
Fig. 4shows he diffe en fibe accumula ed ene gies o he uniax-
ial ensile/comp ession s e ch-s ess es s o he baseline single fibe .
Fig. 4addi ionally plo s he diffe en fibe accumula ed ene gies ela-
Compu e s and S uc u es 310 (2025) 107690
4
J.A. Sanz-He e a, A. Apolina -Fe nandez, A. Jimenez-Ai es e al.
Fig. 3. Pa ame ic analysis o he mechanical beha io o elas ic single andomly cu ed fibe . Model pa ame e s a e a ied wi hin he ange ound in he li e a u e
(see Table 1), om a baseline case ha wi h pa ame e s 𝐿=20μm, 𝑑=0.2μm, 𝐸= 100 MPa and 𝐿𝑓∕𝐿=1.1.
Fig. 4. To al, axial, bending and o sional componen s o he accumula ed ene gy (le ) and ela i e ene gy o he o al ene gy ( igh ), o a single andomly cu ed
fibe wi h pa ame e s 𝐿=20μm, 𝑑=0.2μm, 𝐸= 100 MPa and 𝐿𝑓∕𝐿=1.1. Inse fibe s in he figu e on he igh ep esen he de o ma ion s a e o he fibe o
he selec ed s e ch poin s.
i e o he o al ene gy. Acco ding o Fig. 4, he fi s ensile egion o
he cu e ( ep esen ed by poin 1) is go e ned mos ly by fibe bending
(wi h mino con ibu ion o he axial and o sional ene gies), up o a
ce ain s e ch (poin 2) in which he fibe is almos s aigh excep a
small po ion o he fibe . Indeed, he de o ma ion poin 2 ep esen s
he onse o s aigh ening o he fibe , acco ding o he o e all leng h
o he fibe gi en by pa ame e 𝐿𝑓∕𝐿(equal o 1.1 in his case). Once
he fibe is comple ely s aigh (poin 3) he mos impo an con ibu-
ion o ene gy is due o axial de o ma ion. In he comp ession egime
(poin s 4 and 5) he mos impo an con ibu ion is he bending ene gy
as he fibe is ge ing wa y (cu ly). The o sional ene gy is almos neg-
ligible in he diffe en egimes o de o ma ion o he analyzed fibe . I
was checked (da a no shown) ha he o al accumula ed elas ic in e nal
ene gy o he fibe , ob ained om Eqs. (7)-(8), is equal o he ex e nal
ene gy de eloped by he applied o ce in he fibe . No e ha he poin
𝜆=1, a which all he ene gies a e null, was no compu ed in Fig. 4 o
a oid his singula i y.
Keeping in mind ha he fi s ensile/comp ession egions o he
s e ch-s ess cu e o he fibe beha io a e go e ned by bending, and
he s ain s iffening egion by he axial ene gy, he impac o model pa-
ame e s shown in Fig. 3can be easily unde s ood. Consis en ly, we can
see a s iffe beha io (bo h a ensile and comp ession) o inc easing
fibe s iffness (𝐸) and diame e (𝑑); and dec easing leng h (𝐿). Indeed,
in quali a i e e ms, he fibe beha io depends as ∝𝐸𝑑4∕𝐿3in he
Compu e s and S uc u es 310 (2025) 107690
5
J.A. Sanz-He e a, A. Apolina -Fe nandez, A. Jimenez-Ai es e al.
bending egion; and as ∝𝐸𝑑2∕𝐿in he axial (s ain s iffening) egion.
The effec o fibe wa y/cu ly (pa ame e 𝐿𝑓∕𝐿) is seen in Fig. 3as
a s iffe beha io o dec easing 𝐿𝑓∕𝐿(less wa y fibe ) wi h sho e
bending egion and, consequen ly, wi h as e de elopmen o he s ain
s iffening egion. Indeed, he onse o he de elopmen o he s ain s iff-
ening egion is co ela ed wi h pa ame e 𝐿𝑓∕𝐿.
3. Fibe ed ma ix
This sec ion desc ibes he mul iscale mechanical analysis o a fibe ed
ma ix, composed o andomly cu ed fibe s, om he mechanical in e -
ac ion o fibe s in a ep esen a i e olume elemen (RVE). In pa icula ,
he impac o model pa ame e s on he o e all ma e ial beha io , in-
cluding aniso opy, is analyzed a he ma ix mac oscopic le el. I is
assumed ha he RVE is composed solely by fibe s, neglec ing he effec
o su ounding aqueous subs ances in he ma ix.
3.1. Ma ix gene a o
The gene a ion o he RVE is es ablished unde he basis o he defi-
ni ion o RVE, ha is a selec ed olume o he fibe mic os uc u e ha
s a is ically ep esen s he he e ogenei y o he ma ix [49]. Ou RVE
is composed o isola ed fibe s 𝑁𝑓𝑖and c osslink fibe s 𝑁𝑓𝑥, such ha
he o al numbe o fibe s o he RVE is 𝑁𝑓 =𝑁𝑓𝑖+𝑁𝑓𝑥. The e o e,
he olume concen a ion o fibe s is defined as 𝑉𝑐=𝑁𝑓∕𝑉𝑅𝑉 𝐸 , wi h
𝑉𝑅𝑉 𝐸 being he olume o he RVE. The algo i hm o gene a e fibe ed
ma ices p oceeds as ollows:
Box 1: Algo i hm o gene a e fibe ed ma ices.
0. Se he olume concen a ion o fibe s 𝑉𝑐and he olume o
he RVE, and consequen ly, 𝑁𝑓, 𝑁𝑓𝑖and 𝑁𝑓𝑥.
1. Gene a e 𝑁𝑓𝑖isola ed fibe s (see sec ion 2.1).
2. FOR 𝑚=1..𝑁𝑓𝑖
2.1 Randomly selec 𝑁𝑚𝑥 c osslinking nodes o fibe 𝑚.
2.2 FOR 𝑗=1..𝑁𝑚𝑥
i. C ea e a se wi h po en ial connec ing candida e
c osslinking nodes 𝑘, o fibe s 𝑛=1..𝑁𝑓𝑖(𝑛≠𝑚)
such ha 𝐿⋅(1 − 𝜖)≤𝐿𝑗−𝑘≤𝐿⋅(1 + 𝜖).
ii. Randomly selec node 𝑗 om he se .
iii. Gene a e a c osslinking fibe om nodes 𝑗−𝑘(wi h
leng h 𝐿) as in sec ion 2.1.
END FOR.
END FOR.
We se 𝑁𝑚𝑥 =1 in Box 1, meaning ha a c osslinking fibe is
gene a ed pe isola ed fibe . The e o e, 𝑁𝑚𝑥 ep esen s he deg ee o
c osslinking o he ma ix. On he o he hand, he endpoin s leng h o
he c osslinking fibe is se o he same leng h 𝐿o he isola ed fibe ,
acco ding o i em iin Box 1, up o a ole ance 𝜖which is se o 0.01 in
ou code.
3.2. Mul iscale o mula ion
The o e all mac oscopic mechanical beha io o a fibe ed ma ix and
i s e olu ion, is ob ained om he mic omechanical in e ac ion o fibe s
wi hin he RVE in a mul iscale ashion (see Fig. 5). In he mac oscale,
we define he (Cauchy) s ess enso and he (loga i hmic) s ain enso
as 𝝈𝑀
𝐼,𝑡(𝐗𝐼,𝑡),E𝑀
𝐼,𝑡(𝐗𝐼,𝑡), espec i ely, in a mac oscopic (Gauss) ma e ial
poin 𝐼o he mac oscale (wi h coo dina es 𝐗𝐼,𝑡), o load ime 𝑡(see
Fig. 5). The ( ue) Cauchy’s s ess componen s a e used in he o mu-
la ion o con enience and subsequen alida ion wi h expe imen al
ou comes. On he o he hand, in he mic oscale, we define he a i-
ables F𝑚
𝑡(𝐱𝑖,𝑡),u𝑚
𝑡(𝐱𝑖,𝑡)associa ed o an Eule -Be noulli beam in a ma e-
ial poin (node) 𝑖o he mic oscale (wi h coo dina es 𝐱𝑖,𝑡) o load ime
Fig. 5. Schema ics o he p oposed mul iscale o mula ion o fibe ed ma ices.
𝑡. F𝑚
𝑡(𝐱𝑖,𝑡) ep esen s he ec o o nodal s uc u al o ces ( o ces and mo-
men s), whils u𝑚
𝑡(𝐱𝑖,𝑡)is he ec o ha con ains nodal displacemen s
and o a ions (see Fig. 5).
As s a ed in he analysis o he single fibe , we ollow an upda ed
lag angian app oach such ha he configu a ion (fini e elemen mesh)
is upda ed a each load ime 𝑡. Consequen ly, he mac oscopic beha io
o he ma ix is geome ically nonlinea and es ablished a fini e s ains
h ough he loga i hmic s ain enso . The loga i hmic s ain enso is
ob ained by in eg a ing he s ain a e, nume ically, in a ma e ial ame
o e e ence [50]:
E𝑀
𝐼,𝑡+Δ𝑡=ΔR⋅E𝑀
𝐼,𝑡
⋅ΔR𝑇+ΔE𝑀
𝐼(9)
whe e E𝑀
𝐼,𝑡+Δ𝑡and E𝑀
𝐼,𝑡 a e he o al s ains a inc emen s 𝑡+Δ𝑡and 𝑡, e-
spec i ely; ΔRis he inc emen al o a ion enso ; and ΔE𝑀
𝐼is he o al
s ain inc emen om inc emen 𝑡 o 𝑡+Δ𝑡. Fu he mo e, he mac o-
scopic s ain ene gy densi y om inc emen 𝑡 o 𝑡+Δ𝑡is defined as:
Δ𝑈𝑀
𝐼=𝝈𝑀
𝑡+Δ𝑡,𝐼
⋅ΔE𝑀
𝐼(10)
Analogously, he kinema ics a he mic oscale is defined as:
u𝑚
𝑡+Δ𝑡=u𝑚
𝑡+Δu𝑚(11)
whe e u𝑚
𝑡+Δ𝑡and u𝑚
𝑡a e he o al displacemen s/ o a ions ec o s a in-
c emen s 𝑡+Δ𝑡and 𝑡, espec i ely; and Δu𝑚is he o al displacemen
inc emen ec o om inc emen 𝑡 o 𝑡+Δ𝑡. The displacemen s inc e-
men ec o is defined a nodes 𝑘o he mic oscale h ough he s ain
inc emen a he mac oscale as:
Δu𝑚
𝑘=ΔE𝑀
𝐼
⋅x𝑘,𝑡 (12)
whe e x𝑘,𝑡 a e he (upda ed) mic oscopic coo dina es, a load ime 𝑡, o
he poin s 𝑘in which he mac oscopic s ain is affinely ansmi ed o
he mic oscale. We assume ha poin s 𝑘a e he endpoin s (kno s) o
he fibe s. A e disc e iza ion and FE assembly, he global nodal o ces
ec o is upda ed as ollows,
F𝑚
𝑡+Δ𝑡=F𝑚
𝑡+ΔF𝑚(13)
Using Eq. (5) in (13) yields:
F𝑚
𝑡+Δ𝑡=F𝑚
𝑡+𝐊𝑡⋅Δu𝑚(14)
Finally, he o al mic oscopic s ain ene gy densi y om inc emen
𝑡 o 𝑡+Δ𝑡is defined as:
Δ𝑈𝑚=1
𝑉𝑅𝑉 𝐸
𝑁
∑
𝑖=1
𝐹𝑚
𝑖,𝑡+Δ𝑡
⋅Δ𝑢𝑚
𝑖(15)
Compu e s and S uc u es 310 (2025) 107690
6
J.A. Sanz-He e a, A. Apolina -Fe nandez, A. Jimenez-Ai es e al.
whe e 𝑁is he numbe o nodes o he mic os uc u e. 𝐹𝑚
𝑖,𝑡+Δ𝑡and Δ𝑢𝑚
𝑖
ep esen he 𝑖componen s o ec o s F𝑚
𝑡+Δ𝑡and Δu𝑚, espec i ely. On
he o he hand, we use he Hill’s equali y o ene gies om he ansi ion
om he mac o o he mic oscale [49,51]
Δ𝑈𝑀
𝐼=Δ𝑈𝑚(16)
Thus, using (10) and (15) in (16) yields,
𝝈𝑀
𝑡+Δ𝑡,𝐼
⋅ΔE𝑀
𝐼=1
𝑉𝑅𝑉 𝐸
𝑁
∑
𝑖=1
𝐹𝑚
𝑖,𝑡+Δ𝑡
⋅Δ𝑢𝑚
𝑖(ΔE𝑀
𝐼)(17)
No e ha in Eq. (15), Δ𝑢𝑚
𝑖(ΔE𝑀
𝐼) ep esen s he inc emen al nodal
displacemen s/ o a ions o he mic os uc u e, as a esul o he solu ion
o he disc e e fini e elemen sys em ΔF𝑚=K𝑡⋅Δu𝑚wi h p esc ibed
displacemen s a disc e ized nodes o kno s 𝑘acco ding o Eq. (12).
The e o e, he diffe en 𝑘𝑙 componen s o he s ess enso 𝜎𝑀
𝑡+Δ𝑡,𝐼,𝑘𝑙 a
mac oscopic poin 𝐼(see Fig. 5) a e ob ained using Eq. (17) as ollows:
𝜎𝑀
𝑡+Δ𝑡,𝐼,𝑘𝑙 =1
𝑉𝑅𝑉 𝐸
𝑁
∑
𝑖=1
𝐹𝑚
𝑖,𝑡+Δ𝑡
⋅Δ𝑢𝑚
𝑖(Δ𝐸𝑀
𝐼,𝑘𝑙)⋅
1
Δ𝐸𝑀
𝐼,𝑘𝑙
(18)
Eq. (18) is equi alen o,
𝝈𝑀
𝑡+Δ𝑡,𝐼 =1
𝑉𝑅𝑉 𝐸
𝑁
∑
𝑖=1
𝐹𝑚
𝑖,𝑡+Δ𝑡
⋅Δ𝑢𝑚
𝑖(𝕀)(19)
wi h 𝕀being he uni enso :
𝕀=1
2(𝛿𝑖𝑘𝛿𝑗𝑙 +𝛿𝑖𝑙 +𝛿𝑗𝑘)(20)
Δ𝑢𝑚
𝑖(𝕀) ep esen s he inc emen al nodal displacemen s/ o a ions o
he mic os uc u e, as a esul o he solu ion o he disc e e fini e
elemen sys em ΔF𝑚=K𝑡⋅Δu𝑚wi h p esc ibed displacemen s a dis-
c e ized nodes o kno s 𝑘 o uni s ain enso s 𝕀, ollowing Eq. (12).
Using he decomposi ion o he nodal o ces ec o (13) in (18), he
inc emen al s ess enso can be ob ained as:
Δ𝝈𝑀
𝐼=1
𝑉𝑅𝑉 𝐸
𝑁
∑
𝑖=1
(𝐾𝑡,𝑖𝑗 ⋅Δ𝑢𝑚
𝑗(ΔE𝑀
𝐼))⋅Δ𝑢𝑚
𝑖(𝕀)(21)
The mac oscopic angen s iffness enso can be hen ob ained as:
ℂ𝑡=𝜕Δ𝝈𝑀
𝐼
𝜕ΔE𝑀
𝐼
=1
𝑉𝑅𝑉 𝐸
𝑁
∑
𝑖=1
(𝐾𝑡,𝑖𝑗 ⋅Δ𝑢𝑚
𝑗(𝕀))⋅Δ𝑢𝑚
𝑖(𝕀)(22)
Box 2 summa izes he main s eps o he de eloped mul iscale algo-
i hm.
Box 2: Algo i hm o compu e mac oscopic s ess and s ain quan-
i ies om mic oscopic fibe ed ma ices.
1. Mac oscale: E𝑀
𝐼,𝑡, 𝝈𝑀
𝐼,𝑡 a ailable a load inc emen 𝑡.
2. Mac oscale: Se nex s ain inc emen ΔE𝑀
𝐼.
3. Mic oscale: Ob ain Δ𝑢𝑚
𝑗(ΔE𝑀
𝐼)and Δ𝑢𝑚
𝑖(𝕀) om esul o
he solu ion o he sys em ΔF𝑚=K𝑡⋅Δu𝑚wi h p esc ibed dis-
placemen s a nodes 𝑘 o s ain enso ΔE𝑀
𝐼and uni s ain
enso s 𝕀, espec i ely, ollowing Eq. (12).
4. Mic oscale: Upda e mic oscopic coo dina es as x𝑖,𝑡+Δ𝑡=x𝑖,𝑡 +
Δu𝑚
𝑖(ΔE𝑀
𝐼).
5. Mac oscale: Compu e Δ𝝈𝑀
𝐼acco ding o Eq. (21). Com-
pu e he mac oscopic angen s iffness enso ℂacco ding o
Eq. (22).
6. Mac oscale: Upda e mac oscopic quan i ies as 𝝈𝑀
𝐼,𝑡+Δ𝑡=
𝝈𝑀
𝐼,𝑡 +Δ𝝈𝑀
𝐼, E𝑀
𝐼,𝑡+Δ𝑡=E𝑀
𝐼,𝑡 +ΔE𝑀
𝐼
7. 𝑡←𝑡+Δ𝑡. GO TO 1.
Finally, he inc emen al elas ic in e nal mic oscopic ene gy o a fibe
𝑓o he mic os uc u e o he ma ix, Δ𝑈𝑚,𝑓 , is compu ed as,
Δ𝑈𝑚,𝑓 =
𝑁𝐸𝑓
∑
𝑒=1 {F𝑚,𝑒
𝑡
⋅Δu𝑚,𝑒}(23)
whe e F𝑚,𝑒
𝑡and Δu𝑚,𝑒 a e he o al s uc u al o ces and displace-
men s/ o a ions nodal ec o s a elemen 𝑒, espec i ely. 𝑁𝐸𝑓is he
numbe o elemen s o he fibe 𝑓. The o al accumula ed elas ic ene gy
o each fibe 𝑓o he mic os uc u e o he ma ix is ob ained as:
𝑈𝑚,𝑓
𝑡+Δ𝑡=𝑈𝑚,𝑓
𝑡+Δ𝑈𝑚,𝑓 .(24)
On he o he hand, he o al accumula ed mic oscopic s ain ene gy
densi y is compu ed as,
𝑈𝑚
𝑡+Δ𝑡=𝑈𝑚
𝑡+Δ𝑈𝑚(25)
wi h Δ𝑈𝑚ob ained om Eq. (15).
3.3. Resul s: ensile/comp ession
In his sec ion, diffe en fibe ed mic os uc u es—wi h diffe en mi-
c os uc u al cha ac e is ics—a e subjec ed o a mac oscopic ensile and
comp ession s ess s a e. The s ess-s e ch cu es o he i ual es s a e
hen ob ained ollowing he mul iscale app oach exposed in sec ion 3.2.
Fo his pu pose, a e ical s e ch eloci y
𝜆(𝑌-di ec ion) is p e-
sc ibed. The mac oscopic inc emen al s ain enso (poin 2 o Box 2 o
he mul iscale o mula ion) is ob ained om he solu ion o :
Δ𝝈𝑀=ℂ𝑡⋅ΔE𝑀(26)
ΔE𝑀is ob ained om he solu ion o he sys em in Eq. (26) wi h
uniaxial (mixed bounda y condi ions) s ess s a e. In his sense, Δ𝐸𝑀
𝑌=
Δ𝑡⋅
𝜆 o he mac oscopic inc emen al s ain enso , and Δ𝜎𝑀
𝑋=Δ𝜎𝑀
𝑍=
Δ𝜎𝑀
𝑋𝑌 =Δ𝜎𝑀
𝑋𝑍 =Δ𝜎𝑀
𝑌𝑍 =0.
Fig. 6shows he mechanical beha io o he conside ed ma ices in
uniaxial s e ch-s ess es s. In his figu e, diffe en andomly cu ed
fibe ed ma ices, a diffe en fibe e sus olume concen a ion (pa am-
e e 𝑉𝑐), we e gene a ed o he ange o analyzed model pa ame e s.
These pa ame e s we e a ied wi hin he o de o magni ude o eal
biological fibe ed ma ices (see Table 1), om a baseline case wi h
pa ame e s 𝐿=30μm, 𝑑=0.4μm, 𝐸= 150 MPa, 𝐿𝑓∕𝐿=1.15 and
𝑉𝑐=2.05 ⋅106fibe s∕𝑚𝑚3. The mic os uc u es we e compu ed using
he ma ix gene a o conside ing ha fibe s can be dis ibu ed and o i-
en ed eely along he RVE wi hou any p e e en ial di ec ion, such ha
he mic os uc u es a e assumed iso opic.
Simila ly o he single fibe , he ensile egion (𝜆>1) o he cu es
in Fig. 6shows a linea beha io a small s ains ollowed by a nonlinea
beha io . In his case, he nonlinea egime is he consequence o bo h
fibe s aigh ening and fibe alignmen wi h he di ec ion o load. Be-
sides, he loss o s iffness due o he ins abili y (buckling) phenomenon
de eloped by he fibe s a comp ession, can be seen in he comp es-
sion egion (𝜆<1) in Fig. 6. Acco ding o Fig. 6, and analogously o he
single fibe beha io , he o e all ensile/comp ession beha io o he
ma ices is s iffe o inc easing fibe s iffness (𝐸), fibe diame e (𝑑)
and concen a ion o fibe s. On he o he hand, he o e all ensile/com-
p ession beha io o he ma ices is so e o inc easing fibe leng h (𝐿)
and fibe cu a u e (𝐿∕𝐿𝑓).
The o al accumula ed s ain ene gy densi y de eloped by he mi-
c os uc u e, compu ed om Eq. (25), can be seen in Fig. 7 o he uniax-
ial ensile/comp ession s e ch-s ess es . Resul s a e he e o e shown
as he (accumula ed) mic oscopic con ibu ion o all he fibe s o he
ma ix along he de o ma ion pa h, o he conside ed baseline ma ix
defined abo e. No e ha his quan i y is equal o he o al accumula ed
s ain ene gy densi y de eloped a he mac os uc u e (Eq. (10)), by
i ue o Eq. (16). In Fig. 7, he axial, bending and o sional compo-
nen s o he o al accumula ed s ain ene gy densi y a e also shown.
Compu e s and S uc u es 310 (2025) 107690
7
J.A. Sanz-He e a, A. Apolina -Fe nandez, A. Jimenez-Ai es e al.
Fig. 6. Tensile/comp ession es s. Pa ame ic analysis o he mechanical beha io o elas ic andomly cu ed fibe ed ma ices. Model pa ame e s a e a ied wi hin he
ange ound in he li e a u e (see Table 1), om a baseline case wi h pa ame e s 𝐿=30μm, 𝑑=0.4μm, 𝐸= 150 MPa, 𝐿𝑓∕𝐿=1.15 and 𝑉𝑐=2.05 ⋅106fibe s∕𝑚𝑚3.
These quan i ies we e ob ained om Eq. (25), a e selec ing he ax-
ial, bending and o sional componen s simila ly o he single fibe case
(sec ion 2.3). As in he single fibe case, i can be seen in his figu e he
asymme y o he o al s ain ene gy densi y a ensile (𝜆>1) e sus he
comp ession egime (𝜆<1). Indeed, ma ix comp ession induces fibe
ins abili ies (buckling) whe eas ma ix ensile p oduces fibe un olling
and fibe alignmen wi h load (s ain s iffening). Fig. 7addi ionally plo s
he diffe en accumula ed s ain ene gy densi ies ela i e o he o al
s ain ene gy densi y.
On he o he hand, Fig. 8shows he accumula ed mul iscale ene gies
in he fibe s o he mic os uc u e o he ma ix o diffe en ensile/com-
p ession de o ma ion poin s, as defined in Fig. 7. These quan i ies we e
compu ed om Eq. (24) a e selec ing he axial, bending and o sional
componen s. The ensile/comp ession beha io o he baseline ma ix,
he e olu ion o i s mic os uc u e, and bo h he axial, bending and o -
sional accumula ed ene gies in he fibe s; can be seen in ideos S2, S3
and S4 o he supplemen a y ma e ial. I can be seen in Fig.
7 ha bend-
ing ene gy domina es in all he analyzed de o ma ion egime, excep he
s ain s iffening egion. A he ensile egion (𝜆>1) his ene gy is em-
ployed in fibe un olling and fibe alignmen wi h load, as co obo a ed
in he mic os uc u e e olu ion by Fig. 8and ideo S3. A he comp es-
sion egion (𝜆<1) he bending ene gy is employed in fibe buckling, as
co obo a ed in he mic os uc u e e olu ion (see also Fig. 8and ideo
S3). Mo eo e , he o sional ene gy is almos negligible in all he ana-
lyzed de o ma ion ange, acco ding o Figs. 7and 8.
3.4. Resul s: simple shea
Analogously o he p e ious sec ion, diffe en fibe ed mic os uc-
u es (assumed iso opic) a e subjec ed o a mac oscopic simple shea
s ess s a e. In his case, a shea s ain eloci y (𝛾) is p esc ibed, such
ha he mac oscopic inc emen al s ain enso (poin 2 o Box 2 o he
Compu e s and S uc u es 310 (2025) 107690
8
J.A. Sanz-He e a, A. Apolina -Fe nandez, A. Jimenez-Ai es e al.
Fig. 7. Mul iscale ene gies: Mac oscopic. To al, axial, bending and o sional accumula ed s ain ene gy densi ies (le ), and ela i e o he o al s ain ene gy densi y
( igh ), de eloped by he mic os uc u e o ensile/comp ession es s o he baseline ma ix.
Fig. 8. Mul iscale ene gies: Mic oscopic (fibe scale). Plo in he mic os uc u e o he axial (fi s ow), bending (second ow) and o sional ( hi d ow) o he
accumula ed ene gies in he fibe s o he baseline ma ix o ensile/comp ession es s. Resul s a e shown a diffe en de o ma ion poin s (1,2,3) defined in Fig. 7.
Compu e s and S uc u es 310 (2025) 107690
9
J.A. Sanz-He e a, A. Apolina -Fe nandez, A. Jimenez-Ai es e al.
Fig. 9. Simple shea es s. Pa ame ic analysis o he mechanical beha io o elas ic andomly cu ed fibe ed ma ices. Model pa ame e s a e a ied wi hin he ange
ound in he li e a u e (see Table 1), om a baseline case wi h pa ame e s 𝐿=30μm, 𝑑=0.4μm, 𝐸= 150 MPa, 𝐿𝑓∕𝐿=1.15 and 𝑉𝑐=2.05 ⋅106fibe s∕𝑚𝑚3.
mul iscale o mula ion) is ully p esc ibed as a null alue o all compo-
nen s excep o Δ𝐸𝑀
𝑋𝑌 =Δ𝑡⋅𝛾∕2.
Fig. 9shows he mechanical beha io o he conside ed ma ices in
simple shea s ess es s. In his figu e, we assumed he same ange o
a ia ion o he pa ame e s, as well as he same baseline case, han he
ensile/comp ession es s. The e o e, he same mic os uc u es (RVEs)
defined o he ensile/comp ession es s in he p e ious sec ion, a e
used he e o simula e simple shea es s.
The shea beha io o he ( i ually) es ed ma ices, Fig. 9, shows
also a nonlinea p ofile. Fi s , a linea egion a small s ains ollowed by
a nonlinea egion, also as a consequence o bo h fibe s aigh ening (un-
olling) and fibe alignmen wi h he di ec ion o load. Analogously o
he ensile es s, he shea beha io o he ma ices is s iffe o inc eas-
ing fibe s iffness (𝐸), fibe diame e (𝑑) and concen a ion o fibe s. On
he o he hand, he o e all ensile/comp ession beha io o he ma ices
is so e o inc easing fibe leng h (𝐿) and fibe cu a u e (𝐿∕𝐿𝑓).
As in he p e ious sec ion, he accumula ed s ain ene gy densi ies
de eloped by he mic os uc u e o he shea es can be seen in Fig. 10.
In his figu e, esul s a e shown as he (accumula ed) mic oscopic con-
ibu ion o all he fibe s o he ma ix along he de o ma ion pa h, o
he conside ed baseline ma ix. The axial, bending and o sional compo-
nen s o he s ain ene gy densi ies we e compu ed analogously o he
p e ious sec ion. Fig. 10 addi ionally plo s he diffe en accumula ed
s ain ene gy densi ies ela i e o he o al s ain ene gy densi y.
On he o he hand, Fig. 11 shows he accumula ed mul iscale ene -
gies in he fibe s o he mic os uc u e o he ma ix o diffe en shea
de o ma ion poin s, as defined in Fig. 10. The shea beha io o he
baseline ma ix, he e olu ion o i s mic os uc u e, and bo h he axial,
bending and o sional accumula ed ene gies in he fibe s; can be seen
