NONLINEAR COMPOSITES STRUCTURAL ANALYSIS WITH VISCOELASTIC FIBERS

This paper presents a formulation for material and geometrical nonlinear analysis of composite materials by immersion of truss finite elements into triangular 2D solid ones using a novel formulation of the finite element method based on positions. This positional formulation uses the shape functions to approximate some quantities defined in the Nonlinear Theory of Elasticity and proposes to describe the specific strain energy and the potential of the external loads as function of nodal positions which are set from a deformation function. Because the nodal positions have current values in each node, this method naturally considers the geometric nonlinearities while the nonlinear relationships between stress and strain may be considered by a pure nonlinear elastic theory called hyperelasticity which allows to obtain linearised constitutive laws in its variational form. This formulation should be able to include both viscoelastic and active behavior, as well as to allow the consideration of nonlinear relations between stresses and deformations. It is common to adopt hyperelastic constitutive laws. Few are the works that use the strategy of approaching the problem such as fibers immersed in a matrix. The immersion of fibers in the matrix makes it possible to include both a viscoelastic behavior in a simple and direct way. The examples are simple cases, some of them even with analytical solutions, mainly for validation purposes of the presented formulations. By modeling a structure, the examples show the potentialities of the concepts and proposed formulations.


INTRODUCTION
According to Tang et al. [1], the modeling of the mechanical behavior of a muscle tissue may follow a microscopic approach [2] or macroscopic [3].In general, the macroscopic model is used to describe the passive behavior that corresponds to soft tissue and muscle fibers, whereas the microscopic model is used to model the active behavior that corresponds exclusively to muscle fibers [4].The mechanical properties of the active and passive parts are the basis for the development of muscular numerical models.
According to Chagnon et al. [5], elastic behavior is generally what dominates the response of biological tissues.Thus, nonlinear passive behavior of fibers and soft tissue is usually modeled through Continuous Mechanics using an energy approach known as hyperelasticity.The hyperelasticity consists of writing the specific energy of deformation as a function (mainly) of the deformations tensor invariants, and it is a very general and flexible approach with respect to the behavior of the material, since it allows to consider, besides great deformations, also phenomena of damage, Viscoelastic effects, incompressibility and even anisotropy [6].Deformation energies expressed as a function of the invariants I 1 , I 2 and I 3 result in isotropic behaviors, which are only justified for the matrix of muscle tissue.
The physical behavior resulting from the presence of the fibers, which is a kind of anisotropy, can be taken into account considering the preferential orientation of the fibers through a vector oriented in the direction of the fibers [7].In this context, the simplest type of anisotropy, transverse isotropy, is present in many biological tissues [8] and can be considered simply by introducing directly into the deformation energy equation a portion of additional deformation energy that depends on the preferential direction of the material in the reference configuration, as implemented by some authors [1,9,10].The definition of new invariants as a function of the directions in which the immersed fibers are oriented allows to define additional plots of deformation energy in a compact form, even allowing to consider materials with fibers orientated in different directions [11].Thus, for a transversely isotropic muscle model, the deformation energy also depends on other two additional invariants (I 4 and I 5 ), expressed as a function of a vector oriented in the preferred direction of the fibers.
The active behavior, in turn, is the very tension or deformation of the fibers written explicitly through an activation function.Since muscle fibers have both active and passive behavior, their behavior can be given by the parallel superposition of their active behavior with their own passive behavior, according to the model proposed by Hill.Finally, according to Calvo et al. [12], one has to superimpose the effects of the Active muscle (muscle fiber) behavior with passive muscle behavior (matrix and muscle fiber) to describe the behavior of muscle tissue as a whole.
In finite-element analyses, passive muscles are usually modeled by solid three-dimensional elements, whereas active muscles are represented by single-dimensional single-bar models [13].Therefore, it is reasonable to propose the modeling of muscle tissue by immersing fibers represented by one-dimensional elements of simple bar, which follow the Hill muscular model and transmit forces when activated during the simulation, in a discretised matrix through solid elements with nonlinear relation between tensions and deformations, also including possible viscous effects.
Although muscle tissues also exhibit viscoelastic properties, few authors consider viscous effects in their models [13,14,15], and even when they do, they generally adopt the Maxwell Model instead of the Kelvin-Voigt model.Although the Maxwell model allows to express the phenomenon of relaxation simply through a written rigidity as a function of time, the creep phenomenon described by the Kelvin-Voigt model is much more representative of the viscous muscular effect, because through this it is possible to simulate an effect of muscle "fatigue", in the sense that a muscle contracted for a certain time interval no longer has the same muscular capacity as that of the beginning of the contraction.
The present work is based on three pillars of biomechanics, which according to [8], are: Continuous Mechanics (non-linear elasticity), tooling Mathematical methods (finite element method, for example).

