b i o s y s t e m s e n g i n e e r i n g 1 8 6 ( 2 0 1 9 ) 3 0 7 e3 2 2
Available online at www.sciencedirect.com
ScienceDirect journal homepage: www.elsevier.com/locate/issn/15375110
Research Paper
On the nonlinear viscoelastic behaviour of fresh and dried oil palm mesocarp fibres Ahmad Tarmezee Talib a, Chu Chang Jie a, Mohd Afandi P Mohammed a,*, Azhari Samsu Baharuddin a, Mohd Noriznan Mokhtar a, Minato Wakisaka b a
Department of Process and Food Engineering, Faculty of Engineering, Universiti Putra Malaysia, 43400 UPM Serdang, Selangor, Malaysia b Graduate School of Life Science and Systems Engineering, Kyushu Institute of Technology, 2-4 Hibikino, Wakamatsu-ku, Kitakyushu, 808-0196, Japan
article info
Investigation on viscoelastic behaviour of oil palm fibres was reported through experi-
Article history:
mental and finite element modelling study. The experimental work through tensile stress
Received 13 January 2019
relaxation and cyclic tests revealed time-dependent behaviour and damage within the oil
Received in revised form
palm fibres. From the former test results, stresses of fresh fibres reduced more than the
7 August 2019
dried ones after 1 second relaxation, whereas increasing damage was observed under
Accepted 17 August 2019
larger deformations from the latter test results. Finite element modelling results using Prony series viscoelastic model with damage function only agreed with small deformation test, whereas Parallel Rheological Framework viscoelastic model agreed with large defor-
Keywords:
mation test. The former model can be used for biodegradation study, which does not
oil palm fibres
involve large deformation, whereas the latter model is suitable for biocomposites study
finite element analysis
under large deformation.
mechanical testing
1.
Introduction
Oil palm fibres are lignocellulosic biomass by-product produced after oil extraction process in oil palm mill. Considering that oil palm fibres are the most produced biomass in Malaysia, it is important to understand in detail the behaviour of this biomass. The fibres can be utilised in biocomposites and biocompost applications, where detailed mechanical understanding of the fibre-matrix interface for the development of biocomposites (Hanipah, Mohammed, & Baharuddin,
* Corresponding author. E-mail address:
[email protected] (M.A. P Mohammed). https://doi.org/10.1016/j.biosystemseng.2019.08.010 1537-5110/© 2019 IAgrE. Published by Elsevier Ltd. All rights reserved.
© 2019 IAgrE. Published by Elsevier Ltd. All rights reserved.
2016) and cellulose-hemicellulose-lignin interaction for biodegradation (Omar et al., 2016) are crucial to be understood. An example of complex mechanical behaviour of the fibres is viscoelasticity, which was reported through experimental studies on oil palm fibres by (Hanipah, Xiang, Mohammed, & Baharuddin, 2017; Sreekala, Kumaran, & Thomas, 1997, 2001; Xiang, Hanipah, Mohammed, Baharuddin, & Lazim, 2015). A viscoelastic material exhibits both viscous and elastic characteristics when undergoing deformation. Viscoelastic fibres influenced the overall performance of a biocomposites system under different loading conditions (sudden impact or small
308
b i o s y s t e m s e n g i n e e r i n g 1 8 6 ( 2 0 1 9 ) 3 0 7 e3 2 2
Nomenclature Knowles model shear modulus (MPa) Knowles model hardening parameter Knowles model material constant Inverse Bulk modulus for Knowles model (MPa1) B volume preserving part of the left Cauchy strain tensor (mm mm1) trðÞ Tensor trace function I1 First invariant of strain tensor Deformation gradient tensor (mm mm1), F, F Distortional gradient tensor J Volume change (mm3) Kronecker delta dij t, S Kirchoff stress, Second Piola Kirchoff stress (MPa) DEVðÞ Deviatoric function C Volume preserving part of the right Cauchy strain tensor (mm mm1) e Material tangent modulus for the Knowles C potential (MPa) ex dex a ðaÞ, db ðbÞ Un-continuous and continuous damage variables aðtÞ, bðtÞ Un-continuous and continuous damage parameters WðsÞ Strain energy function (MPa) Un-continuous and continuous damage ha , hb constants ∞ Un-continuous and continuous damage d∞ a , db constants signðÞ Absolute value of function Viscoelastic term from convolution integral Hvd j Relaxation constant gj Fixed time constant (e.g. 0.1, 1 10s …) (s) xj ith volume fraction of the network si Ogden model material constant aogden li The deviatoric principal stretches ogden Ogden model bulk modulus (MPa) D1 mogden Ogden model shear modulus (MPa) cr ε_ Equivalent creep strain rate (s1) cr ε Equivalent creep strain (s1) q~ Equivalent deviatoric Kirchoff stress (MPa) A, mv , nn PRF power law material parameters Equivalent stress (plastic stress) (MPa) sp Equivalent plastic strain (mm mm1) εp h Stress softening function ~m Deviatoric part of the strain energy density of U dev the finite strain viscoelastic model r, mm , bm Stress softening material parameters Maximum stress before relaxation occur (MPa) st¼0 Stress at any relaxation time (MPa) st¼n Normalised stress (MPa) snormalised t εe , εp Elastic and plastic true strain at zero stress (mm mm1) mknowles n Bknowles Dknowles 1
loads applied over time), as well as long-term biocomposites mechanical performance (Sreekala, Kumaran, Joseph, & Thomas, 2001). On the other hand, understanding the effect of viscoelasticity towards fibre's biodegradation can provide an accurate in-situ biodegradation study (as conducted by Omar, Hafid, Baharuddin, Mohammed, and Abdullah (2017)), which involves very small loads applied by microbial agent to degrade/reduce the volume of a single fibre (observed using Xray micro-tomography). Two important mechanical tests involving the study of viscoelasticity are stress relaxation and cyclic deformation tests, where the former test is conducted by allowing a fibre sample to deform to a specified deformation under tension, before reduction of force is recorded (while the sample is held fixed over time). The cyclic test is conducted by loading and unloading (stretching and un-stretching) a fibre sample for a few cycles at a constant or increasing deformation. From the cyclic test, damage due to deformation can be identified from the unloading-reloading stressestrain curves, which can be modelled using finite element method. Damage function within a viscoelastic finite element model can be studied either through continuum-based or microstructure-based approach. The former approach includes a damage function/parameter into a constitutive stressestrain model that can be implemented using finite element analysis (FEA) (Kaliske, Nasdala, & Rothert, 2001). The latter approach considers micromechanics approach by taking into consideration the microstructure of a material studied (e.g. material geometry from microscopy analyses). For example, Mohammed, Tarleton, Charalambides, and Williams (2013) modelled damage in terms of cohesive damage between the interface of viscoelastic matrix and viscoplastic filler of bread dough, with geometry of the sample obtained from cryo-SEM image. This is quite expensive to implement, considering large computational power and time needed especially for complex multiphase materials. Therefore, the continuum-based damage model has been a popular approach by researchers due to its simplicity and easiness to implement. An example of such model is the parallel rheological framework (PRF) by Hurtado, Lapzcyk, & Govindarajan (2013), which considers the inclusion of stress softening and plasticity functions into a nonlinear viscoelastic model. However, Omar et al. (2016) showed that there is a need to understand the contribution of different PRF model parameters, i.e. viscoelastic, stress softening and plasticity functions contribution to the unloading-reloading curves of a cyclic tensile result. On the other hand, Hanipah et al. (2016) modelled oil palm fibres using stress-softening function (Ogden & Roxburgh, 1999) within a Prony series finite viscoelastic model. However, the modelling results only considered the unloading part of a single cyclic test. This suggested that the oil palm fibres damage due to deformation can be modelled either through stress softening and/or plasticity function. It should be noted that the oil palm fibres mechanics works reported before (Hanipah et al., 2016; Omar et al., 2016; Omar,
b i o s y s t e m s e n g i n e e r i n g 1 8 6 ( 2 0 1 9 ) 3 0 7 e3 2 2
Hafid, Mohammed, Baharuddin, & Abdullah, 2017; Omar, Mohammed, & Baharuddin, 2014; Xiang et al., 2015) were all obtained from dried oil palm fibres that went through milling, washing and drying processes, which changed the structural behaviour of the oil palm fibres. It is important to compare the behaviour of dried and fresh oil palm fibres (before any mechanical process) so that further understanding on the underlying damage of the fibres can be provided. In this work, experimental work on fresh and dried oil palm mesocarp fibres is conducted through mechanical tests and microscopy observation. Comparison between dried and fresh fibres is provided to investigate the effect of processing towards integrity and nonlinear behaviour of oil palm fibres. In addition, we attempt to simulate oil palm fibres cyclic deformation using viscoelastic and nonlinear viscoelastic models.
2.
Materials and methods
2.1.
Experimental work
Oil palm mesocarp fibres (OPMF) (as illustrated in Fig. 1) from the variety of Elaeis guineensis Jacq. were obtained from Seri Ulu Langat Palm Oil Mill in Dengkil, Selangor, Malaysia. The fresh fibres were taken directly from a mature fresh oil palm fruitlet. The fibres were carefully removed from the fruit using a scalpel to minimise damage. On the other hand, processed mesocarp fibres were obtained from the oil palm mill and were kept in a freezer (at 20 C) to avoid growth of fungae on the surface of the fibres. For tensile test preparation, the processed fibres were thoroughly washed with detergents to eliminate remnant of palm oil and dirt. The fibres were rinsed with tap water and final rinse was performed with distilled water. The fibres were then oven dried at 105 C overnight to remove excess moisture (dried fibres) and the moisture content for dried fibres were kept below 5% by placing them inside desiccators for storage. The determination of moisture content for mesocarp fibres were determined following Khalil, Ismail, Rozman, and Ahmad (2001).
309
The morphology of the fibres was observed using a portable digital USB microscope (Dino-Lite 4113 Series, Taiwan). Mechanical tests were conducted using texture analyser (TA-XT, Stable Micro System Ltd., UK) according to ASTM D3379-75 standard (Standard Test Method for Tensile Strength and Young's Modulus for High Modulus Single Filament Fibers) with modification following the method by Hanipah et al. (2017), as follows. A 20 mm length of fibre with average diameter between 0.1 mm and 0.25 mm was glued on a Cshaped paper before the sample was loaded to the texture analyser and the paper was cut prior to the mechanical test where the samples were clamped using a vice clamp (TA.XT A/TG tensile grips). Uniaxial tension tests (under stress relaxation and cyclic modes) were performed by pulling each end of the sample in an opposite direction. Based on the uniaxial tension results by Hanipah et al. (2017) (where the stressestrain curve was approximately linear at 5%, became plateau at 10%, before fractured at 15%), as well as our preliminary tests where most of the fibres used failed before 20%, the cyclic and stress relaxation tests in this experiment then were conducted at 5%, 10% and 15% deformations. Cyclic tests were performed for five loading-unloading cycles at each strain deformation. All tests were commenced at a constant crosshead speed of 1 mm s1 with at least five replications for each deformation. To investigate possible anisotropic behaviour, tensile tests were conducted by twisting the fibre first by up to 1800 in clockwise and anti-clockwise positions (different samples) before commencing with tensile tests at 1 mm s1. This is following the method proposed by Placet Cisse, & Boubakar (2014). The twisting procedure was conducted manually after the sample was attached to the top and bottom grips (by rotating the top grip by 1800 while fixing the bottom grip). At least four samples were used for each test conditions.
2.2. Viscoelastic model (Prony series) with damage function A finite viscoelastic model is used to simulate deformation of oil palm fibres, but anisotropic behaviour is not considered
Fig. 1 e (a) Transverse section and (b) longitudinal section of oil palm fruitlet. Fibres from mesocarp region of oil palm fruitlet were used in this study.
310
b i o s y s t e m s e n g i n e e r i n g 1 8 6 ( 2 0 1 9 ) 3 0 7 e3 2 2
Fig. 2 e Prony series represented by Maxwell springs and dampers.
here, since to the best of the authors’ knowledge, no experimental work on this is available for oil palm fibres. However, considering lignocellulosic materials like wood showed anisotropy due to cellulose micro-fibrils arrangement (Baskin, € nsson, 2016), it is 2005; Baskin & Jensen, 2013; Braybrook & Jo most likely oil palm fibres would have similar behaviour (as will be shown later in the experimental results). This will be the interest for further improvement on oil palm fibres characterisation, which requires the use of a very short fibre (to prevent buckling during tensile testing) and a specialised micro-mechanical tester to conduct tests under compression, shear and torsion (or pre-rotation before tensile) modes. The viscoelastic model can be represented through the Prony series consisting of a series of Maxwell elements (springs and dampers) connected in parallel with a spring, as shown in Fig. 2. Elasticity part of the finite viscoelastic model is represented through the Knowles strain energy function (Knowles, 1977; Suchocki, 2011): W¼
mknowles 2Bknowles
n Bknowles ðI1 3Þ 1 þ U 1þ n
(1)
where mknowles is the shear modulus, n is a “hardening” parameter, Bknowles is a material constant and Uis volumetric stored1 ðJ 1Þ2 with Dknowles is energy function defined as U ¼ Dknowles 1 1 the inverse of bulk modulus. When Bknowles and n are equal to 1, the neo-Hookean strain energy density is recovered. Note that the neo-Hookean model is not used in this work due to the model is not able to capture nonlinear elastic behaviour in the range larger than 5% strain (oil palm fibres can deform up to 25% strain (Hanipah et al., 2017). The advantage of this model (hereafter referred to as the Knowles model) is that it uses only four material constants and only one invariant (first invariant), which simplifies the material parameters identification process. The first invariant in Equation (1) is defined as:
I1 ¼ trðBÞ
(2)
where B is the volume preserving part of the left Cauchy strain tensor: B ¼ FF
T
(3)
and the volume change J is defined as: J ¼ detðFÞ
(4)
Bis related to the left Cauchy-Green tensor through: B ¼ Bij ¼
Bij J2=3
(5)
and similarly, the relationship between deformation gradient, F and distortion gradient, F can be defined as: F ¼ J1=3 F
(6)
Equation (1) can be expressed in terms of the Kirchhoff stress as: n1 Bknowles 1 2 ðI1 3Þ Bij I1 dij þ JðJ 1Þdij te ¼ mknowles 1 þ 3 D1 n (7) which can then be used to obtain the second Piola Kirchoff stress: Se ¼ F1 te FT
(8) e
Deviatoric stress, DEVðS Þ , which is needed to model damage in hyperelastic model is then obtained through: 1 1 DEVðÞ : ¼ ðÞ ½C:ðÞC 3
(9)
The material tangent modulus for the Knowles model is (Knowles, 1977; Suchocki, 2011):
311
b i o s y s t e m s e n g i n e e r i n g 1 8 6 ( 2 0 1 9 ) 3 0 7 e3 2 2
Ce ¼ Ceijkl
n1
2 1 mknowles Bknowles 1 dik Bjl þ djl Bik þ dil Bjk þ djk Bil þ 1þ ðI1 3Þ I1 dij dkl Bij dkl dij Bkl 2 3 3 J n n2
12 m Bknowles ðn 1Þ Bknowles 1 1þ ðI1 3Þ Bij Bkl I1 Bij dkl þ dij Bkl þ I1 dij dkl þ2 knowles n 3 3 J n ¼
(10)
2 þ knowles ð2J 1Þdij dkl D1
Damage and undamaged stress relationship can be described through the following (Schmidt & Kaliske, 2004): DEVSd ¼ ð1 dÞDEVSe
(11)
where d is a damage variable ranging from value 0 which indicates undamaged to 1 for fully damaged material. Note that the damage is applicable to the isochoric part only. The damage variable is expressed as: ex
ex
d ¼ da ðaÞ þ db ðbÞ
(12)
ex da ðaÞ
ex db ðbÞ
and represent un-continuous and continwhere uous damage variables, respectively. The un-continuous damage function: ex ∞ da ðaÞ ¼ da
a 1 exp ha
(13)
where the un-continuous damage parameter, aðtÞis described as: aðtÞ ¼ max WðsÞ
(14)
s2½0;t
which can be expressed numerically as: anþ1 ¼
8 <
Wnþ1 ; Wnþ1 > an ; otherwise : an
(15)
On the other hand, the continuous damage part in Equation (12) is as follows: b ex ∞ db ðbÞ ¼ db 1 exp hb
where Ceiso is the isotropic material tangent modulus, which in this case is given by Equation (10). In addition, the derivative in Equation (19) is given through the chain rule:
(16)
vd vda va vdb vb ¼ þ vW va vW vb vW where the partial derivatives are given as: ∞ vda da a ¼ exp ha va ha ∞ vdb db b ¼ exp hb vb hb 8 va < 1 ¼ vW : 0
; anþ1 ¼ Wnþ1 > an ¼ Wn ; otherwise
vb jWnþ1 Wn j ¼ ¼ signðWnþ1 Wn Þ vW Wnþ1 Wn
Zt vWðsÞ ds bðtÞ ¼ vs
d DEV Svd nþ1 ¼ DEV S0;nþ1 þ
Equation (17) can be expressed in numerical form as:
Hvd j;nþ1
(22)
where the term Hvd is evaluated through the convolution j integral:
0
(17)
N X j¼1
Ztnþ1
0
(21)
where the sign in Equation (21) is used to provide an absolute value of the function evaluated. Deviatoric damaged stress in Equation (11) is used in the finite viscoelastic deviatoric stress (Kaliske & Rothert, 1997), which is expressed as:
Hvd j;nþ1 ¼ gj
where:
(20)
tnþ1 s d DEV Sd ðsÞ ds xj ds
1 exp Dt h i xj Dt zexp DEV Sd0;nþ1 DEV Sd0;n Hvd j;n þ gj Dt xj x j
(18)
(23)
Tangential modulus in terms of damage energy potential can be expressed as:
which is based on the linear rate equation of the internal stress,
bnþ1 ¼ bn þ jWnþ1 Wn j
Cdiso ¼ ð1 dÞCeiso
vd DEVSe 5DEVSe vW
(19)
_ j þ 1 ¼ gj DEV S_ 0 H xj
(24)
312
b i o s y s t e m s e n g i n e e r i n g 1 8 6 ( 2 0 1 9 ) 3 0 7 e3 2 2
5. Calculate damage function. d. 6. Obtain damage stress, DEVSd and damaged material tangent modulus, Cdiso . vd 7. Calculate Hvd j;n and DEVSnþ1 using damaged stress. vd d 8. Calculate Cnþ1 using Ciso . vd 9. Update svd nþ1 and Cnþ1 as STRESS and DDSDDE in UMAT subroutine. 10. Update variable Hvd j;nþ1 and bnþ1 from current increment.
It should be noted that the summation of relaxation constants, gj need to be equal to 1. Using the tangent modulus given in Equation (19), viscoelastic tangent modulus for the damaged potential is given by: 9 8 Dt > > > > 1 exp < N xj = X e Cdiso;nþ1 1 þ ¼ C þ g Cvd i nþ1 vol;nþ1 Dt > > > > xj j¼1 ; :
(25)
A push forward transformation of the second PiolaKirchhoff stress tensor from Equation (22) yields the Kirchhoff stresses in the current configuration: T vd DEV tvd nþ1 ¼ Fnþ1 DEV Snþ1 ðFnþ1 Þ
(26)
which then gives the total stress tensor through: 0 vd tvd nþ1 ¼ Jnþ1 Unþ1 þ DEV tnþ1
(27)
and finally, the Cauchy (true) stress tensor: 1 vd svd nþ1 ¼ J tnþ1
(28)
The constitutive model described before was implemented in Abaqus (Abaqus, 2015) using User Material (UMAT) Fortran subroutine. This subroutine requires updated Cauchy stresses from Equation (28) (STRESS) and Jacobian matrix (DDSDDE) of material tangent modulus (Equation (25)), respectively for each increment step. The steps under the UMAT subroutine are provided as follows:
2.3.
Parallel rheological framework (PRF) model
An alternative model to those presented in the previous section (Prony series with damage) is the Parallel Rheological Framework (PRF) model available in the material library of Abaqus (Abaqus, 2015). The model is represented in Fig. 3, which is almost similar to the Prony series model, except that the first network on the left of the figure now consists of an elastoplastic network, as well as the power law viscous dampers response. The elastic response in each network is given by a hyperelastic strain energy potential, which is similar to those used in the Prony series model in the previous section. It is assumed that all networks have the same functional form of the strain energy potential and that the total strain energy, WT is equal to the weighted sum of the strain energies of all the networks: WT ¼
N X
si W Cei
(29)
i¼0
1. Get input from deformation gradient at the beginning of increment, Fn (DFGRD0) and deformation gradient at the end of the increment, Fnþ1 (DFRGRD1). Obtain state variables Hvd j;n and bn from previous increment. 2. Calculate volume preserving left Cauchy strain T tensor B ¼ FF and J ¼ detðFÞ. 3. Calculate neo-Hookean strain energy, Wn and Wnþ1. Calculate corresponding neo-Hookean Cauchy stress, sen and sen . 4. Calculate hyperelastic tangent modulus, Ce .
where si denotes the ith volume fraction of the network. The strain energy used is the Ogden hyperelastic model (Ogden, 1972, 1997): Wogden ¼
2mogden aogden 1 a a l1 þ l2ogden þ l3ogden 3 þ ogden ðJ 1Þ2 2 aogden D1
(30)
where li are the deviatoric principal stretches, which is 1 related to principal stretch through li ¼ J3 li , mogden is the ogden is the shear modulus, aogden is a material constant and D1 bulk modulus. Note that the Ogden model is used here
Fig. 3 e Nonlinear viscoelastic model (Parallel Rheological Framework) illustration.
b i o s y s t e m s e n g i n e e r i n g 1 8 6 ( 2 0 1 9 ) 3 0 7 e3 2 2
(Ogden, 1972, 1997) instead of Knowles model in the previous section due to the latter model is not available in Abaqus material library. However, it should be noted that both models are capable to simulate large deformation materials. The stiffness ratios in Equation (29) must satisfy the following relation: N X
si ¼ 1
(31)
i¼0
In addition, the remaining viscous networks can be modelled as power law, hyperbolic sine law or BergstromBoyce law. In this work, the strain hardening power law function is used due to its simplicity, which can be represented as:
m 1=ðmþ1Þ cr ε_ ¼ Aq~nn ðmv þ 1Þεcr
(32)
cr where ε_ is the equivalent creep strain rate, εcr is the equiva~ is the equivalent deviatoric Kirchoff stress, lent creep strain, q and A, mv and nn are material parameters. The range of parameter, mv must be: 1 < mv 0 for the modelling in Abaqus. The plasticity response of the PRF model (the first network in Fig. 3) is defined using equivalent stress (plastic stress), sp and equivalent plastic strain, εp in Abaqus. The software allows a direct tabular entry of the plastic stress and plastic strain obtained from stressestrain experimental cyclic tests data. Finally, the stress softening function is defined using the Ogden-Roxburgh model (Ogden & Roxburgh, 1999), using a damage variable, h , which varies with deformation:
m ~m 1 Udev U dev h ¼ 1 erf r mm þ bm Um dev
(33)
~ is the deviatoric part of the strain energy density of where U dev the finite strain viscoelastic model, and Um dev is the maximum ~ m at a material point during its deformation history. value of U dev The other parameters, r, mm and bm are material parameters, and erfðxÞ is an error function. The model can be implemented m
313
in Abaqus through modification of the input file (.inp) as shown in Fig. 4.
3.
Results
3.1.
Experimental results
The mean diameter of fresh fibres is 0.22 ± 0.02 mm, and the mean diameter for dried fibres is 0.19 ± 0.02 mm, while the moisture content for fresh OPMF was 20.04 ± 1.41% (w/w). Observation of oil on the surface of fresh OPMF (transverse cross-section of oil palm fruitlet) is shown in Fig. 5. Oil droplets (yellow in colour) are mostly observed in the region of parenchyma cells, specifically on the outside of the vascular bundle of OPMF (sclerenchyma cells) as shown in Figure 5(c),(d). During processing involving mechanical load and heat, most of the oil is removed from the fruitlet. It is worth noting that Scanning Electron Microscopy of dried mesocarp fibres was reported by Hanipah et al. (2017), whose focused on the contribution of silica bodies towards integrity of the fibres. For the stress relaxation and cyclic tests, normalised stress and strain were used for comparison from both fresh and dried fibres test results, as shown in Fig. 6. Normally stress relaxation and cyclic test are presented as stress and strain (directly) as shown by Omar et al. (2016) and Hanipah et al. (2017). However the problem with this is that it is difficult to observe the difference between the relaxation curves at different deformation (also similar to cyclic test for elastic and plastic strains) which is why the normalised values are presented here, so that this allows the normalised stress to overlap each other for results at different deformations in order to obtain a master curve. By referring to Fig. 6, the following normalised stress is used: ¼ snormalised t
st¼n st¼0
Fig. 4 e Material model script included into the Abaqus input file for PRF model.
(34)
314
b i o s y s t e m s e n g i n e e r i n g 1 8 6 ( 2 0 1 9 ) 3 0 7 e3 2 2
Fig. 5 e Observation of oil on the surface of fresh OPMF. a) and b) longitudinal cross-section; c) and d) transverse crosssection; e) and f) fractured fibre and tensile tests.
Fig. 6 e Illustration for (a) stress relaxation and (b) permanent set illustration.
where st¼0 represent maximum stress before relaxation occur, and st¼n is the stress at any time between t ¼ 0 s to t ¼ 100 s relaxation time shown in Fig. 6(a). This allows the stress to be valued between 0 and 1. Likewise, the normalised elastic strain is defined as follows: εne ¼
εe 100% εp
(35)
and the normalised plastic strain: εnp ¼
εp 100% εp
(36)
where εe and εp represent elastic and plastic true strain at zero stress, as shown in Fig. 6(b). The comparison results for stress relaxation are shown in Fig. 7 and Table 1. Master curves of stress relaxation results for both dried and fresh fibres are obtained, suggesting a viscoelastic model that is independent of deformation (viscoelastic behaviour is sensitive to moisture due to solideliquid behaviour causing time-dependent behaviour). However, it is observed that the stress relaxation of the fresh fibres reduced more than the dried ones after 1 s relaxation, which suggested that the fresh fibres contained more moisture than the dried
315
b i o s y s t e m s e n g i n e e r i n g 1 8 6 ( 2 0 1 9 ) 3 0 7 e3 2 2
Fig. 7 e (a) Comparison of relaxation stress results between dried and fresh fibres. Cyclic loading-unloading results for (b) fresh mesocarp fibres and (c) dries mesocarp fibres. (d) comparison of permanent set (from cyclic test) between dried and fresh fibres. Note that the true strain has dimension of mm mm¡1 (dimensionless).
Table 1 e Normalised stress values at different relaxation time for fresh and dried fibres. The parentheses refer to standard deviations. Material Fresh fibres
Dried fibres
Deformation (mm)
Normalised stress at 1s
1 2 3 1 2 3
0.773 0.778 0.784 0.817 0.798 0.823
Normalised stress at 10s
(0.003) (0.009) (0.009) (0.006) (0.011) (0.018)
0.628 0.632 0.636 0.699 0.685 0.713
fibres. This relaxed stress is related to the viscous stress in spring-dashpot components of the Prony series. Clearly, different viscoelastic model parameters are needed to model both the fresh and dried fibres. Since viscoelastic behaviour is sensitive to moisture (due to solideliquid behaviour causing time-dependent behaviour), this is the reason why such difference is noticed. To do this, a one dimensional viscoelastic model was used (see Goh, Charalambides, and Williams (2004) and Mohammed (2014) for details of model derivation), where the equation used is a follows:
(0.016) (0.009) (0.011) (0.004) (0.016) (0.017)
i
zi
0.502 0.497 0.493 0.594 0.590 0.628
(0.030) (0.006) (0.012) (0.0001) (0.021) (0.015)
where the true stress, so , is related to the nominal stress, Po , through so ðtn Þ ¼ Po ðtn Þlðtn Þ, and the term, Z t dPo ðsÞ hi ðtÞ ¼ gi exp ts ds ds. The one dimensional equaxi 0
tion (Equation (37)) is comparable to the three dimensional forms (Equations (22) and (23)) (Mohammed, 2012, 2014). The advantage of this model is that it can be implemented directly into a spreadsheet with calibration conducted using the least squares method, as demonstrated by Goh et al. (2004). This can reduce the time needed to find the best model parameters
0 1 △t 1 exp N B X C zi Bexp △t hi ðtn Þ þ gi sðtnþ1 Þ ¼ g∞ so ðtnþ1 Þ þ lðtnþ1 Þ ½Po ðtnþ1 Þ Po ðtn ÞC @ A z △t i¼1
Normalised stress at 100s
(37)
316
b i o s y s t e m s e n g i n e e r i n g 1 8 6 ( 2 0 1 9 ) 3 0 7 e3 2 2
to be used in a more complex 3D viscoelastic model. Figure 6 shows the results of fitting the tests with the Prony series relaxation constants (first and second rows in Table 2) by making the stress in Equation (37) to become dimensionless through Equation (34), where the difference of parameters between the fresh and dried fibres are the short and long term constants (g1 and g∞ , respectively). In particular, the difference of g∞ values shows significant changes of the residual stress for both dried and fresh fibres. On the other hand, it is worth noting that Xiang et al. (2015) reported changes of Prony series constants (stress relaxation tests) from oil palm fibres pre-treated with alkali under different concentrations, which is due to changes of the lignocellulosic contents within the fibres. Both results in this work and Xiang et al. (2015) suggested that moisture and lignocellulosic contents are responsible for the viscoelastic behaviour of oil palm fibres. The cyclic test results for dried and fresh fibres are shown in Figure 7(b),(c). Higher stressestrain curves are shown from the fresh fibres compared to the dried ones, which can be due to damages within the dried fibres after the milling, washing and drying processes. On the other hand, the results for permanent set (elastic and plastic strains) comparison are shown Fig. 7(d). Notice that while both dried and fresh fibres show increasing plastic strain at increasing deformation (more damage at larger deformation), there is no clear elastic and plastic strains pattern difference between the dried and fresh fibres at all deformation tests. This suggested that damage occurring in the fibres can be due to lignocellulosic interaction inside fibre rather than moisture (e.g. effect of oil content). Lignocellulosic interaction can cause fibre damage, where Adler and Buehler (2013) showed elasticeplastic behaviour due to slippage between cellulose-hemicellulose-lignin interface. Experimentally this behaviour has been reported by Keskes et al. (2003), who refer the phenomena as the “stick slip” behaviour at the interfaces. This is the ability of the cellulose-hemicellulose-lignin interface to recover part of its elasticity due to damage from the interface slippage. On the other hand, Xiang et al. (2015) found correlation between lignocellulosic content with stress relaxation and cyclic test. For example, for relationship between elastic and plastic strains with cellulose content of alkali pre-treated oil palm fibres (Xiang et al. (2015); Omar et al. (2016)), even though increases of cellulose content (23%e58%) was reported at higher concentration (0e40%), this does not translate to
stronger fibres. This can be due to the weakening of the interface of the cellulose, hemicellulose and lignin interface from the pre-treatment conditions. Note that Sreekala et al. (1997) reported the following composition of oil palm mesocarp fibres: Lignin: 11%, cellulose: 60%, ash: 3%. It should be noted that a few fibres failed before ~0.15 strain, where observation of the fractured area along the fibres length is shown in Figure 5(e),(f). It is suggested that two mechanisms of damage occurred in the fibre due to deformation; one is due to splitting of fibre at the cross section, and the other is due to small fractures along the fibre length leading to total fracture of the fibre. Experimental results on the anisotropic behaviour of oil palm fibres (dried fibres) are shown in Fig. 8. The results shown are for clockwise (CW), counter clockwise (CCW) and not twisted fibres (uniaxial) as control. All tests were conducted on the same batch of fibres to minimize the effect of stressestrain variation, as mentioned by Omar et al. (2014). The results show higher results for CW and CCW fibres compared to the uniaxial fibres after 0.03 strain. However, considering the standard deviation of the tests are high, it is quite difficult to draw conclusion from these results. It is worth noting that Sreekala et al. (1997) reported theoretical values of microfibril angle between 42 and 46 . From the previous x-ray microtomography results by Omar et al. (2016) and Omar, Hafid, Mohammed, Baharuddin, & Abdullah (2017), the fibres bundle in the microscale can be considered as a tranverse isotropic material, due to the cellular structure of the fibre which is continuous along the fibre length. However, it is possible that anisotropic behaviour is dominant at microfibril level, as suggested by Sreekala et al. (1997). This highlighted the complexity of oil palm fibres, which require further detailed analysis under different loading condition. An example would be testing of fibres under compression mode for more detailed analysis on anisotropic behaviour. The study by Ueda, Saito, Imahori, Kanazawa, and Jeong (2014) on carbon fibre can be referred, where in their case they set the height of a carbon fibre sample to be ~10 mm and diameter of ~5 mm (diameter-height ratio of 0.5) inside a SEM microcompression tester. This would minimize buckling of a fibrous sample. For our case here, we need to set an oil palm
Table 2 e Prony series viscoelastic model with a damage function parameters used. Model constant gj (for xj ¼ 0.1, 1, 10, 100, 1000 and ∞s)
mknowles (MPa) Bknowles n Dknowles (MPa) 1 d∞ a and ha (MPa) d∞ b and hb (MPa)
Value 0.47, 0.12, 0.06, 0.05, 0.03, 0.27 (dried fibres) 0.57, 0.12, 0.06, 0.05, 0.03, 0.17 (fresh fibres) 3000 0.5 0.2 1 109 5 and 600 1 and 10 000
Fig. 8 e Experimental results on twisted-tensile fibre samples to investigate possible anisotropic behaviour.
b i o s y s t e m s e n g i n e e r i n g 1 8 6 ( 2 0 1 9 ) 3 0 7 e3 2 2
fibre sample dimension at ~250 mm height and diameter of ~500 mm to obtain the same ratio as Ueda et al. (2014). SEM results of dried fibres was provided in the previous works (see Hanipah et al. (2017)), whereas no SEM was conducted the fresh fibres due to the existence of oil and moisture on the surface of the fibre. Considering that this paper focussed on viscoelastic behaviour, a comprehensive study on anisotropic behaviour of these lignocellulosic fibres and on the morphological characterisation of the surface of dried and fresh fibre need to be conducted separately in the future.
3.2.
Modelling results
3.2.1.
Prony series viscoelastic model with damage function
The finite viscoelastic model was fitted to cyclic tensile test results on dried fibres to demonstrate the capability of the model. The fresh fibres results were not included here due to similarity of cyclic tests results pattern when compared with the dried fibres (except that the fresh fibres are stiffer than the dried fibres). From the experimental results the true strainε, is calculated as: ε ¼ lnðl =lo Þ, whereas true stress, s, is s ¼ 4Fl=pD2 lo , where F, D, l and lo are force, fibre diameter, deformed length, and original length, respectively. A three-dimensional single element model (C3D8H) was loaded under cyclic uniaxial tension (at 1 mm/s). This is a 8 nodes element with 8 integration points. The modelling results shown in Fig. 9 are for (a) viscoelastic model with no damage function, (b) viscoelastic model with continuous and un-continuous damage function, (c) viscoelastic model with un-continuous damage function only, and (d) viscoelastic model with continuous damage function only. The following calibration method is used. Firstly, the model without damage was fitted to small deformation cyclic results (1 mm), where it was assumed that at this deformation (at 5% strain), damage
317
was minimal. The parameters used are shown in Table 2. The viscoelastic constants (first row in Table 2) were set to fit the relaxation curves in Fig. 7(a) (which made the values unique). The other Knowles model parameters (mknowles Bknowles , n and ) were set to fit the small deformation results, where Dknowles 1 the effect of the model parameters is shown later in Fig. 10. The damage function was then included, and the model was fitted further with larger deformation tests, at 2 mm (10% strain) and 3 mm (15% strain). One can see that the un-continuous damage influenced the primary loading curve more than the continuous damage function. Even though the model with damage in Fig. 9(b) fits the cyclic results at 1 mm deformation, it only agrees to the initial loading curves of the larger deformations (2 and 3 mm deformations). The unloading-reloading curves do not agree with the model, which can be due to plasticity, as well as large stress softening of the unloading-reloading curves. To further investigate the influence of the damage function, a parametric study was conducted by separating the Knowles model parameters, as well as the un-continuous and continuous damage functions at large deformation (3 mm). Figure 10(a)e(d) show results with different Knowles model parameters without including any viscoelastic and damage effect. The Knowles model parameters in Table 2 were used initially before the other model parameters were varied. Increasing the value of Bknowles and n caused lower stressestrain results, whereas increasing the value of mknowles increase the stressestrain results. On the other hand, (related to compressibility of the increasing the value of Dknowles 1 material) does not cause any significant changes to the revalues are shown to improve the sults, and larger Dknowles 1 convergence of the modelling results. Figure 10(e),(f) shows results of different un-continuous damage parameters, da and ha , respectively, whereas Fig. 10(g),(h) for continuous damage
Fig. 9 e (a)e(d) Cyclic modelling results using the Prony series viscoelastic model with a damage function.
318
b i o s y s t e m s e n g i n e e r i n g 1 8 6 ( 2 0 1 9 ) 3 0 7 e3 2 2
Fig. 10 e Sensitivity analysis for Prony series viscoelastic model with a damage function: (a) to (d) Knowles model parameters; (e) and (f) un-continuous damage function; (g) and (h) continuous damage function.
parameters, db and hb , respectively. The Knowles model parameters, viscoelastic constants and damage function parameters values in Table 2 were initially set before the damage function parameters were varied. One can see that the uncontinuous damage function influenced the primary loading curve, whereas the continuous damage function affects the unloading-reloading curves more than the primary loading curve, similar to the finding in Figure 9(c),(d). In all cases, the damage function is not capable of simulating the unloadingreloading experimental curves under 2 mm (10% strain) and 3 mm (15% strain) in Fig. 7(b). To improve modelling under
these larger deformations, PRF model is used, as shown in the next section.
3.2.2.
PRF model
A similar model geometry and element (C3D8H) used in the previous modelling section was applied to the PRF model. The modelling results are shown in Fig. 11 for large deformation (3 mm) test. The model parameters used are shown in Table 3, which were set as follows. The shear modulus was set to be the same as the one used in previous section (3 GPa). The viscoelastic network constants were also set to be the same as
b i o s y s t e m s e n g i n e e r i n g 1 8 6 ( 2 0 1 9 ) 3 0 7 e3 2 2
319
Fig. 11 e (a) Comparison between PRF model by removing plasticity and stress softening functions; (b) to (d) the effect of different power law viscoelastic constants; (e) to (g) the effect of different stress softening (Ogden-Roxburgh) parameters.
the Prony series constants used in the previous model. The plastic stresses and plastic strains were estimated from the permanent set strain values and corresponding stress values from the experimental results in Fig. 7(b). Figure 11(a) shows a significant improvement of the model which is now closer to the experimental results in Fig. 7(b) when both stress softening and plasticity function are used. This can be clearly seen when the plasticity function is removed from the model, which is consistent to the modelling
results using the Prony series viscoelastic model only in Fig. 8. Furthermore, the stress softening function by Ogden and Roxburgh (1999) only influenced the unloading-reloading part of the cyclic tests. Figure 11(b)e(d) shows the effect of the power law viscoelastic constants, where the values in Table 3 was initially used without the stress softening and plasticity function included. Increasing the value of A reduced the stressestrain curves of the cyclic results. In contrast, reducing the constant nn slightly reduces the cyclic results. It
320
b i o s y s t e m s e n g i n e e r i n g 1 8 6 ( 2 0 1 9 ) 3 0 7 e3 2 2
Table 3 e PRF model parameters used to simulate cyclic tests under large deformation. Parameter mogden (MPa) aogden ogden D1 si (ratio i ¼ 1 to 6) A mv nn r mm bm sp (MPa) εp
Values 3000 1 1 109 0.47, 0.12, 0.06, 0.05, 0.03, 0.27 0.003 1 0.5 2 2 0.3 100, 180, 200, 220 0, 0.02, 0.05, 0.1
can be seen that when mv is 0.1 and 10, the cyclic unloadingreloading is significantly changed compared to the other two constants. A parametric study was then conducted using different stress softening (Ogden & Roxburgh, 1999) constants. The hyperelastic, viscoelastic, stress softening and plasticity function values in Table 3 were set initially, before any one of the stress softening constant was varied. The results are shown in Figure 11(e)e(g), which shows that increasing the values of r, mm and b shifted the unloading-reloading curves further from the zero stress. This also confirmed that the stress softening function by Ogden and Roxburgh (1999) only influenced the unloading-reloading part of the cyclic tests, as mentioned before in Fig. 11(a).
4.
Discussion
The results in previous section highlighted the complex behaviour of oil palm fibres under large deformation (larger than 5% strain). It was mentioned earlier that the modelling work on oil palm fibres can be used in biodegradation and biocomposites applications. The Prony series model with damage is valid for small deformation, which can be useful for biodegradation application, since it does not involve large deformation of the fibres. In addition, it is worth mentioning the cyclic test result of an isolated sclerenchyma tissue (Aristolochia macrophylla) by Kohler and Spatz (2002) is shown in Fig. 12(a). Comparison with our model at 1 mm deformation
in Fig. 12(b) shows quite similar hysteresis (from unloadingreloading) behaviour. To relate the biodegradation study (for future work) with the model proposed here, we first show a microscopy image in Fig. 13(a) of a network of microbial (fungi) filaments and spores on a fibre surface. The microbes, in this case fungi, releases a group of enzymes (cellulase, hemicellulose and lignolytic enzymes) which then degrade the lignocellulosic content within the fibre. This degradation will then cause volume reduction of the fibre, where the volume changes can be compared between the fresh and degraded fibre. An X-ray micro-tomography (mCT) analysis can be conducted by scanning a fresh fibre prior to experiment and another scan after degradation test, as shown in the illustration in Fig. 13(b). An example of a mCT scan of an oil palm fibre can be seen in Fig. 13(b), which shows cellular structure within the fibre cross section. Omar, Hafid, Mohammed, Baharuddin, & Abdullah (2017) recently reported m-CT results that revealed volume reduction after degradation process of oil palm fibres due to delignification process. The modelling aspect of natural fibre degradation is not well developed compared to experimental studies. Microbial activities that cause degradation and volume reduction can be modelled using FE micromechanics geometry developed from measured data such as mCT scan (as demonstrated by Omar et al. (2016)). The damage function proposed here then need to be modified such that it will describe degradation over time after exposure to microbial agent on the surface of the fibre. The idea of using a damage function to describe degradation has been reported by previous researchers in biomedical stents (Muliana & Rajagopal, 2012) and biodegradable plastic (Viera, Guedes, & Tita, 2014) applications. On the other hand, for biocomposites application, the PRF model can be used considering the good agreement of the model with the experimental results under large deformation (3 mm cyclic tests). However, it should be noted that the surface of the fibres consists of silica bodies (observed through SEM and FTIR analysis by Omar et al., 2016), where the results from previous modelling studies by Omar et al. (2014) and Omar et al. (2016)) showed that silica bodies (which are partly embedded on the outer layer of oil palm fibres) did not contribute to the strength of the fibre as a whole. However, in biocomposites application, Hanipah et al. (2016) showed that the silica bodies contributed to the mechanical behaviour of fibre-matrix interface by preventing sliding of the fibres
Fig. 12 e (a) Cyclic test results of a plant tissue; (b) modelling results of oil palm fibres without damage function.
b i o s y s t e m s e n g i n e e r i n g 1 8 6 ( 2 0 1 9 ) 3 0 7 e3 2 2
321
Fig. 13 e (a) SEM image of a network of filaments of fungi with spores on the surface of an oil palm fibre, (b) geometric representation of mCT scan on fibre before and after degradation for fibre volume reduction analysis.
during deformation of biocomposites. This phenomenon was previously reported by Nascimento, Ferreira, Monteiro, Aquino, and Kestur (2012) on piassava fibres for their biocomposites application. Both works on oil palm fibres biodegradation and biocomposites are currently underway. Finally, it is hoped that the understanding of the mechanical behaviour of these fibres in this paper can be used in higher value-added applications, such as filler in 3D printer's filament materials and material in biodegradable biomedical device.
5.
Conclusion
Viscoelastic behaviour and damage occurring in oil palm fibres were reported from stress relaxation and cyclic tests. A Prony series finite viscoelastic model was then plotted to the experimental results, which agreed only with small deformation test. For large deformation test, parallel rheological framework (PRF) model that consists of nonlinear viscoelasticity, plasticity and stress softening is needed. It is therefore proposed that the Prony series viscoelastic model with damage function to be used for biodegradation, whereas the PRF model is suitable for biocomposites application involving large deformation.
Acknowledgement This work was supported by Kyushu Institute of Technology, Japan and Universiti Putra Malaysia, and Putra Berimpak
Grant Scheme 2017 (Matching Grant 9300433, UPM/700-2/1/ GPB/2017/9521400 and UPM/700-2/1/GPB/2017/9416700).
references
Abaqus. (2015). User's manual version 6.14. Providence, USA: Karl, Hibbit and Sorensen. Adler, D. C., & Buehler, M. J. (2013). Mesoscale mechanics of wood cell walls under axial strain. Soft Matter, 9(29), 7138e7144. Baskin, T. I. (2005). Anisotropic expansion of the plant cell wall. Annual Review of Cell and Developmental Biology, 21(1), 203e222. Baskin, T. I., & Jensen, O. E. (2013). On the role of stress anisotropy in the growth of stems. Journal of Experimental Botany, 64(15), 4697e4707. € nsson, H. (2016). Shifting foundations: The Braybrook, S. A., & Jo mechanical cell wall and development. Current Opinion in Plant Biology, 29, 115e120. Goh, S. M., Charalambides, M. N., & Williams, J. G. (2004). Determination of the constitutive constants of non-linear viscoelastic materials. Mechanics of Time-dependent Materials, 8(3), 255e268. Hanipah, S. H., Mohammed, M. A. P., & Baharuddin, A. S. (2016). Nonlinear mechanical behaviour and bio-composite modelling of oil palm mesocarp fibers. Composite Interfaces, 23(1), 37e49. Hanipah, S. H., Xiang, L. Y., Mohammed, M. A. P., & Baharuddin, A. S. (2017). Study of non-linear mechanical behaviour of oil palm mesocarp fibers. Journal of Natural Fibers, 14(2), 153e165. Hurtado, J. A., Lapczy, I., & Govindarajan, S. M. (2013). Parallel rheological framework to model non-linear viscoelasticity, permanent set, and Mullins effect in elastomers. In Constitutive models for rubber VIII (pp. 95e100). London: Taylor & Francis Group.
322
b i o s y s t e m s e n g i n e e r i n g 1 8 6 ( 2 0 1 9 ) 3 0 7 e3 2 2
Kaliske, M., Nasdala, L., & Rothert, H. (2001). On damage modelling for elastic and viscoelastic materials at large strain. Computers & Structures, 79, 2133e2141. Kaliske, M., & Rothert, H. (1997). Formulation and implementation of three-dimensional viscoelasticity at small and finite strains. Computational Mechanics, 19, 228e239. Keskes, J., Burgert, I., Fruhmann, K., Muller, M., Koll, K., Hamilton, M., et al. (2003). Cell-wall recovery after irreversible deformation of wood. Nature Materials, 2, 810e814. Khalil, H. P. S. A., Ismail, H., Rozman, H. D., & Ahmad, M. N. (2001). The effect of acetylation on interfacial shear strength between plant fibres and various matrices. European Polymer Journal, 37(5), 1037e1045. Knowles, J. K. (1977). The finite anti-plane shear field near the tip of a crack for a class of incompressible elastic solids. International Journal of Fracture, 13, 611e639. Kohler, L., & Spatz, H. C. (2002). Micromechanics of plant tissues beyond the linear-elastic range. Planta, 215, 33e40. Mohammed, M. A. P. (2012). Mechanical characterisation, processing and microstructure of wheat flour dough. Phd Thesis. College London, UK: Imperial. Mohammed, M. A. P. (2014). Visco-hyperelastic model for soft rubber-like materials. Sains Malaysiana, 43(3), 451e457. Mohammed, M. A. P., Tarleton, E., Charalambides, M. N., & Williams, J. G. (2013). Mechanical characterization and micromechanical modelling of bread dough. Journal of Rheology, 57(1), 249e272. Muliana, A., & Rajagopal, K. R. (2012). Modeling the response of nonlinear viscoelastic biodegradable polymeric stents. International Journal of Solids and Structures, 49, 989e1000. Nascimento, D. C. O., Ferreira, A. S., Monteiro, S. N., Aquino, R. C. M. P., & Kestur, S. G. (2012). Studies on the characterization of pissava fibers and their epoxy composites. Composites Part A Applied Science and Manufacturing, 43(3), 353e362. Ogden, R. W. (1972). Large deformation isotropic elasticity-on the correlation of theory and experiment for incompressible rubberlike solids. Proceedings of the Royal Society of London, 326, 565e584. Ogden, R. W. (1997). Non-linear elastic deformations (1st ed.). New York: Dover Publications (Chapter 4). Ogden, R. W., & Roxburgh, D. G. (1999). Pseudo-elastic model for the Mullins effect in filled rubber. Proceedings of Royal Society A, 455, 2861e2877. Omar, F. N., Hafid, H. S., Baharuddin, A. S., Mohammed, M. A. P., & Abdullah, J. (2017). Oil palm fiber biodegradation: Physico-
chemical and structural relationships. Planta, 246(3), 567e577. Omar, F. N., Hanipah, S. H., Xiang, L. Y., Mohammed, M. A. P., Baharuddin, A. S., & Abdullah, J. (2016). Micromechanical modelling of oil palm empty fruit bunch fibres containing silica bodies. Journal of the Mechanical Behavior of Biomedical Materials, 62, 106e118. Omar, F. N., Mohammed, M. A. P., & Baharuddin, A. S. (2014). Effect of silica bodies on the mechanical behaviour of oil palm empty fruit bunch fibres. Bioresources, 9(4), 7041e7058. Placet, V., Cisse, O., & Boubakar, M. L. (2014). Nonlinear tensile behaviour of elementary hemp fibres. Part I: Investigation of the possible origins using repeated progressive loading with in situ microscopic observations. Composites Part A Applied Science Manufacturing, 56(1), 319e327. Schmidt, J., & Kaliske, M. (2004). Introduction of a viscoelastic damage approach into ANSYS. In 22nd CAD-FEM Users’ Meeting 2004, Dresden, Germany. Sreekala, M. S., Kumaran, M. G., Joseph, R., & Thomas, S. (2001b). Stress-relaxation behaviour in composites based on short oilpalm fibres and phenol formaldehyde resin. Composites Science and Technology, 61(9), 1175e1188. Sreekala, M. S., Kumaran, M. G., & Thomas, S. (1997). Oil palm fibers: Morphology, chemical composition, surface modification, and mechanical properties. Journal of Applied Polymer Science, 66, 821e835. Sreekala, M. S., Kumaran, M. G., & Thomas, S. (2001a). Stress relaxation behaviour in oil palm fiber. Materials Letters, 50(4), 263e273. Suchocki, C. (2011). A finite element implementation of Knowles stored-energy function: Theory, coding and applications. Archive of Mechanical Engineering, 3, 319e346. Ueda, M., Saito, W., Imahori, R., Kanazawa, D., & Jeong, T. K. (2014). Longitudinal direct compression test of a single carbon fiber in a scanning electron microscope. Composites Part A Applied Science Manufacturing, 67(12), 96e101. Vieira, A. C., Guedes, R. M., & Tita, V. (2014). Constitutive modeling of biodegradable polymers: Hydrolytic degradation and time-dependent behaviour. International Journal of Solids and Structures, 51, 1164e1174. Xiang, L. Y., Hanipah, S. H., Mohammed, M. A. P., Baharuddin, A. S., & Lazim, A. M. (2015). Microstructural, mechanical, and physicochemical behaviours of alkali pre-treated oil palm stalk fibers. BioResources, 10(2), 2783e2796.