Constitutive modelling of fibre reinforced nonhomogenous hyperelastic materials

In the paper two classes of poly-convex stored energy functions for transversally isotropic hyper-elastic materials have been proposed. Within these classes, a simplified model of hyper-elastic materials with an isotropic matrix reinforced with fibre family was derived and implemented in finite element software ABAQUS/Standard. The two numerical examples proving correctness of the implementation and showing its positive physical features are presented.


Introduction
Large elastic deformations are typical among others in mechanics of composite materials in which polymeric matrix is reinforced with different types of fibres and in biomechanics of soft tissues [1][2][3]. Therefore, in this work, the constitutive relations of these materials within the theory of hyper-elasticity are considered. In case of materials with isotropic matrix reinforced with one fibre family the considered relationships belong to the group of transversally isotropic materials. In the paper, the results of the theory of scalar anisotropic tensor functions representation are used [4]. The material heterogeneity result from the fibre distribution in considered region. Authors aim to indicate such constitutive relationships for which the stored/strain energy (SE) function is a poly-convex function [5,6]. The simplest possible examples of poly-convex SE functions and related constitutive relationships are given altogether with presentation of the basic equations needed for implementation in FEM software ABAQUS/Standard [7][8][9]. The issues presented here may be treated as a continuation and generalization of concepts presented by authors in papers [10][11][12][13] related to modelling of fibre reinforced materials.

Constitutive relationships for transversal isotropy
The SE function W is considered as a function of deformation tensor is a deformation tensor gradient and det 0 J ! F , while "det" stands for the tensor determinant) with additional tensor parameter M m m. Unit vector m X is tangent to specific chosen fibre direction, and X is a vector determining the location of material point in reference configuration. Analysed material in reference configuration is a transversally isotropic one [4]. When the SE function will be postulated in the following form ˆ, , tr , tr Cof ,det , tr , tr Cof then from basic relationships of continuum mechanics and hyper-elasticity the following constitutive relationship in Lagrange description for the second Piola-Kirchhoff stress tensor may be determined: In (1)-(3) symbols "tr" and "Cof" stands respectively for tensor trace and tensor cofactor (transposition of adjugate tensor) [6]. From natural state condition results that: ˆ3 ,3,1,1,1 0 W and 3,3,1,1,1 0, 1, ,5 After substituting (2) hold.

Poly-convex stored energy functions
In the literature on transversally isotropic hyper-elastic materials the constitutive relationships in the form resulting from application of the results presented in Spencer's monograph [1] are used. In that case, a slightly different than in (1) invariant basis of tensor C are applied. It is worth underliynig that invariants used in SE function (1) have clear geometrical interpretation (connected with deformation of material line, surfece and volume elements) and are very convinient for polyconvexity testing of SE function.
The SE function , W X F is a poly-convex one if for every X the extended function are poly-convex when 0 i a ! and 1 j t D , etc. The function assumed as convex function tending to infinity for J tending to zero and infinity as well. A poly-convex SE functions other than (10) and (11) may be found in [14][15][16][17][18][19].
where 1 2 , , M I I J < is SE function of matrix (isotropic function), Zn < are SE functions of fiber families, while n p stands for the volume fraction of the fiber family in unit of material volume. Functions of SE accumulated in n-th fiber family that is working axially in direction determined by n m X is approximated as: and parametric tensors n n n M m m are defined in reference configuration. Parameter zn E may be interpreted as the initial Young's modulus of n-th fiber family.
Among others in works [10,12,13] the special case of SE function, which corresponds to the Ciarlet model for incompressible materials discussed in [6]: was considered. For

UMAT procedure
Here we are proceeding analogically as shown in [10][11][12][13]. The model presented in section 4 was programmed as the UMAT procedure, which requires the formulation of constitutive relations (5) in incremental form. After counting Kirchhoff's tensor increment (as defined by the objective Jaumann derivative), In case of function (14) the forth order tangent tensor operator M C have the following form: The UMAT material procedure was tested among others on boundary value problems with homogenous stress and strain fields and known analytical solutions. The main purpose of this type of tests was to demonstrate the correctness of the implementation and the equivalence of the relationship (5) (with potential (11)) and the incremental relationship (15) with the operator (16).