POSITIONAL FINITE ELEMENT METHOD
The main references consulted about the finite element method were the following: [20÷26].In relation to the positional formulation of the method, the reference [16] were the main sources of information.In this context it is also important to highlight the work of [27], since it is one of the precursors in the positional approach of the finite element method.

Finite element of plate (2D)
In the positional finite element method, three configurations are used to formulate a finite element: the final deformed configuration, the initial undisturbed configuration, and the dimensionless reference element.Figure 1 illustrates the map- ping between the dimensionless reference domain (Ω 0 ), the undeformed initial configuration (Ω 1 ) and the deformed final configuration (Ω) using the configuration change functions (f, f0 and f1) and their respective gradients (F, F0 and F1).
By adopting a cubic approximation for the configuration change function f1 (ξ 1 , ξ 2 ) that maps all the particles of the continuum in the domain of a triangular element with respect to the dimensionless coordinates ξ 1 , ξ 2 and at a given instant, we have For directions 1 and 2 in the plane the following equation presented in matrix form: (1) Where [ϕ 1 ] and [ϕ 2 ] are complete polynomials of degree 3, whereas the terms a 1 n and a 2 n are vectors that contain the constants of the function used to approximate the positions.Viscoelastic behavior is a phenomenon of non-reversible internal energy dissipation that results in the existence of unbalanced states that evolve over time.The relaxation and creep phenomena are two typical viscoelastic effects, where the return to a (new) equilibrium position after a disturbance depends on the time.Relaxation describes the decrease in tension under a constant deformation regime, while creep describes the increase of deformation under the application of a constant tension.The models of Maxwell and Kelvin-Voigt are two typical mechanical models used in the description of the phenomena of relaxation and fluency, respectively [16].
In order to take into account the viscous behavior (relaxation or fatigue phenomenon), this work proposes a viscoelastic constitutive model for the simple bar element based on the Kelvin-Voigt model.This model, shown in Figure 2, consists of an elastic member in parallel with a damper element.Due to the effect of the damper on the model, the deformations do not arise immediately after the application of a load.This is because at the moment of application of the load, the entire force is initially supported by the damper, which, unlike the elastic element, does not respond with immediate deformations upon application of a stress.However, over time, this damper tends to "relax", gradually transferring the load to the elastic element.
According to Figure 2, it is additionally considered an active element of contraction series with the viscoelastic model of Kevin -Voigt.The positioning of the element was chosen strategically so as to allow immediate contraction of the fiber element.Thus, the expected behavior is that a prescribed deformation in the contraction element immediately generates a contraction in the single rod element (instant deformation) without the contraction element interfering with the viscoelastic model or vice versa.
In the proposed model, the traction P in the fiber element is equal to the traction P c in the contraction element, which is equal to the viscoelastic traction P ve , which in turn is the sum of the elastic traction P e with the viscous traction P v , that is, P = P c = P ve = P v + P e .Then, adopting a linear elastic model for the behavior of the elastic element such that P ve = Eɛ̌v e , where E is the modulus of elasticity (Young), and still c is the viscous constant of the model, ɛ̌t +1 ve and ɛ̌t ve , respectively, the deformations at the current instant and at an earlier instant, have by the finite difference method for a time interval ∆t: (2) Since the total deformation ɛ̌t +1 of the model, at time t + 1 and at the previous instant t, is equal to the deformation ɛ̌c prescribed in the contraction element plus the viscoelastic deformation ɛ̌v e , ɛ̌ = ɛ̌c + ɛ̌v e : (3) In Equations 2 and 3, the deformation ɛ̌c prescribed in the contraction element, both at the current instant t + 1 and at the previous instant t, is known.The total deformation, ɛ̌ of the fiber element, in turn, is iteratively approximated by the method.

