A nonlinear finite viscoelastic formulation relative to the viscous intermediate configuration applied to plants

In the contribution at hand, a new formulation for finite strain viscosity relative to the viscous intermediate configuration is presented. The evolution of the viscous deformations is based upon a new numerical approach, which allows for a consistent consideration of anisotropic finite strain viscoelasticity, according to the authors knowledge. A standard Maxwell model is used to describe viscous behaviour at finite deformations. Furthermore, the orthotropic Yeoh material model is extended to include a distinction between behaviour under tensile and compression loading. The proposed formulation is validated, and parameters of the model are identified by material tests on Sorghum bicolor plants. Subsequently, numerical examples are shown to demonstrate the capabilities of the model. In general, the proposed Yeoh material formulation is shown to accurately represent the inability of fibres to carry compression loading. Furthermore, the viscoelastic approach, developed relative to the viscous intermediate configuration, is demonstrated to be capable of producing plausible results. Additionally, the mechanical behaviour of Sorghum bicolor plants is simulated using the introduced formulation. The results show that the contribution at hand describes a novel methodology to simulate the viscoelastic behaviour of plant materials reliably.

fibres for instance.Furthermore, the mechanical behaviour of plants depends on the amount of water included in the cells.If drying becomes significant, the material characteristics are subjected to change.In the contribution at hand, a first modelling approach, which addresses viscoelasticity and the difference in tensile and compression loading is presented.
Orthotropic Yeoh material models are general elastic material formulations at finite strains.Originally, it was developed for the isotropic case, see Reference 1. Subsequently, extensions have been introduced to enable the description of anisotropic behaviour, see References 2-4.Nevertheless, an incorporation of a distinction between tensile and compression within the Yeoh model is a novelty to the authors knowledge.
Plants often exhibit deformations at finite strains.Therefore, a material description in the finite strain framework is utilised.Furthermore, the material response is dominated by viscoelasticity.Thus, the aforementioned orthotropic Yeoh material model is extended to allow for the simulation of viscoelasticity at finite deformations.For modelling these characteristics, multiple approaches exist in the literature.One commonly used scheme is based upon the formulation of the viscous evolution equations in the eigenspace, see References 5-9.In Reference 5, it is noted that the assumption that the inelastic right Cauchy-Green strain C i and the elastic left Cauchy-Green strain b e being in the same eigenspace does not hold true in general.This aspect is also shown in Appendix D, where a simplified version of the model described in Reference 5 is used to compare the eigenvectors of the tensors b e and C i .Other approaches are also established in the literature, see References 10-12.A more recent formulation is depicted in Reference 13.Here, the viscous part of the right Cauchy-Green tensor is evolved in time.In References 14 and 15, it is mentioned that this evolution equation is only valid for isotropic tensor functions.In the contribution at hand, an anisotropic material, which exhibits anisotropic viscous deformations is investigated.In order to work around the solution of the viscous evolution equations in the eigenspace, a more straightforward solution scheme in the viscous intermediate configuration is introduced.The concept of an intermediate configuration is based on the split of the deformation gradient.In a multiplicative split, all parts exhibit transformation properties.Therefore, the intermediate configuration is introduced, and the viscous evolution equations are solved in this reference frame.To the authors knowledge, the formulation of the time evolution equations with respect to the intermediate configuration and the corresponding time integration algorithm are novelties introduced in this manuscript.The multiplicative split of the deformation gradient into elastic and viscous parts is broadly available in the literature, see References 16-18 among others.The formulation relative to the viscous intermediate configuration put forth in this manuscript, resulting from the multiplicative split of the deformation gradient, has certain advantages.Firstly, the consistent and stable time integration of Reference 19 is adopted for the evolution equation.Also, the evolution of these viscous deformations is based upon the work conjugate stress measures, as shown in Section 3.4.The formulation is consistent, since multiplicative splits of the deformation gradient and work conjugated stresses and strains based upon the dissipation inequality are utilised.Additionally, anisotropic viscous behaviour can be modelled by the proposed formulation, see Section 5.2.2.This represents another novelty according to the authors' knowledge.
For validating the model, material tests are performed on plants.Generally, significant uncertainty is associated with the determination of material parameters of biological structures.Material tests are conducted on Sorghum bicolor plants specifically.This plant is known to show significant viscoelastic material behaviour. 20However, for the tests depicted in the contribution at hand, uncertainty is not significant.The proposed formulation is capable of capturing the behaviour of other plants as well.The data obtained from the material tests are utilised for validation of the proposed formulation.
To summarise, in the present contribution, material tests are conducted on plant material, see Section 2. Subsequently, the principles of continuum mechanics are utilised to describe finite strain viscoelasticity with respect to the viscous intermediate configuration, described in Section 3. The constitutive equations are specified in Section 4.Then, numerical examples are shown to demonstrate the capabilities of the proposed formulation in Section 5. Additionally, validation of the formulation is carried out based on material tests on Sorghum bicolor plants, and material parameters are identified.Finally, a conclusion and an outlook are given in Section 6.

