INTRODUCTION
Because thermodynamic variables such as entropy and temperature are assumed
constant in hyperelastic materials, which are defined as a sub-group of elastic
materials, theory of such materials is viewed as mathematical abstraction. However,
if application range is correctly determined, hyperelasticity theory can be
successfully applied in various branches of engineering to solve problems in
a wide range. Composites such as fiber reinforced rubber like materials and
biologic tissues demonstrate a rather complex non-linear behaviour. Most of
the composites in this category can be regarded as anisotropic non-linear hyperelastic
materials. Car tires, lots of biologic tissue types and fiber reinforced rubber
panels belong to this class of materials. Generally, due to its simplicity and
usefulness, hyperelasticity is used to model constitutive behaviour of rubber
like materials. Previous studies used in the incompressible hyperelasticity
model have been mostly restricted to isotropic materials[1,2].
For hyperelastic materials undergoing finite deformation, constitutive behaviour has been characterized by W, strain energy density function dependent on the local deformation gradient tensor. Because scalar values dependent on strain are invariant under coordinate transformation, they would provide an advantage if used as arguments of W. According to the reference configuration, three scalar parameters are required for materials demonstrating isotropic behaviour. Typically, W is either defined through principal invariants of deformation tensor, like Rivlin[3] did or in terms of stretches, like Ogden[4-6] did. A physically meaningful constitutive formulation has been developed dependent on three invariants of natural strain to correct the constitutive definition of an isotropic hypereslastic material[7].
Browsing though sources on the subject, one can notice that Hoger et al.[8] in their studies on anisotropic non-linear elastic materials, have developed a constitutive equation for finite deformation of a surplus stress transverse isotropic hyperelastic material. Meschke and Helnwein[9] have used rebar elements to represent behaviour of fiber reinforced composites. Rebar element has been used in the finite element analysis of finite deformation of a fiber reinforced rubber panel. By adding extra terms to the function of isotropic strain energy, Srinivasan and perucchio[10] have presented an anisotropic constitutive equation. However, their studies are restricted to a special class of materials and results in some problems related with reinforcements in general curvilinear coordinates. Due to the mathematic difficulties they involve, curvilinear coordinates are not widely used in constitutive analysis. However, curvilinear material coordinates can be effectively used to model the existence of reinforcements at a certain material point made into random curvilinear directions. Constitutive equation for anisotropic non-linear hyperelastic materials in general curvilinear material coordinates has been developed by Park and Youn[11]. Mixed finite element method has been used in this study, selecting a suitable system of coordinates. The model developed, based on the result expressions of these research studies, can be effectively used to analyze biologic tissues and composite materials similar to fiber reinforced rubber[11].
In this study, in the framework of continuum mechanics axioms, a constitutive equation has been obtained characterizing non-linear mechanical behaviour of an arbitrary fiber reinforced hyperelastic material. Medium has been assumed homogenous and incompressible while fiber family has been assumed in extensible. After presenting information on the fiber geometry and kinematics equilibrium equations have been summarized, proceeding with the constitutive formulation[12,13]. While constitutive functional was found based on the axioms of the constitutive theory, especially the axioms of objectivity and material symmetry and on concepts concerning the symmetry group of a material, findings of the theory of invariants have been used as a routine method to materially determine arguments of the stress potential[14]. As a result, constitutive equation of stress has been obtained in both material and spatial coordinates. The constitutive equation obtained can help determine mechanical behaviour of fiber reinforced rubber like materials.
As a scientific application, considering three different case, ANSYS, a powerful finite element software has been used. In the first case, two hyperelastic materials of Money-Rivlin type have been examined under a suitable geometry and load conditions. In the second case, a hyperelastic material was reinforced using steel fibers in axial direction and resolved. In the third case, solution has been made for an A 36 steel material under same geometry and load conditions and results have been compared. In this last case the geometry was constituted completely of steel material.
FIBER DEFORMATION GEOMETRY AND KINEMATICS
A material curve representing the fiber family passes through every point of the composite material considered in this study. This fiber family is represented by a continuous vector field A(X) before deformation and by a continuous vector field a(x) after it (Fig. 1). Transition of the fiber field before deformation to the positions assumed after deformation occurs in the way of dragging along with the medium.
For a differential fiber length dL after deformation, the following relationship can be written in terms of dL before deformation.
|
| Fig. 1: |
Fiber curve before and after deformation[15] |
Expanding the right side of the relationship around the point X using Taylor series, neglecting terms with a high grade, will lead to the relationship between fiber lengths before and after deformation:
Substituting the following expression (3),
in the Eq. 2 will yield the following expressions:
Deformation geometry pertaining to fiber family is expressed by relationship (4). Since A and a vectors here represent unitized fiber vectors before and after deformation, operations,
are valid. (Expression (4) have been considered here.) Accordingly, expression,
can be written.
EQUILIBRIUM EQUATIONS
During interaction with the surroundings of the medium we have taken into consideration
that occupies volume V(t) at the moment t, if integral equations governing the
global balance equations are localized using “updated” Gauss and Stokes
theorems, the following differential expressions are obtained[2,16]:
Here, ρ0, ρ, Vk, t, J, Vp, fp,
tkl, εkrp, qk, h, ε, η, θ
and ργ, respectively represent, density before deformation, density
after deformation, velocity vector, time, jacobian, acceleration vector, body
force density per unit volume, stress tensor, permutation tensor, heat flux
vector, heat source per unit mass, internal energy density per unit mass, entropy
density per unit mass, absolute temperature and entropy generation per unit
mass any time t.
Local energy Eq. (11) and entropy inequality (12)
are then suitably combined and, using a Legendre transformation such as ψ
≡ ε-θη for the free energy, entropy inequality is obtained
as follows in the material form:
Relationships between material and spatial forms of the values appearing in this inequality are given as follows:
Here Σ stands for thermodynamic stress potential, Ψ: for generalized free energy density, CKL: Green deformation tensor, dkl: deformation rate tensor, xk,K: deformation gradient, XK,k: deformation gradient of inverse motion, TKL: symmetric stress tensor in material coordinates and QK: heat vector in the material coordinate system.
CONSTITUTIVE MODEL
To enable the use of the inequality (13), the general expression of entropy production, independent variables should be known, on which thermodynamic potential Σ depends as well as the way of such dependence. Accordingly, selection of arguments of Σ would formally mean selection of a certain material. In this study, arguments of Σ and variables it depends on have been found using constitutive axioms, according to the material selected. Assuming that, in a considered material, stress potential Σ of a material point X at a point in time t depends on the motion and temperature background of all material points constituting the object, arguments determining the stress potential of a fiber reinforced material under consideration and the constitutive equation of the stress have been obtained as follows:
Placed among the arguments of Σ, A here represents the fiber vector field before deformation. Tkl is a component of the stress tensor in material coordinates. Equations (20) clearly shows us that the stress tensor is derived from the stress potential Σ. In this situation, appearing as a constitutive function with clear arguments, open form of Σ should be obtained.
However, firstly we should revise the restraints brought upon the constitutive
function of the material in consideration by the material symmetry axiom. Because
any isotropic material belongs to the fully orthogonal group, the symmetry operation
would yield the following property:
Accordingly, each material coordinate conversion matches an orientation. The said conversion would be as follows for each
leaving the constitutive functional form-invariant. Mathematically this means that the following conversion would be valid:
On the other hand, inextensibility of the fibers and incompressibility of the
medium enjoys broad recognition in practice in terms of formulation. Provided
that the medium is assumed to be incompressible and fiber family-inextensible,
the following conditions would be respectively satisfied:
and
Under these conditions the constitutive equation of stress would be obtained as follows in material and spatial coordinates:
p and Ta in these expressions are Lagrange coefficients defined by field equations and boundary conditions. Deriving J =1 from the incompressibility condition, spatial form of the stress tensor can be expressed as follows, using expression (15):
STRESS DETERMINATION
Because matrix medium is assumed to be isotropic, as understood from the theory
of invariants, the stress potential Σ should depend on its arguments in
terms of invariants[15,16]. Using the methods of the theory of invariants
6 invariants of a symmetric matrix
and a polar vector
independent from each other have been expressed as follows:
In lieu of the first three invariants of the Green deformation tensor C, the primary invariants given below can be used[16,17].
Incompressibility of the medium, inextensibility of the fiber family and A fiber distribution vector’s being, at the same time, a unit vector make the invariants III, I4 and I5 in expressions (26) and (27) equal to 1, thus canceling Σ’s dependence on these invariants. As a result, the following arguments will remain in the list of invariants Σ depends on.
Taking the derivative of expression (28) according to CPR and substituting it in the Eq. 23 we obtain:
If derivatives of the related invariants involved in this equation according to CPR are extracted from expressions (26) and (27) and substituted, constitutive equation of the elastic stress is obtained in the non-linear form as follows:
A more concrete form of these constitutive equations can be obtained through Lagrange coefficients p and Ta and Σ’s derivatives according to its invariants. It has previously been stated that it is possible to calculate p and Ta from field equations and boundary conditions. To determine Σ’s derivatives according to its invariants, we need to know how Σ, whose dependence invariants are clear from the expression (28), is dependent on these invariants. If Σ is an analytic function of these invariants, it can be represented by a polynomial. In this study, Σ has been represented by a second-degree polynomial according to the invariants it depends on. On the other hand, since internal energy has been defined as positive, for this polynomial to be positively defined and for the sequence of invariants not to affect Σ, it should have a symmetric coefficient. Accordingly, if a polynomial approximation is to be selected, it can be expressed as follows in terms of invariants for the stress potential Σ.
Derivatives of Σ, which appears in the Eq. (30), according
to its invariants can be obtained as follows from the expression (30).
At this stage, disregarding whose function the invariants are, derivatives of the stress potential have been taken according to invariants. Expressions (26) and (27) show us the dependence of invariants in the expression (32). Dependent invariants in the Green deformation tensor (CKL) can be expressed in terms of the strain tensor (EKL), which is a more useful parameter. For this purpose, invariants dependent on the Green deformation tensor, which have been provided in expressions (26) and (27), have been calculated as follows in terms of the strain tensor using expressions generally known in the continuum mechanics.
Terms following the second term on the right side of the Eq.
(30) have been calculated using the partial derivatives provided in the
expressions (32) and invariants provided in the expressions (33). Because mechanical
interactions during these operations have been assumed non-linear in this study,
only terms of the EKL tensor until the second grade have been taken
into consideration. On the other hand, because matrix material is expected not
to be affected by the change of direction along the fiber, expressions of exterior
multiplications in even numbers of the vector field representing the fiber dispersion,
as an argument, have been included into the operations. In this situation the
stress is initially obtained as follows:
without any stress or load (assuming the term α1δPR to be zero), taking the terms into common coefficient parentheses. Coefficients αk(k=1,2,3,...,8) in this equation have been defined as follows and depend on the medium temperature θ.
Equation (34) is the constitutive equation of the stress obtained for the medium considered in scope of this study under the said assumptions in material coordinates. First and second terms of this equation represent contributions of the hydrostatic pressure and fiber stretch to the stress respectively while third and seventh terms, if combined, stand for the linear contribution of the elastic deformation; fourth, fifth and seventh terms-non-linear contribution of the elastic deformation; eighth term-contribution of the fiber dispersion; ninth, tenth and twelfth terms-stress caused by the linear interaction of the fiber dispersion and the elastic deformation; eleventh, thirteenth, fourteenth, fifteenth and sixteenth terms-stress caused by the non-linear interaction of the fiber dispersion and the elastic deformation; seventeenth term-contribution of the non-linear interaction of the powerful fiber dispersion and the elastic deformation. This equation can be reduced to simpler forms for incompressible isotropic media without fiber reinforcement.
SPATIAL FORM OF STRESS
Equation (34) has been substituted in the expression (25), methods generally known in continuum mechanics have been applied, yielding the spatial form of elastic stress as follows:
Here, tensors c-1km and ekl defined as c-1km
and ekl=EKLXK,k XL,l are known as finger deformation
tensor and spatial strain tensor, respectively.
Applying relations generally known in the continuum mechanics on Eq.
(36), the stress dependent on the displacement gradient is obtained as follows[19].
Coefficients in this equation have been defined below:
ANALYSIS OF A HYPERELASTIC MATERIAL UNDER PRESSURE
Strain energy potential (or strain energy function) is represented by W. Strain energy potential can be expressed as follows as a function of either principal stretch ratios or strain invariants.
Various hyperelastic models exist such as Polynomial, Mooney-Rivlin, Neo-Hookean, Yeoh, Arruda-Boyce, Gent, Ogden, Hyperfoam and Blatz-Ko. W strain energy function will be required for parameter types entered as material constants.
| • |
Number of the material constants will depend on the W strain
energy function selected. |
| • |
Despite some certain general rules provided to help the user
select W, selection of W depends on the user. |
| • |
Using the selected W and the material constants entered, ANSYS
software calculates stress and strain behaviour. |
Polynomial form of the strain energy function depends on the first and second strain invariants and can be expressed as follows:
Bulk module and initial cutoff module are defined here as follows:
ANSYS software contains Mooney-Rivlin models with 2, 3, 5 and 9 terms. At the same time these can be viewed in the polynomial form:
| • |
A two-term Mooney-Rivlin model is in the polynomial form when
N=1[19]. |
Mooney-Rivlin model is a suitable model to represent stress and strain behaviour
of almost incompressible natural rubber. For theoretic simplification, hyperelastic
material is assumed to be incompressible. Hyperelastic constants in the strain
energy density function of the material are defined as the mechanical reaction
of the material. Therefore, to obtain correct results during hyperelastic analysis,
Mooney-Rivlin constants of a material should be correctly detected though a
good examination. Mooney-Rivlin constants are generally determined using empirical
stress and strain data. During modeling Mooney-Rivlin constants with two, five
and nine parameters can be used.
During simulation of a hyperelastic material using the finite element method,
geometry of elements before loading can get seriously corrupted. This situation
can lead to the material error. Distortion of mesh formation and structural
instability is not new for researchers and may occur in all traditional materials.
On the other hand, material instability is a new concept for researchers not
dealing with the hyperelastic theory. This concept leads to extra confusion
and difficulties in the analysis of finite elements. In fact, when structural
instability and mesh distortion get combined, the actual cause of deviating
results obtained cannot be clearly explained.
| Table 1: |
Maximum values of materials |
|
|
| Table 2: |
Extremum values of nodal points for different materials |
|
|
|
| Fig. 2: |
Geometry of structure |
After long non-linear studies conducted using ANSYS software researchers are often
likely to get error messages on different solutions. The deviations obtained cannot
be the cause of different practice during the solution operation. Nature of non-linear
analyses and such deviating results in the framework of acceptable error tolerance
levels of the digital solution techniques applied should be carefully considered.
This study analyzes an empty cylinder fixed on one side with pressure applied on the other side. As illustrated on the Fig. 2 and 3, 1000 kPa pressure was applied on the upper surface of cylinder-shaped rubber with empty inside while the other surface was fixed. Inner diameter of the cylinder was 14 cm and the outer diameter 20 cm, while the cylinder was 5 cm deep.
1000 kPa pressure was applied on the upper surface of the hyperelastic material (Fig. 2 and 3). while the other surface was fixed.
RESULTS OF THE ANALYSIS CONDUCTED ON HYPERELASTIC MATERIAL
Mooney-Rivlin model was used for the hyperelastic material with material constants
C1=550 kPa and C2=138 kPa[20]. Poisson ratio
for the isotropic material under consideration has been used at the value of
0.49, which is almost incompressible.
|
| Fig. 3: |
Loads and boundary conditions |
Suitable for the Mooney-Rivlin type modeling, three-dimension HYPER158 element
has been used as the element type. Of hyperelastic material analysis properties,
large deformation impacts have been enabled. Mooney-Rivlin constants were entered
into the data Table. The following results were obtained after the problem was
solved, are expressed in Fig. 4-7.
RESULTS OF THE ANALYSIS CONDUCTED ON FIBER REINFORCED HYPERELASTIC MATERIAL
For a case with similar geometry and hyperelastic properties, 8 fibers were placed in the direction of force applied, after which the analysis was repeated. Properties of the fiber material: an isotropic A36 steel with 200 GPa elasticity module and 0.32 Poisson ratio[22]. Just as in the mathematical analysis, multi-optional way was selected for fibers. Distribution of fibers in the material is as shown in Fig. 8-12.
RESULTS OF ANALYSIS CONDUCTED ON THE A 36 STEEL
To compare the results obtained in the problem another analysis has been made
for an isotropic A36 steel, steel with elasticity module 200 GPa and Poisson
ratio 0.32.
|
| Fig. 4: |
Isometric view of node point displacement in hyperelastic material analysis |
|
| Fig. 5: |
Side view of node point displacement in hyperelastic material analysis |
|
| Fig. 6: |
Isometric view of node point stress in hyperelastic material analysis |
|
| Fig. 7: |
Isometric view of node point strain in hyperelastic material analysis |
|
| Fig. 8: |
Distribution of fibers in the metarial |
RESULTS
In this study, behavior of an isotropic matrix medium of a fiber reinforced
hyperelastic material has been systematically observed against mechanical loads
such medium is exposed to in the external environment, in the framework of continuum
mechanics.
|
| Fig. 9: |
Isometric view of node point displacement in fiber reinforced hyperelastic material analysis |
|
| Fig. 10: |
Side view of node point displacement in fiber reinforced hyperelastic material analysis |
|
| Fig. 11: |
Isometric view of node point stress in fiber reinforced hyperelastic material analysis |
|
| Fig. 12: |
Isometric view of node point strain in fiber reinforced hyperelastic material analysis |
|
| Fig. 13: |
Isometric view of node point displacement in A 36 steel analysis |
|
| Fig. 14: |
Side view of node point displacement in A 36 steel analysis |
|
| Fig. 15: |
Isometric view of node point stress in A 36 steel analysis |
|
| Fig. 16: |
Isometric view of node point strain in A 36 steel analysis |
|
| Fig. 17: |
Selected nodal points on the hyperelastic materials |
|
| Fig. 18: |
Total displacement of between 18 and 98 nodes in hyperelastic material |
|
| Fig. 19: |
Von Mises stress of between 18 and 98 nodes in hyperelastic
material |
|
| Fig. 20: |
Von Mises strain of between 18 and 98 nodes in hyperelastic material |
|
| Fig. 21: |
Total displacement of between 146 and 98 nodes in hyperelastic material |
|
| Fig. 22: |
Von Mises stress of between 146 and 98 nodes in hyperelastic material |
|
| Fig. 23: |
Von Mises strain of between 146 and 98 nodes in hyperelastic material |
|
| Fig. 24: |
Selected nodal points on the Fiber-reinforced materials |
|
| Fig. 25: |
Total displacement of between 26 and 122 nodes in fiber reinforced hyperelastic material |
|
| Fig. 26: |
Von Mises stress of between 26 and 122 nodes in fiber reinforced hyperelastic material |
|
| Fig. 27: |
Von Mises strain of between 26 and 122 nodes in fiber reinforced
hyperelastic material |
|
| Fig. 28: |
Total displacement of between 186 and 122 nodes in fiber reinforced
hyperelastic material |
|
| Fig. 29: |
Von Mises stress of between186 and 122 nodes in fiber reinforced
hyperelastic material |
|
| Fig. 30: |
Von Mises strain of between 186 and 122 nodes in fiber reinforced
hyperelastic material |
|
| Fig. 31: |
Selected nodal points on the steel materials |
|
| Fig. 32: |
Total displacement of between 18 and 98 nodes in A 36 Steel |
|
| Fig. 33: |
Von Mises stress of between18 and 98 nodes in A 36 Steel |
|
| Fig. 34: |
Von Mises strain of between 18 and 98 nodes in A 36 Steel |
|
| Fig. 35: |
Total displacement of between146 and 98 nodes in A 36 Steel |
|
| Fig. 36: |
Von Mises stress of between 146 and 98 nodes in A 36 Steel |
|
| Fig. 37: |
Von Mises strain of between 146 and 98 nodes in A 36 Steel |
Matrix material has been normally assumed to be an isotropic medium,
however, only due to fiber dispersion, it gained the property of a directed
medium, turning into anisotropic looking structure. Because the matrix material under consideration was isotropic, findings of
the theory of invariants have been suitably used to materially determine arguments
of the stress potential. To later obtain the non-linear stress constitutive
equation provided in the expression (30) in a more material way, derivates of
Σ have to be known according to its arguments, which it depends on. Therefore,
stress potential Σ was represented by a second degree polynomial and its
derivatives have been calculated according to its invariants. During these operations
mechanical interactions have been assumed non-linear. Furthermore, considering
that the material was insensitive to direction change along the whole fiber,
an even number of fiber vector components was included in the operation. Later,
constitutive Eq. 34 of stress was expressed in material coordinates
in terms of its components. This expression can be easily induced to simple
mediums where fiberless and mechanical interaction is assumed linear. Form of
the constitutive Eq. 34 in spatial coordinates has been provided
in the expression (37), depending on the displacement gradient.
Some biologic tissues and rubber like solid polymeric materials belong to the
class of fiber reinforced hyperelastic composite materials. Biologic materials
and elastomeric materials used in practical engineering applications are examples
of incompressible composite materials. Non-linear stress behaviour of such materials
can be expressed by the Eq. 37 we obtained.
Numerical solutions have been made for three different types of materials, the first being pure hyperelastic, the second-hyperelastic reinforced by steel fibers (composite material) and the third-pure steel. Besides, while examined area between node points 18-98 and 26-128 have shown axial changes parallel to the external force applied, changes observed between nodal points 146-98 and 186-122 have demonstrated radial changes in vertical direction to the external force applied. Number of node points has been changed for fiber reinforced hyperelastic material analysis because, due to the steel fiber reinforcement, the number of nodes and elements has increased. However, points considered in all three observations occupy the same geometric position. Node points 186, 128 and 26 have been selected as points where two steel fibers were fixed to the hyperelastic material.
Observation of diagrams brings us face to face with interesting and exemplary
situations. For example, comparison of Fig. 19 and Fig.
33 lets us clearly see the behaviour of a rubber like material and a steel
material with high rigidity against a pressure load. General characteristics
of shock absorbers can be seen on Fig. 19. After transferring
the external load to its support, the material shrinks and re-expands itself.
No such behaviour was observed for the steel material with high rigidity. Besides,
Figure 26 provides Von-Mises (average) stress dispersion
in a fiber reinforced elastic material. In this diagram consisting of three
regions the first region is where hyperelastic material was exposed alone. The
second region is where local pressure load was acting on both materials together
(transition region) while the third region is where steel fibers were exposed
alone. Furthermore, at the table presented in the end, values of displacement,
stress and strain have been listed for regions between every node point.
Purpose of the analysis conducted in this study is to consider behaviours of fiber reinforced hyperelastic materials under a pressure load in both analytically and numerically.