NUMERIC PROCESSING
In this section examples are presented, two with the purpose of verifying the proposed formulations and one to demonstrate the applicability of these formulations in the modeling of muscular tissues.
The first example is a validation of the single bar finite element formulation that uses the socalled nonlinear engineering deformation measure and considers viscoelastic behavior and activation capacity (contractions) of these single bar elements.The second example is a verification of the simple bar finite elements immersion strategy in the plate elements as proposed by Vanalli (2004) and presented by Sampaio, Paccola and Coda (2013) and Sampaio (2014).
The third example, besides a comparison with the example proposed by Baiocco, Coda and Paccola (2013), illustrates some of the advantages of the formulation adopted in the present research for representation of muscular tissues when compared to the results obtained in the work of the mentioned authors.

Bar with active viscoelastic constitutive model
The validation of the implementation of the active viscoelastic constitutive model in simple bar elements using the nonlinear engineering deformation and the Piola-Kirchhoff stress, respectively as measures of deformation and stress, was performed from a structure composed of three elements of single bar: The supported properties for all bar elements are as follows: Modulus of elasticity E = 2.1 × 10 6 , viscous constant c = 2.1 × 10 3 and cross-sectional area S = 1.0 × 10 -3 .In the proposed model, the cross-sectional constraint of the bar is neglected (a reduction of the cross-sectional area could be estimated by an incompressibility hypothesis).Thus, in the suggested formulation, the Cauchy stress and the Piola-Kirchhoff stress coincide.The first Piola-Kirchhoff stress is a very interesting stress measure, because it is based on the initial area, which is known.This example was solved in 23,000 steps and adopting a constant time interval T = 1.0 × 10 -6 .The viscous updating is done through equation 4 and equation 5 using the finite difference method.
In Figure 3, node 1 is prevented from moving both in the x-direction and the y-direction, while node 2 is restricted only in the y-direction.In node 2 a horizontal force is still to be determined so that the vertical displacement of node 3 is equal to 0.25 (down), according to Figure 4: The purpose of this example is to present the graph of elastic and viscous stresses for the bar element 1, as well as the horizontal and vertical displacements of the node 3 and also the horizontal displacement of the node 2. This example is particularly interesting to verify that, in the proposed formulation, the bars 2 and 3 are free of stresses throughout The loading process, since these only suffer pure rigid body movement.The expected final length for bar 1 is: As viscosity is a transient phenomenon, which tends to disappear over time, the force required to cause the proposed elongation of 1.1180 in bar 1 (after the disappearance of the viscoelastic phenomenon) is: (5) The Piola-Kirchhoff stress (elastic and viscous plots) obtained in bar 1 with the developed program is shown in Figure 3, where the elastic The graph of Figure 5 matches the expected result, i.e., the viscous stress part of the stress initially applied P = 247,80 / 1,0 × 10 -3 = 2,478 × 10 5 and tends to zero.On the other hand, the elastic stress starts from zero and tends to 2,478 × 10 5 , which is the applied stress.As the area remains unchanged, the initial viscous Cauchy stress and the final elastic Cauchy stress coincide.The 1st Piola-Kirchhoff stress (elastic) obtained in the bars over the 23,000 time steps is as follows:Thus, it can be seen in Figure 6 that the bars 2 and 3 are free of stresses throughout the process, since they have pure rigid body movement.The displacements of the nodes in the horizontal (x) and vertical (y) directions are shown in Figure 7: Thus, as expected, the horizontal displacement of node 3 is the half of the horizontal displacement of node 2. The vertical displacement observed at node 3 was 0.245 (down), which cor-responds to a value of 2 % below the expected value of 0.25.Note that part of this error is due to rounding of the elongation and force calculated in Equations 2 and 3, respectively.It is also noted that, just like the elastic tension, and due to the viscoelastic behavior, the displacements at the nodes happen gradually.Applying the force of 247.80 calculated in Equation 2 in the structure of Figure 3, the final elastic (non-linear engineering) strain of bar 1 obtained by the program was of 0.118, that is, exactly the expected value of A contraction is now introduced into the bar 1 of 0.118 at the 8,001 ° time step such that the structure returns to the initial position given in Figure 3.This shrinkage will then be removed at the 15,001 ° time step, so that the structure returned to the configuration of Figure 4. Most likely, the displacement of 0.245, slightly below the exact 0.25, ensured that the structure returned to Even activating and deactivating bar element 1, the tensions obtained in the bars were exactly those of Figure 6.The contraction of the bar element, therefore, does not alter neither the elastic tensions acting nor interfere in the exponential behavior characteristic of viscoelastic materials or vice versa.The displacements obtained at nodes 2 and 3 were as follows: Figure 4 again reinforces that the contraction applied on the member is immediate.It can be said therefore that the displacements induced contraction element are "hard" and immediate.
For a single bar element, equation 2 proposed for the model is immediately programmable in spreadsheets, allowing even to consider the active behavior.Using the same parameters of the present example both considering and disregarding the bar element activation, the displacements and stresses observed with Eq. 3 programmed in a spreadsheet were the same as those observed for the bar element 1 of this example, which validates the implemented implementations for the model.