MATERIALS AND METHODS
This section describes the plant material used for the validation tests, and the methodology on which the tests are carried out.

Plant material
For validation tests, Sorghum bicolor L. Moench, a cereal crop plant, is used.The plants are grown in the Botanical Garden of Technische Universität Dresden under field conditions from June to September 2022.Only fully grown and visually healthy plants are harvested and, then, tested on the same day.Photographs of the growth over time are given in Figure 1.F I G U R E 2 Experimental setup.

Biomechanical testing
Three-point bending relaxation and cyclic loading tests are performed using a universal testing machine, and corresponding software (Zwick/Roell Allround Line Z005, testXpert II V3.5, Zwick/Roell GmbH & Co. KG).Displacements are recorded from the traverse movement and resulting forces with a 5 kN load cell (Zwick/Roell xforce P, Zwick/Roell GmbH & Co. KG).Supporting width is set at 180 mm.Samples are cut from the whole length of the harvested culms to a length of 200 mm just before testing to minimise drying and, thus, prevent an influence of water loss during testing.
The samples are positioned in such a way on the culm that they could be kept node-free or with the nodes positioned symmetrically, and as close to the sample ends as possible.The diameter is measured in the sample centre using digital calipers (Precise PS 7215, Burg-Wächter).For relaxation tests, maximum displacement is set at 5 mm with a testing speed of 2 mm/min and a holding time of 600 s.For cyclic tests, maximum displacements are set at 8 mm with a testing speed of 10 mm/min for a total of 10 cycles.The lower reversal point is set to occur at 0 N of measured force to ensure that the loading pin is always in contact with the sample.A graphical depiction of the experimental setup is given in Figure 2. The results of the test are further depicted in Section 5.4.An individual plant during growth is shown in Figure 3.

CONTINUUM MECHANICAL FRAMEWORK
This section briefly describes the continuum mechanical principles of viscoelasticity.The formulation shown in this work is mainly based upon References 5,14,21-23.

Configuration and movement
The configuration of the solid body  under observation is a subset of the Euclidean space  ⊂ R 3 .At time t = 0, the configuration of the body is considered as the reference configuration  0 , and is assumed to be stress free.For any time t > 0, the configuration is referred to as the current configuration  t .The deformation mapping  t ∶ X → x describes the motion of the body under observation.Based upon this mapping , the deformation gradient is defined as the relation between the deformation gradient and the gradient of the displacement vector u and the identity matrix .Furthermore, the determinant of the deformation gradient J ∶= det(F) is introduced, which transforms an infinitesimal volume element from the reference to the current configuration.The starting point of modelling viscoelasticity at finite deformations is the multiplicative split of the deformation gradient.This separation is carried out as into a viscous part F v , and an elastic part F e .The main property of the deformation gradient, that is, the transformation between different configurations, is preserved by the individual parts.Therefore, a viscous intermediate configuration is introduced.Based upon References 14,22-26, a consistent evolution law for the viscous deformation can be formulated relative to the intermediate configuration.The deformation gradient can be further split multiplicatively into volume preserving F and volume changing F vol parts as where  is the second-order identity tensor.Moreover, in plants, different material behaviour for purely volumetric and purely isochoric deformations can be observed as shown in the experimental results in Section 5.4.

Viscous deformations
Subsequently, the viscous deformations are introduced based on an evolution law.The fundamental quantity utilised in the contribution at hand is the velocity gradient in the intermediate configuration An additive decomposition of this quantity can be carried out into an elastic part L e and a viscous part L v , reading The evolution of the viscous deformation is achieved by utilising the definition of the viscous part of the velocity gradient in the intermediate configuration given by Here, the derivation from References 14, 23, and 27 is followed.

Elastic deformations
The elastic deformations are considered to be stress inducing.The evaluation of the induced stress requires an adequate strain measure.The deformation gradient is not an admissible strain measure, since it contains rigid body rotations.Therefore, the right Cauchy-Green strain tensor is utilised, given by Additionally, the elastic part of the right Cauchy-Green strain tensor reads Also, this quantity is the metric tensor of the intermediate configuration, see References 14 and 28.The elastic part of the velocity gradient in the intermediate configuration is given by