Numerical testtension test of a plane strain quarter of circular disk
In this section the problem of hydrostatic tension of the plane strain circular disc (with radius a) is studied. The disc is reinforced with fibres symmetrically placed in relation to two planes. Due to the symmetry of the problem, the modelled region is quarter-circle shape. Heterogeneity of the disc is due to different orientation of the preferred specific fibre direction in material. There is seven variants of placement of reinforcement fibre families for which the problem is solved: i) fibers constantly inclined to 1 x axis with an angle 4 S , The region of the quarter of circular disk was modelled with 9516 CPE4 elements (plane strain, four nodes with linear shape functions). The constitutive model for material is described with (15) and (16) where the following parameters are substituted: The hydrostatic loading was assumed as multiple value of the initial shearing modulus equal to q =100.0 0 P . On the basis of the analysis of the examples presented in Figure 1 it can be stated that only three qualitatively different solutions were obtained. The solutions for variants iii) and iv) are analogous: the circle has deformed into an ellipse (there is a homogeneous orthotropic material in the plane), except that the extreme deformation values in the first case are observed in the direction of the 2 x axis while in the second case in 1 x direction. With respect to variants v) and vi) the solutions are qualitatively identical and very similar to the solution obtained for non-reinforced material. In these three cases the values of stresses, strains and displacements in some particular locations are not the same. In the case of solution for a radially-reinforced disk (i.e. variant vi) the deformation shown in Figure 1 corresponds to a load value of 0.748 q . This is the last value at which the integration algorithm with the same steering paremeters for all variants was able to find a solution with a certain assumed accuracy. In tasks i) and ii) there is solution for the orthotropic circular non-homogeneous disk. In the paper [21], the boundary problem of hydrostatic tension test for sphere made of transversally isotropic hyper-elastic material was considered. For this type of material available in ABAQUS [9], non-physical results were obtained. The sphere deformed into a sphere. The reason for this is the application of too simple constitutive relationships for material components (matrix and fibres) in case of hyper-elasticity of fiber reinforced materials, in which additive decomposition of SE function is assumed.

Numerical testtension test of a flat bar with heterogeneity
Next numerical example refers to the tension of the thin flat bar with reinforcement placed in plane 1 2 OX X according to sinus function, as a one fibre family. The flat bar is modelled as 3D problem, with thickness equal to g, length equal to l=50g and height h=12.5g, see applying displacement 1 u =25g on that boundary. On all other surfaces, zero stress boundary conditions are assumed. The flat bar region is modelled with 5000 C3D8 elements (eight-node elements with linear shape functions) 100x50x1 (length x height x thickness).The constitutive model of the material is the same like in previous example and the following material parameters are taken: In (18) symbol 2 p X stands for the offset of the function in the direction of the 2 X variable.
It is worth noting that parameter 2 p M m m. In initial configuration 1 1 = X x and 2 2 = X x . The examplary results obtained for the analysed problem in the form of a contour graphs are presented in Figures 3-4. h=12.5g

Final remarks
Two classes of poly-convex SE functions for transversally isotropic materials have been proposed. Within these classes, a simplified model of hyper-elastic materials with isotropic matrix reinforced with fibre families was proposed. In particular, the proposed SE functions are simplifications of the SE functions applied in our earlier works [10,12,13]. These models do not lead to non-realistic results, as shown in [21], where the transversally isotropic material models available in the ABAQUS program [9] were used. The heterogeneity of the material result from the different arrangement of the fibres in the matrix material. The fibres orientation were defined by the parametric tensors, which are used as arguments in the constitutive relation. This facilitates the numerical implementation in the finite element method programs. For example, in the ABAQUS/Standard software, for programming the material model the UMAT user procedure was used.