Reinforced matrix with fibers
This example was proposed for validation of the coupling model implementation between fi-ber and matrix.Figure 9 consists of a unit sides square where 150 bar elements are distributed (50 elements distributed in the vertical direction and each of these 50 elements divided into 3 other elements of equal size) with 2 nodes each, totaling 200 nodes of bar elements.These bar elements are distributed evenly spaced in a discretised array in 200 plate elements of 10 nodes each, totaling 961 nodes of sheet metal elements.
In Figure 9, the two horizontal faces are prevented from moving in the y-direction, while the left vertical face is prevented from moving in the x-direction.In the right vertical face, a displacement of 1.0 × 10 -3 is applied.The purpose of this example is to calculate the horizontal force F on the right side necessary to cause the prescribed displacement of 1.0 × 10 -3 and to compare with the analytical result.
The constitutive model chosen was that of Saint Venant-Kirchhoff in the form (6) for both matrix and fibers.Therefore, the relationship between stress and strain is linear.The supported properties for all bar elements are as follows: modulus of elasticity E = 1,05 × 10 11 and cross-sectional area A f = 0.4 × 10 -3 .The properties used for the matrix material (sheet metal elements) were chosen arbitrarily.Are they: The parameters c 10 , c 20 of the Saint Venant-Kirchhoff model are the constant of Lamé, which can be written as a modulus function of elasticity to shear G: (8) The numerical result obtained was 4.13964 × 10 6 , which means a difference of less than 2% over the calculated analytical value.Figure 10 shows the displacements and tensions obtained for the present example.
In Figure 10 we observe increasing displacements along the x-axis, in addition to almost zero displacements in the y-direction (on the order of 10 -10 units in length), since the Poisson constant ν was adopted with a value equal to zero, so there is no "coupling" between the x and y directions.Furthermore, the normal stress along the x-axis varies between 1.8657x10 8 and 2.1323x10 8 , while the stress obtained in the ydirection were of the order of 10 4 .Thus, the stresses in the y-direction can be considered zero if compared to magnitude of the stresses in the x-axis direction.
Figure 11 illustrates the Cauchy stress obtained in simple bar elements.In the case of a single rod element with linear approximation for the positions, the tension is constant along the same.
With the results presented in agreement with the expected behavior, everything indicates that the procedure of immersion of the fibers has been correctly implemented.