Strain energy function, dissipation inequality and stress measures
The strain energy function  arises from the partial Legendre transformation of Gibbs equation, see Reference 29.It mainly characterises the elastic response of the material under observation.In the current configuration, the Cauchy stress tensor can be defined as where g refers to the metric tensor in the current configuration.The Kirchhoff stress tensor can then be stated as in relation to the Cauchy stress tensor.Additionally, the transformation of the Kirchhoff stress tensor to obtain the second Piola-Kirchhoff stresses is given by In Equation ( 2), the viscous intermediate configuration is introduced.Based on the findings outlined in References 14 and 27, the viscous deformations evolve relative to this configuration.Therefore, the dissipation inequality can be formulated relative to the viscous intermediate configuration.This inequality is based upon the derivative of the strain energy function, which is expressed as Consequently, the dissipation inequality relative to the intermediate configuration can be formulated as compare References 14 and 18.Also, in the above Equation (13), which describes the dissipation , the Mandel stress tensor is defined as work conjugate to the viscous deformations, reading The partial derivative of the strain energy function with respect to the right Cauchy-Green strain tensor in the intermediate configuration from Equations ( 13) and ( 14) can be expressed as the second Piola-Kirchhoff stress tensor in the intermediate configuration, reading Equations ( 7) and (15) are further described in References 14 and 18. Utilising Equations ( 7), (10), (11) and (15), the Mandel stress tensor in the intermediate configuration is then expressed as This proposed formulation relative to the viscous intermediate configuration offers important advantages.Firstly, a consistent time integration, following Reference 19, is utilised in Section 4.2.Furthermore, the viscous deformations are evolved based upon a work conjugate stress measure, leading to a consistent formulation, see References 14, 18, 22, and 30.Additionally, there are no assumptions made regarding the eigenspaces of strain measures, compare References 5-9.These assumptions are shown in Appendix D to not be valid for general cases.Subsequently, the viscous part of the deformation gradient does not need to be symmetric.Therefore, more general formulations for the specific choice of the evolution operator for the inelastic deformation can be utilised.Furthermore, the formulation relative to the viscous intermediate configuration allows for an anisotropic evolution of viscous deformations, compare Section 5.2.2.

CONSTITUTIVE FORMULATION TO DESCRIBE VISCOELASTICITY OF PLANTS
In this section, the constitutive formulations are specified, which are utilised to represent the viscoelastic behaviour of plants.

Strain energy function
In Section 5.4, it can be observed that the mechanical response of plants is subdivided into two main features.Firstly, a contribution independent of time can be identified, since the stress converges asymptotically to a value not equal to zero.Secondly, a time dependent part influences the stress response, which declines over time.Therefore, the strain energy function is split additively into two parts as In Equation (17), the elastic  e and the viscous  v part of the strain energy function are given.The resulting rheological model is depicted in Figure 4. To obtain a generalised Maxwell model, more viscous branches can be added.Therefore, the material behaviour can be specified for different time scales.In the contribution at hand, one branch is sufficient for the conducted material tests, see Section 5.4.In general, the elastic response is characterised by a combination of volumetric and isochoric deformations.Therefore, the elastic contribution to the strain energy is split as In Equation ( 18), the volumetric  e V and isochoric part of the strain energy  e are introduced.A common assumption is that no volumetric portion exists in the viscous deformation. 5-7,31Therefore, the viscous part takes the form The elastic volumetric response is obtained by utilising the energy density function following References 5-7, and 31.The isochoric parts of the elastic and viscous strain energy density functions are chosen based on invariants as In Equation ( 21),  □ can refer to the elastic isochoric  e or the viscous part  v of the strain energy density.The formulation is taken from Reference 3. The original formulation for the isotropic case is stated in Reference 1. Furthermore, the principles to extend this kind of free Helmholtz energy to orthotropic material behaviour is presented in References 2 and 32.In general, n terms are considered in each sum.The values of n = 3 for the isotropic case, and n = 6 for the case with fibre contributions, have been established in the literature, see Reference 3. In Equation ( 21), the invariants Īi of the right Cauchy-Green strain tensor and the shear modulus equivalents a i , b j , c k , d l , e m , f n and g o are introduced.Therefore, an orthotropic Yeoh material model is obtained for both elastic springs in the rheological model.This material model is chosen to yield a general formulation, which can represent different kinds of plants.The invariants and their derivatives are specified in Appendix A. Utilising Equations ( 9), ( 10) and ( 20), the volumetric elastic Kirchhoff stress tensor takes the form In a similar manner, the isochoric Kirchhoff stresses are given by Here, an isochoric stress measure is obtained, since the isochoric part of the free Helmholtz energy is considered.Based upon Equations ( 17)-( 19), the total Kirchhoff stresses are decomposed into To model the small fibres observed in plants, an additional modelling assumption is made.In reality, fibres in plants or other bio-materials are typically too small to carry any load under compression.Therefore, in Equations ( 21) and ( 23), the Macaulay brackets are included for the invariants associated to anisotropic material behaviour These are used to indicate that the terms are only active under tension.The modelling strategy based upon the Macaulay brackets is taken from Reference 33.The free Helmholtz energy proposed in Reference 33 is formulated for collagen fibres.In Sorghum bicolor, the fibres consist of lignified cell walls.Consequently, a more general energy formulation is chosen.Nevertheless, a similarity of these two materials is that the fibres only act under tensile loads.Due to utilisation of Macaulay brackets, the fibres can be deactivated at compression.Thus, the modified isochoric free Helmholtz energy function reading and the isochoric part of the Kirchhoff stresses reading can be obtained.The above definition of the stresses is employed in all elastic springs in the rheological model.The difference of this material formulation compared to a standard formulation, where the fibres are not deactivated during compression, is investigated in Section 5.3.For the description of viscoelasticity, the specification of the viscous deformations remains and is addressed in the next section.

Evolution of viscous deformations
The starting point of the analysis is the assumption that only the viscous stresses  v induce viscous deformations.Furthermore, the time evolution of the viscous deformations depends on the Mandel stresses in the intermediate configuration, since this stress measure is work conjugate to the viscous deformations.The viscous part of the Mandel stress tensor C e is transformed to the intermediate configuration by Equation ( 16).As a consequence, the viscous part of the velocity gradient is a function of the viscous Mandel stress tensor The above equation is chosen as the constitutive equation, according to Reference 5.In Equation ( 27), the effective creep rate ̇ and the directional tensor N are introduced.In plants, a stress dependent creep rate is observed.Therefore, an approach similar to that in Reference 7 is followed.The creep rate is given by In Equation ( 28), the viscosity  t and the dimensional consistency parameter Σ are introduced.The dimensional consistency parameter is needed for obtaining the correct units of the creep rate.Utilising Equations ( 5) and ( 27), the time evolution equation of the viscous deformations takes the form For the application within a numerical simulation, Equation ( 29) is solved via the algorithm proposed by Reference 19, reading In Equation (30), □ n+1 refers to quantities at the next time step and □ n refers to the values of the quantities at the last converged time step.Here, the viscous Mandel stress tensor is utilised.This quantity is obtained from Equation ( 26) by applying Equation (16).As a consequence, viscous behaviour is obtained, which is isotropic during compression and anisotropic under tension.The exponential mapping of the tensor is further elaborated in Reference 27.Equation ( 30) is non-linear, and needs to be solved by an iterative procedure.Therefore, a local Newton-Raphson iteration scheme is employed at each Gauss integration point within the domain of a finite element.It is based upon reformulating Equation ( 30) as a residual expression, reading Subsequently, the Newton-Raphson method is applied by performing the consistent linearisation of Equation ( 31) stated as In Equation (32), the increment of the viscous deformations ΔF v for each iteration step of the local Newton-Raphson iteration scheme is introduced.Rearranging the linearised equation leads to the conditional equation for the increment In Equation (33), the computation of the inverse of a non-symmetric fourth order tensor is required.This quantity is obtained by representing the tensor as 9 × 9 matrix and calculating its inverse, see Appendix B. Finally, the update of the viscous deformations is performed by Here, the index k represents the Newton-Raphson iteration step, see Appendix C.

Algorithmic tangent
For the solution of the boundary value problem, the numerical framework of the finite element method is considered.In order to solve the nonlinear set of equations, the Newton-Raphson method is applied.Therefore, the algorithmic tangent is required.It is obtained by applying the consistent linearisation In Equation (35), the elastic tangent C = 2  g is introduced.Additionally, the computation of the total derivative dF v dF is required.It should be noted that the deformation gradient can also generally be replaced by any other convenient deformation measure.Here, the deformation gradient F is chosen, since its linearisation with respect to the nodal displacements is straight forward. 22,34The conditional equation for the total derivative is obtained by employing the implicit function theorem 24,25,27 to the residual expression for the viscous deformations, reading Subsequently, Equation ( 36) is solved for the total derivative dF v dF , given by

NUMERICAL EXAMPLES
Subsequently, numerical examples are presented to demonstrate the capabilities of the proposed formulations.Additionally, a physical validation is performed with material tests on plant material.