Biomechanical Arm
This example demonstrates the applicability of the positional formulation of the finite element method according to subsection 2 and the fiber immersion technique presented in subsection 2 in the modeling of biological structures.The present example demonstrates the potential of the active viscoelastic model (contraction / elongation) for single bar elements using nonlinear engineering deformation measurement as developed in subsection 2.
In order to prove the superiority of other hyperelastic models with respect to the hyperelastic Saint Venant-Kirchhoff model, a limitation of the Saint Venant-Kirchhoff model is presented here.In this way, we intend to illustrate very clearly the improvement in the quality of the analyses that can be obtained with the use of hyperelastic constitutive laws in a positional finite element code.
Figure 12 shows a model of an arm where are embedded 7350 bar elements of 2 nodes each, totaling 7420 nodes of bar elements.These bar elements are distributed in a discrete array in 640 plate elements of 10 nodes each, totaling 3031 plate element nodes.The blue sheet members (elements 1-200) represent hard tissue, while the red sheet members (elements 201 -640) represent soft tissue.The 7350 single bar elements are equally distributed in 3675 upper fibers (elements number 1-3675) and 3675 lower fibers (elements number 3676-7350).In Figure 12, the left vertical face is prevented from moving in both the x and y direction.The model of Figure 12 is the same as that used in Baiocco, Coda and Paccola (2013).The purpose of this example is to encompass all concepts studied throughout this work and to demonstrate the applicability of formulations developed.
The constitutive model chosen was that of Saint Venant-Kirchhoff in the form for both soft tissue (red elements) and hard tissue (blue elements).Therefore, the relationship between stress and strain is linear.The modulus of elasticity adopted for the upper fibers is E = 0.3 × 10 6 and for the lower fibers is E = 0.3 × 10 3 .The cross-sectional area of all fibers is S = 1.0 × 10 -1 .The properties used for hard tissue (blue elements) are: (9) The properties used for soft tissue (red elements) are: (10) The parameters c 10 , c 20 of the Saint Venant-Kirchhoff model are the Lamé constant, which can be written as a function of the modulus of elasticity to shear The horizontal force F of Figure 12, which is 320, is uniformly distributed on the face of the two hard tissue elements at the end of the forearm.In Fig. 13 the results obtained with fibers formulated from single bar elements using the Green-Lagrange strain and the 2nd Piola-Kirchhoff stress tensor and single bar elements using non-linear engineering and the 1st Piola-Kirchhoff stress tensor.
The result obtained using single-bar elements with Green-Lagrange deformation and the 2 nd Piola-Kirchhoff tensor stress for the fibers is practically identical to that obtained by Baiocco, Coda and Paccola (2013).In the case of small displacements, this result is close to the solution obtained with single bar elements using nonlinear engineering deformation and the 1st Piola-Kirchhoff tensor stress (about 0.2% smaller), as shown in Figure 13, because, in the regime of small displacements (and deformations), such measures of deformation and tension are confused.

Neo-hookeno model
In order to observe large displacements and moderate deformations, the modulus of elasticity E of the upper fibers is reduced from E = 0.3 × 10 6 to E = 0.3 × 10 4 and the properties of the plate elements in blue (hard ) and red (soft tissue) maintained the same as 17 and 18, respectively.
Based on the plot of energy of deformation ½E∈ 2 , a neo-Hookean hyperelastic model of the form is proposed: The results obtained with the Saint Venant-Kirchhoff hyperelastic constitutive law and with the neo-Hookeian hyperelastic model proposed in Equation 11 are compared in Figure 14.In both cases the fibers are considered simple bar elements formulated from non-linear deformation engineering and the first stress tensor of Piola-Kirchhoff.
In Figure 14 it can be seen that the Saint Venant-Kirchhoff constitutive law permits the self-intersection of material, whereas the proposed neo-Hookeian hyperelastic law, although simple, is physically more consistent (even if the displacement is about 9% higher).The autointersection in the Saint Venant-Kirchhoff model disrupts the convergence of the program, so the greater ease of convergence of the neo-Hookeian model is remarkable.For Figure 14, the displacement tolerance of 1.0 × 10 -8 with the Saint Venant-Kirchhoff model was reached in 13 iterations of the Newton-Raphson method versus 8 iterations of the proposed neo-Hookeian model.
Figure 15 shows the Cauchy stress in the x and y directions for the neo-hookeano hyperelastic model.
Figure 15 shows compression and tensile stresses in the bone, respectively, in the lower and upper portions of the biceps, as well as the left and right side of the forearm.
These results are in agreement with the calculated displacements and the direction of the applied load.Note also that, since the soft tissue stiffness is less than the hard tissue stiffness, the tensions in the muscular matrix are much smaller than the tensions in the bone.