Hysteresis
The first example investigates the hysteresis observed for the inelastic materials.In this case, a cube with side length 5 mm as a model problem is considered.Two sets of boundary conditions are applied.In the first simulation, the top and bottom side are fully clamped.Next, a simulation is performed where the top and bottom side are only restrained in x 1 -and x 3 -direction.To prevent rigid body motions, the nodes on the middle layer along x 2 -direction are fixed along x 2 -direction.The displacement load is applied at the upper surface along x 1 -direction.The material parameters utilised in both simulations are stated in Table 1.The simulated time is 80 s.The first load cycle (0 to 50 to −50 to 0% strain), is applied in the first half of the simulated time.In the second half of the simulation, the amplitude of the first cycle is doubled (0 to 100 to −100 to 0% strain).As a consequence, the strain rate is doubled for the second loading cycle.The material behaviour for this test is chosen to be isotropic.The results of the simulation are shown in Figure 5.The stress is defined as the force measured in the direction of the applied load divided by the area of one side of the cube in the reference configuration.As it can be seen, the aforementioned increase in strain rate results in stiffer material behaviour for both sets of boundary conditions.Since the load is applied in a shorter time, the viscous deformations are smaller compared to the purely elastic ones.Therefore, the material behaves stiffer.Furthermore, the clamping leads to higher norms of stresses.Nevertheless, the two curves coincide for zero stresses, as expected.The convergence of the simulation with the released second degree of freedom is shown in Figure 6.As it can be seen, quadratic convergence is obtained indicating a proper implementation of the proposed formulation.
To highlight the influence of the viscosity  t on hysteresis, a second set of material parameters is investigated.The boundary conditions are similar to the previous example which is free in x 2 -direction.The material parameters of Table 1 are taken, but  t is increased to 150 s.The results of the simulation are given in Figure 7.For higher viscosity, the material  behaviour is more non-linear.Also, the difference between the loading and unloading portions of the characteristics is significantly more pronounced, indicating higher dissipation.

Stress relaxation
In this section, material parameter studies regarding the relaxation behaviour are conducted.First, isotropic material behaviour is studied.Subsequently, fibres, which lead to anisotropy, are considered.

Isotropy
For the isotropic material parameter study, a cube with edges of 5 mm is investigated.

TA B L E 2 Material parameters for the comparison of different values for a v
1 and  t .

Material parameter Value Unit
500 MPa common material parameters, which are used in all the simulations, are given in Table 2.The other material parameters, which are varied between different runs, are given in the resulting diagrams.Firstly, the overall material stiffness is investigated for different values of a v 1 , and the results are displayed in Figure 9. Here, the stress is calculated as the force obtained along the prescribed displacement divided by the initial area of one side of the cube.As it can be seen, the initial slope of the curve increases for increasing values of a v 1 .Additionally, the viscous deformations evolve faster.This can be observed by the steeper decrease of the slope of the curves during relaxation for an increasing value of the material parameter a v 1 .This is directly a result of Equation ( 28).The formulation depends on the stresses in the viscous intermediate configuration, and the behaviour is as expected.Nevertheless, the stress obtained for fully developed viscous deformation is not changed.Secondly, an investigation of the effect of the parameter  t is conducted.The results are displayed in Figure 10.Here, the stress is defined in the same way as in Figure 9.As it can be seen, an increase of  t does not change the initial slope of the stress-time-curve.However, the point in time, where viscous deformations become apparent, increases as expected.Furthermore, the stress obtained at any given moment for constant displacements differs.For higher values of  t , the developed viscous deformations evolve much slower.As a result, the stresses in the elastic spring of the viscous branch of the rheological model are larger.Furthermore, the stress obtained at infinite time is constant for changing  t .This can be seen by observing the curves for  t = 5 s and  t = 20 s, respectively.For larger values of  t , the time to reach that stress state is not within the simulated time.The convergence behaviour of a simulation is shown in Figure 11.As it can be seen, the convergence rate is quadratic, apart from errors due to numerical precision in the last iteration step.

Time (s) Stress (MPa)
F I G U R E 10 Comparison for different viscosities  t .