Viscoelasticity
Maintaining the modulus of elasticity of the upper fibers at E = 0.3 × 10 4 and the lower fibers at E = 0.3 × 10 3 and adopting the properties for sheet metal elements in blue (hard tissue) and red (soft tissue) , respectively, according to 9 and 10, the displacements of the structure of Figure 12 are determined under the action of the horizontal force F of 320 over 500 time steps using a time step ∆t = 0.5 × 10 -4 and a coefficient of viscosity c of 15 for the fiber elements.The constitutive law adopted for all plate elements is the neohookena proposed in Equation 11.The fibers are considered single bar elements formulated from the nonlinear engineering deformation and the 1st Piola-Kirchhoff stress tensor.In Figure 16 it is noted that it was possible to introduce a viscous effect into the structure from the consideration of fibers formulated from viscoelastic single bar members as per subsection 2. At the end of the 500 steps of time, the displacements in the x direction are those shown in the legend of Figure 16, these being in accordance with those obtained in Figure 14.The displacements shown in Figure 16 refer to the time steps 1,10,20,40,70,110,160,220,290,370,460,500. Figure 17 shows the displacements (in modulus) in the x and y direction of node 3031 at the end of the arm.
Note in Figure 16 and Figure 17 that in time step 160 (time 290 × 0.5 × 10 -4 = 0.0145 in Figure 17), the structure has practically reached the final equilibrium configuration, more change of position from this point.

CONCLUSION
The numerical examples presented clearly demonstrate the potentialities of the proposed single bar element with viscoelastic properties and active behavior (with contraction / elongation capacity).The use of nonlinear engineering deformation and the 1st Piola-Kirchhoff tensor in the formulation of this element, besides facilitating the physical interpretation of the results obtained from numerical simulations, also simplifies the introduction of the viscoelastic properties and the consideration of the active behavior of the same.
It has also been shown that this proposed element can be immersed in elements for the simulation of fiber-reinforced composite materials.This immersion strategy allowed the viscoelastic properties and the active behavior present in the simple bar elements immersed would manifest themselves in the flat structure as a whole and not only in the bar members themselves.
The modeling of biological structures or the simulation of composite materials reinforced with are some of the possible applications of this proposed element.
As the whole development of this work was carried out using the proposed formulations are already automatically suitable for simulations involving large displacements, and since the positional formulation of the finite element method allows the adoption of hyperelastic constitutive laws in an immediate way, it is also possible to consider large deformations, nonlinear relations between stresses and deformations, to include incompressibility behavior, besides proposing physically consistent constitutive laws, for example that do not allow degeneration or self-intersection of material.

Fig. 2 .
Fig. 2. Active viscoelastic model of the single bar finite element

3 Fig. 8 .
Fig. 8. Displacements of nodes 2 and 3: a) Displacement in the x-direction, b) Displacement in the y-direction

Fig. 13 .
Fig. 13.Displacements in the x direction: a) Fibers with deformation of Green-Lagrange, b) Fibers with nonlinear engineering deformation

Fig. 15 .Fig. 16 .Fig. 17 .
Fig. 15.Cauchy stresses in the x and y directions: a) Cauchy stress in x-direction, b) Cauchy stress in the y direction