Anisotropy
To investigate the influence of fibres on the stress relaxation behaviour, a cube is simulated.The length of the edges is 5 mm.The boundary conditions are fully clamped opposing sides of the cube (Figure 12).The simulated time is 4 s.The maximum applied displacement results in 50% strain.The material parameters for all three different loading conditions are stated in Table 3.The material parameters are chosen in such a way that the fibres induce elastic and viscous stresses.
To study the influence of the loads, both along and perpendicular to the fibre direction, the material behaviour is chosen to be transversely isotropic.The material parameters are chosen in such a way, that fibres aligned with vector B are not considered.The bulk matrix is described by viscous material parameters.The displacement is applied as a linear ramp until 2 s.Then, it is kept constant.The resulting loading and relaxation behaviour is investigated.In the contribution at hand, anisotropy is induced by the fibres.The results are depicted in Figure 13.The stress is calculated for loading along x 3 and loading along x 2 by dividing the force measured along the prescribed displacements by the initial area of one side of the cube.For a combined loading along x 2 and x 3 the forces of both loaded directions are measured.Subsequently, the Euclidean is calculated and divided by the same area, resulting in the stress.The loading along x 3 results in a displacement of the fixed surface only along x 3 -direction.The fibre direction A is aligned with this coordinate direction.Therefore, it can be observed that the fibres influence the stress response and relaxation behaviour significantly.In the case where loading is along x 2 -direction, the two sides of the cube in the x 1 -x 3 -plane are fully clamped.On one side, the displacement is applied to achieve tension.Since the fibres are not aligned along this loading direction, they have no influence on the stress response.Therefore, the relaxation behaviour is affected, and the stress decrease is significantly slower.This behaviour is caused by the lower stresses inducing viscous deformations, since no fibres are aligned with the x 2 -direction.For the given material parameters, it can be observed that the elastic response differs.The initial slope is significantly higher (material behaves stiffer) due to consideration of fibres.Furthermore, the utilisation of c 3 , c 4 , d 2 , d 3 and d 4 induces an increasing of stiffness, starting at t = 1.5 s.Moreover, the difference can be seen in the stress at infinite time as well.
For a combined loading along x 2 -and x 3 -direction, the boundary conditions are such that the top and bottom side of the cube are fully clamped.The displacement load is applied along x 2 and x 3 simultaneously.The quantity of displacement is kept constant by applying 2.5   √ 2 mm along both directions.As it can be seen, the fibres are loaded for this set of boundary conditions.However, the influence is significantly smaller, since they are not fully aligned along the load.Due to the material parameters c 3 , c 4 , d 2 , d 3 and d 4 being utilised, nonlinear fibre behaviour is obtained.Thus, the stress in this case is significantly lower, compared to that obtained when loading only along x 3 -direction.Nevertheless, higher stresses are reached compared to loading only along x 2 -direction.Furthermore, the relaxation behaviour differs significantly from loading only along x 2 -direction, since fibres are active in the viscous branch.To conclude, the introduced formulation relative to the intermediate configuration is capable of modelling anisotropic viscoelasticity.

Effect of fibre deactivation under compression
To highlight the effects of deactivating the fibres under compression in the proposed formulation, a numerical example is investigated.A cube with edge length 5 mm is considered.The cube is fully fixed at the top and bottom sides.The load is applied as displacements, and both tension and compression of the cube are simulated.The material parameters of the simulations are given in Table 4. Viscous behaviour is neglected for this analysis.Additionally, fibre direction A is activated.The maximum displacement applied is 2.5 mm in both compression and tension.
The results are displayed in Figure 14.The stresses are obtained by dividing the forces measured along the prescribed displacement by the reference area of one side of the cube.It can be observed that the results coincide under tensile loading.This behaviour is expected since the fibres are active for tensile loading.For compression loads, however, the behaviour changes significantly.For Equation (23), the fibres are also loaded while being compressed.Therefore, higher stresses are obtained.This modelling approach does not yield good agreement with reality for the materials investigated in Section 5.4.Additionally, it should be noted that the characteristics utilising Equation ( 23

TA B L E 4
Material parameters for the comparison of Equations ( 23) and (26).

Material parameter
Value Unit from compression to tension.The proposed formulation with fibres deactivated under compression (containing Macauley brackets) is only C 0 -continuous, since the fibres are deactivated for compression loading.Therefore, the material stiffness changes abruptly while transitioning from tension to compression.Nevertheless, the convergence is quadratic as it can be seen in Figure 15.Therefore, the proposed formulation can be considered to be implemented properly.

Application to plants
In this section, the material parameters corresponding to Sorghum bicolor are identified, and the proposed formulation is validated.

Material parameter identification
For the material parameter identification, numerical simulations are carried out.The material parameters are altered until good agreement compared to reality is obtained.The material parameters are tested in a structured way.The parameters associated to invariants with linear dependence on the right Cauchy-Green tensor are tested first.For the sake of brevity, only the first summands for the corresponding invariants of Equation ( 21) are tested.Since a good agreement with reality can already be achieved, no need arises for utilising further material parameters.This highlights the versatility of the model to fit and reproduce even more complicated material behaviours.The geometry of the specimen simulated is a cylinder, which spans 200 mm.The diameter of the idealised geometry is 11 mm.This equals the average value of the samples of the experiments.To simulate the conducted physical experiments, the three nodes at the bottom at x 1 = 10 mm are fully fixed.To obtain a gauge length of 180 mm, the three nodes at the bottom at x 1 = 190 mm are fixed in longitudinal direction of the sample.The displacement is prescribed at the middle of the gauge length at the three nodes at the top of the sample.This set of boundary conditions is chosen, since it represents the three-point flexural test, detailed in Section 2.2, best.A graphical representation of the geometry, the load and the boundary conditions is given in Figure 16.
The discretization consists of 5 elements in radial direction, 20 elements along circumferential direction and 66 elements in longitudinal direction.The top and bottom three nodes resemble the parts of the specimen, which are in contact with the load pin and the supports.The maximum displacement is set to |u 3 | = 5 mm.The load is increased until a simulated time of 30.36 s is reached.Until the end of the simulation, the prescribed load is kept constant.In total, 650 s are simulated.The fibre direction is set to coincide with the primary growth direction of stem.In Sorghum bicolor, the fibres mainly developed around the tissues, which transport water from the roots to the leaves (xylem) and the tissue conducting nutrients (phloem).Therefore, the fibre directions are aligned with the direction of the stem growth.The results of the material parameter identification can be seen in Figure 17.The stresses in this figure are obtained by multiplying the forces measured at the nodes where the displacement is applied by the factor L R 3 .R denotes the radius of the circular cross-section and L is the gauge length.Furthermore, the material parameter set obtained is shown in Table 5.As it can be seen, reasonable agreement to reality is obtained.The characteristics measured in the experiment show a natural difference due to the uncertain nature of the material properties of bio-materials.The simulation is within the range of the experiments.To reduce the natural uncertainty, experimental results of the most similar samples are considered.This is achieved by measuring the distance from the tip of the plant.Only samples in the range of 120-560 mm from the tip of the plant are considered.The experimental results show an immediate relaxation response when no further displacement is applied.The sharp drop also shows that the deformation to this point is mostly elastic.Roughly 60 s after the relaxation start, the stress decrease slows down but continues for the whole experiment time.

Physical validation
Physical validation is performed by using the previously identified material parameters to simulate another experimental test.A bending test at cyclic loading is considered.The geometry is taken from Section 5.4.1.The boundary conditions are TA B L E 5 Material parameters obtained by the material parameter identification process.not changed.The material parameters are taken from Table 5.The displacements, which are prescribed at the aforementioned points, are directly taken from the testing machine.To obtain a continuous function, the values at times between measurements by the testing machine are linearly interpolated.The maximum displacement applied is |u 3 | = 7.999 mm.The results are evaluated until a simulated time of 300 s is reached.The results are shown in Figure 18.The stresses are obtained in the same manner as in Section 5.4.1.It is observed that the agreement to reality is satisfactory.The reduction in stress peaks during every load cycle is captured by the proposed formulation.The larger difference compared to the measurements is due to multiple reasons.Firstly, the biological data is subjected to uncertainty in material parameters.Secondly, the approximation as a circle is not perfect since imperfections occur in naturally grown structures.Additionally, during the tests, plastic deformations occur, clearly visible as a stress decrease during the first cycle.

Material parameter Value Unit
The following cycles show no further pronounced plastic deformation, but a decrease in maximum stress.The upward shift of the numerical simulation compared to the experiment in Figure 18 is caused by plasticity, and needs to be investigated in future research.Furthermore, a plot of the deforming structure over time is presented in Figure 19 for the cyclic test.

CONCLUSION AND OUTLOOK
To summarise, the proposed model is capable of accurately describing the main aspects of the mechanical behaviour of plants.The main findings are that the plants utilised in Section 5.4 show significant viscous deformations, and a fibre stiffness, which is only active under tensile loading.The proposed formulation with respect to the viscous intermediate configuration can represent the significant viscous deformations in a consistent way.Furthermore, the prominent fibre stiffness under tensile loading and resulting anisotropy can be captured by the utilised Yeoh material model.Additionally, the deactivation of the plant fibres under compression can be modelled by the proposed extension to the material model.
Further research can be conducted on the inclusion of multi-physical phenomena.For instance, the viscous behaviour of plants depends on moisture.Additionally, the failure mechanics of the plants can be investigated.Plasticity is also present in this specific material.Therefore, the incorporation of plasticity or viscoplasticity can be subject of further research.The modelling of the failure mechanisms would allow for a more detailed understanding of the mechanical behaviour of plant-like systems. and The elastic tangent, introduced by Equation ( 35), is given by and Here, the additive split C = C vol + C iso is utilised.For the calculation of the elastic tangent, the second-order partial derivatives are stated as and The formulation for the branch subjected to the total deformation is given.For the branch including viscous deformations, the formulas can be utilised as well.The total deformation gradient is replaced by the elastic part F e .The fourth order tensor K is given by For the calculation of matrix K, the deviatoric projection tensor dev( I ) is required.Furthermore, the partial derivative of the Mandel stresses with respect to the viscous deformations is required.The partial derivatives of the elastic deformations with respect to the viscous deformation are stated as and The linearisation term of Equation ( 35) is specified as As a consequence, the specification of the differential of the isochoric part of the Kirchhoff stress tensor with respect to the isochoric elastic part of the deformation gradient remains, since The partial derivatives of the invariants are stated as For the calculation of the derivatives for the Yeoh model with Macauley brackets, the derivatives of the invariants are set to 0, if Īi < 1.Furthermore, the derivative of the elastic isochoric deformations with respect to the total elastic deformation gradient reads (A47)

𝜕F
Finally, the derivative of the elastic deformations with respect to the total deformations is given by

APPENDIX B. INVERSE OF A FOURTH-ORDER TENSOR
In Equation ( 33), the inverse of K, which is a fourth-order tensor with 81 components, is required for computation.Therefore, reordering of the tensor is applied.As a result, a representation of the tensor as a [9×9] matrix is achieved.For the reordering of the tensor, the matrix The resulting [9×9] matrix can be inverted by utilising numerical matrix inversion algorithms.In this case, a non-symmetric residual is differentiated with respect to the non-symmetric part of the deformation tensor.Hence, all 81 entries of Equation (B2) need to be evaluated.

APPENDIX C. ALGORITHMS
In this section, Algorithm 1 provides the update of the viscous deformations.

APPENDIX D. COMPARISON OF THE EIGENSPACES OF B E AND C I
When using the eigenspace to evolve viscous deformations, some modelling assumptions are commonly applied to simplify the update procedure and avoid complicated tensor calculations. 5It is noted that these assumptions render the update procedure of the evolution equations only conditionally stable.The work at hand circumvents these assumptions by avoiding the eigenspace entirely.In this section, the eigenspaces of the elastic left Cauchy Green strain b e and the corresponding associated history variable, that is, the inelastic right Cauchy-Green strain C i , are compared.This is done by considering a simple material model with a Neo-Hookean base, and a similar viscous formulation to the one described in Reference 5.This model used for the comparison of eigenspaces also utilises the same elastic predictor and inelastic corrector type framework to evolve the viscous deformations.A simple simulation is run, with a single linear eight node brick type finite element subjected to combined tension and shear, as shown in Figure D1.The displacement controlled loading is such that tensile strain of 100%, and shear strain of 50% are applied within a duration of 1 s, and held constant thereafter until the total simulation duration of 12 s.Then, at a converged stage in each time step of the simulation, the eigenvectors associated with the tensors b e and C i are stored for the first integration point, and are normalised.Then, the angle between corresponding eigenvectors of the tensors b e and C i is computed.The results are shown in Figure D2.From these results, it is clear that the tensors b e and C i do not share the same eigenspace, because the angle between the vectors is non-zero.

F I G U R E 5 6 7
Comparison of the influence of different boundary conditions on the hysteresis.Convergence of the first time step of the example which is not fully clamped.Influence of alternative material parameters on hysteresis.
The boundary conditions are such that the top and bottom sides are fully clamped in all examples.The applied displacement is 5 mm at the top surface.The boundary conditions of the simulations are depicted in Figure 8.Only tensile loads are considered.The simulated time equals 4 s.The load is applied in the first half of the simulated time.Afterwards, the displacements are kept constant.TheF I G U R E 8 Boundary conditions for isotropic stress relaxation.

11 12
Convergence of the first time step with only holding for a 1 = 10 MPa,  t = 5 s.Boundary conditions for anisotropic stress relaxation (A) Boundary conditions for loading in x 3 and loading in x 2 + x 3 , (B) Boundary conditions for loading x 2 .

) are C 1
U R E 13Comparison of different loading directions.

14 15
Comparison of stress response induced by Equations (23) and (26).Convergence of the first time step of the compression part of the characteristics utilising Equation(26).

F I G U R E 16 17
Boundary conditions of the material parameter identification process and validation.Stress time relation of the material parameter identification process.

F I G U R E 18
Stress time relation of the physical validation.

19
Evolution of the stresses in the cyclic bending test.
The first two indices of K and the last two indices are combined.Each of the pair is input into k and the corresponding index in the resulting [9×9] matrix is obtained.Applying this algorithm, the matrix representation of K is expressed as

Algorithm 1 .
Update of the viscous deformation gradient Require: u, F v n elastic trial

F
I G U R E D1 Single linear brick type finite element subject to tension and shear.
Angle between eigenvectors of b e and C i .
Material parameters for the comparison of different boundary conditions.
TA B L E 1 Material parameters for the comparison of different loading directions.
TA B L E 3 Subsequently, the specification of R F is required for computation.The formulation is stated as