On the objectivity of the nonlinear along-fibre shear stress–strain relationship for unidirectionally fibre-reinforced composites

The nonlinear along-fibre shear stress–strain relationship for unidirectionally fibre-reinforced composites has been investigated in this paper aiming at its applications in general 3D stress conditions in a consistent manner. So far, such relationship has only been addressed in plane stress conditions. In this paper, it has been shown that its straightforward generalisation to 3D stress states lacks objectivity, which is a basic requirement for all theoretical studies of physical problems. A new formulation is proposed based on the stress invariants and the complementary strain energy guided by the rational theoretical framework of nonlinear elasticity. A consistent and objective stress–strain relationship has been obtained and verified through an example of application to a torsion problem.


Introduction
Unidirectionally fibre-reinforced composites exhibit significant nonlinearity in their shear stress-strain relationship in the direction along fibres as one of the generic characteristics of the constitutive behaviour of such materials. The phenomenon can be found in textbooks in terms of in-plane shear stress-strain curves, e.g. that by Jones [1]. It is a typical source of material nonlinearity where significant nonlinear stress-strain relationship can be observed within a few percent of the strain. It falls in the conventional small deformation regimes reasonably, at least for engineering applications. The nature of the material nonlinearity was not revealed until the publication of [2] by Mohseni Shakib and Li, where it was found that whilst there was a significant viscosity involved in the deformation, plasticity was a more significant factor, especially at relatively high levels of strains. This was confirmed through the observation of the unloading behaviour after the material had been loaded to different levels of shear strains. The plasticity was further categorised as being of kinematic hardening (as opposed to isotropic hardening) in the terminology of theory of plasticity. A reasonable literature survey to the date of publication could also be found in [2]. The nonlinear shear behaviour was simulated by Li et al. [3] employing a pseudo-plastic model in conjunction with the nonlinearity due to transverse matrix cracking in laminated composites.
A significant characteristic of the along-fibre shear for unidirectionally (UD) fibre-reinforced composites is that it exhibits pronounced nonlinearity within the regime of small deformation, with typical magnitude of strain being up to a few percent. For most engineering applications, the errors stemming from modelling the deformation up to this level of strain within the framework of small deformation are considered to be acceptable. However, if the along-fibre shear stress-strain relationship had been linearly extended over this range of strains, the outcomes of the analysis would be greatly distorted. A typical consequence is that structures could be designed excessively conservative, undermining the benefits of composites for lightweight structural applications.
As early as in 1973, an attempt was made by Hahn and Tsai [4] to represent the nonlinear shear behaviour mathematically, where the strain-stress relationship was derived from the complementary strain energy (CSE). The applicability of CSE should of course not be taken for granted. Whilst the use of CSE will be adopted in this paper, the rational basis of employing it in the formulation of the concerned strain-strain relationship will be presented guided by rational theoretical framework of nonlinear elasticity, in addition to the assumptions of small deformation and static and monotonic loading so that the effects of viscosity and plasticity could be avoided. Their formulation in [4] was narrowed down to the plane stress condition before it was fully developed in its complete form suitable for 3D stresses, as the demand for genuine 3D stress conditions was not high then. As a result, the general applicability of the nonlinear shear stress-strain relationship was not touched there. In fact, to the best of authors' knowledge, there has never been any serious investigation in the literature where the problem of nonlinear shear stress-strain relationship has been addressed under 3D stress conditions. In the meantime, an intuitive and straightforward extension of the 2D relationship into 3D compromises the objectivity of the model as will be demonstrated in the next section of this paper. The principle of objectivity is one of the basic rules in physics. One of its requirements in the context of the present discussion is that the behaviour of the material should not vary when observations are made in different coordinate systems. Basic rules of this kind should be complied with before improvement of the accuracy of any model or theory becomes relevant.
As the composite technology becomes more and more mature, 3D stress states are encountered a lot more frequently than they used to be. Structural parts or components of significant thicknesses are often employed in real applications, where the through-thickness stresses become non-negligible, leading to genuine 3D stress conditions. Another application is in the study of textile composites, 2D and 3D. When the material is analysed at the mesoscopic scale, fibre tows are all subjected to 3D stresses due to the interlacing of the fibre tows. The present paper is partially motivated by this type of applications although the detailed investigation into the problem will be addressed in a subsequent publication.
The nonlinear stress-strain relationships were involved in all three phases of the World Wide Failure Exercises organised by Hinton et al. [5], Kaddour and Hinton [6] and Kaddour et al. [7] where they were provided under plane stress conditions as they were measured. Their applicability to 3D stress conditions has never been seriously investigated by the participants of the Exercises. The objective of the present paper is to bring forward a systematic account on this subject.

Revelation of the lack of objectivity
Consider transversely isotropic materials as a reasonable representation of unidirectionally fibre-reinforced composites, which are the type of composites most widely applied in engineering structures, e.g. as laminae in conventional laminated structures and fibre tows in textile composites. Practically, the available means of obtaining relationship between shear stress and strain are most likely to be under a plane stress condition [2], given as Under the assumption of transverse isotropy, one would intuitively expect In fact, since the material is transversely isotropic, in any coordinate system sharing the common axis-1, i.e. the axis lying in the direction of fibres, one has To facilitate the discussion in this section, the formulation in [4] is adopted for the time being in the following strain-stress relationship: When the along-fibre shear stresses and strains are transformed from one coordinate system to another sharing the common axis-1, one has ⎧ ⎨ ⎩ τ 13 There are two ways of obtaining γ 12 . From the coordinate transformation of strains γ 12 = γ 12 cos θ + γ 13 sin θ = aτ 12 + bτ 3 12 cos θ + aτ 13 + bτ 3 13 sin θ.
Apparently (9) cannot be satisfied in general. In other words, if along-fibre shear stress-strain relationship under 3D stress conditions as defined in (1) and (2) was satisfied in one coordination system, e.g. 1-2-3 system, the same relationship could not be satisfied in general in a different coordinate system, e.g. the 1-2'-3' system sharing the common axis-1. This contradicts the assumption of transverse isotropy of the material. This is a typical sign of lack of objectivity theoretically. Since the coordinate transformations for stresses and strains are equivalent mathematical manipulations as one of the properties of tensors, they cannot be wrong. The only possible cause for contradiction lies in the strainstress relationship as given in (1) and (2) when 3D stresses are involved. For 2D applications, (1) was effectively a curve fitted from experimental data in a straightforward manner and it could not be wrong there. The problem can then be narrowed down to its generalisation from 2D to 3D stress states. Obviously, contradiction (9) results from the nonlinear term in (1) and (2). In other word, with a linear stress-strain relationship, the objectivity is always observed. Any generalisation from 2D stresses to 3D in the presence of nonlinearity in the stress-strain relationship should be made with due respect of objectivity as will be presented in the next section.

Rational considerations
The problem needs to be address rationally. To incorporate the effects of nonlinearity, a nonlinear theoretical framework must be resorted to before linearisation can be made in a consistent manner. In the general problem of Cauchy elasticity, linearisation can only lead to the conventional linear elasticity according to Rajagopal [8]. A more general theoretical framework is implicit theory of elasticity, as was first introduced by Rajagopal [9], which would allow material nonlinearity to be accommodated in the constitutive relationship after the stress and strain involved to be linearised to be suitable for small deformation applications. When Green elasticity, in which stress and strain can be related through some kind of potential, is considered as a special form of implicit theory of elasticity, it could also deliver a nonlinear constitutive relationship between linearised stresses and linearised strains according to Rajagopal [10]. As was in [4], the complementary strain energy (CSE) could be adopted as the potential to derive the strain-stress relationship. However, strict application of Green elasticity results in constitutive relationship with strain and stress both defined in their conventional nonlinear sense. Although the strain and stress can then be linearised to deliver the desired constitutive relationship for small deformation problems, the process of deriving the constative relationship according to Green elasticity where nonlinearity is involved should not be considered to be based on linearised stresses and linearised strains directly. Linearisation to the stresses and strains has to be made after the constitutive relationship has been obtained from CSE.
Another aspect of rational consideration is in the way of constructing the expression of the complementary energy. It will be automatically objective if it is constructed based on the irreducible integrity bases of invariants. For a transversely isotropic body, these invariants as introduced by Mulhern et al. [11] and their nondimensionalised forms are I 1 = S 11 , where stress is defined in its theoretical framework of nonlinear elasticity, and a bar on the top a stress component signifies appropriate nondimensionalisation, for instance, with respect to the Young's modulus of the material, so that the nondimensionalised stresses will be kept small relative to unity to justify subsequent finite truncation of the series.
The nonlinearity in the material is assumed to be elastic, which should be reasonable if the material is subject to static and monotonic loading conditions, as viscosity will not show any effects in a static problem and plasticity can be treated as nonlinear elasticity under monotonic loading. In terms of these invariants, the complementary strain energy can be expressed as where coefficients a 1 , a 2 , b 1 , . . . are material constants. Expression (11) is a polynomial of stresses but involving only independent terms, given the transverse isotropy of the material. Nonlinearity is signified by the presence of terms of orders higher than the second. According to Green elasticity, strain-stress relationship can be obtained as where strain E i j should also be defined according to the theory of nonlinear elasticity in order to incorporate material nonlinearity consistently. Under along-fibre shear alone, only I 4 survives, thus Practically, (11) or (13) can be truncated in order to formulate the strain-stress relationship in a finite form. As it will be argued later, it will be sufficient for the applications concerned to keep the terms up to the order ofĪ 3 4 . Absorbing the constant factors for nondimensionalisation into their respective coefficients to recover the stresses in their original sense, one has where the involved material constants have been denoted as a, b and c for simplicity after absorbing the respective nondimensionalisation factors and constant factors associated with the orders as opposed to d 1 , d 2 and d 3 in (13) so that the derivatives of U would take their simplest appearances. An argument of nondimensionalisation should also be understood as being present in the background of similar accounts employing stress invariants, e.g. that for the construction of failure criteria for composites by Hashin [12], although it was not presented explicitly. The corresponding along-fibre shear strains can be obtained as where the shear components are presented in the engineering sense. Only after (12) or (15) has been obtained, linearisation can then be made to S i j and E i j so that One can now obtain the following strain-stress relationship incorporating the material nonlinearity whilst having the stress and strain linearised to be suitable for small deformation applications: Apparently, the above expressions are objective if they are re-written into where K = a + b τ 2 12 + τ 2 13 + c τ 2 12 + τ 2 13 2 is coordinate transformation invariant between coordinate systems sharing common axis-1, bearing in mind that τ 2 12 + τ 2 13 = I 4 and therefore stresses on the right-hand side of (18) are transformed in exactly the same way as the strains on the left-hand side, given (5) and (6). Thus, any coordinate system, i.e. 1-2'-3', sharing the common axis-1, one has γ 12 = aτ 12 + bτ 12 τ 2 12 + τ 2 13 + cτ 12 τ 2 12 + τ 2 13 2 , γ 13 = aτ 13 + bτ 13 τ 2 12 + τ 2 13 + cτ 13 τ 2 12 + τ 2 The along-fibre shear strain-stress relationship as given in (17) is therefore objective. The two along-fibre shear stresses involved are coupled in general under a 3D stress condition. There are three material constants, a, b and c, to be determined in order to describe the nonlinear along-fibre shear strain-stress relationship. Since they are material constants, their values should not be affected by stresses applied. Therefore, a plane stress state can be employed to determine them. Under this stress condition, one has which is in line with [4] except that the order of the polynomial has been raised to the 5th. Consider IM7/8551-7 composites whose nonlinear shear strain-stress data could be found in [13,14] by Kaddour et al. as shown in Fig. 1 as the diamond symbols. The data in Fig. 1 have been fitted into two polynomial functions. For one of them, the polynomial has been truncated at τ 3 12 , i.e. the same order as in [4]. Apparently, a significant discrepancy can be observed. In the second case, the polynomial was of the 5th order, i.e. as given in (20), and its agreement with the experimental data is satisfactory. This was the reason why the I 3 4 term has been kept in the truncated complementary strain energy expression (14).
with the unit for stress being MPa. They have thus been determined from available experimental data. Consider the plane stress condition as a special case of 3D stress state, these constants can be readily incorporated in (17) to deliver the proper along-fibre shear strain-stress relationship in general. As a 'sanity check', the inverse of constant a provides the initial slope of the stress-strain curve, i.e. the linear shear modulus. This gives a value of 5.35 GPa. The first 7 experimental stress-strain data points as originally provided in [6] followed a straight line of a slope, i.e. linear modulus, of 5.56 GPa which was rounded up to 5.6 GPa and provided in [13,14] as the linear elastic shear modulus. The disparity between 5.35 and 5.56 GPa offers an indication of the discrepancy between the original data and the fitted 5th-order polynomial curve.
There are a number of considerations associated with the along-fibre shear strain-stress relationship obtained above: (1) In addition to the justification given in [4] regarding the choice of complementary strain energy instead of strain energy, there can be another reason as stated below. The value of a polynomial will always approach infinity towards both extremes. Take the positive side as an example for sake of argument, the derivatives share the same sense as the value of the function since they are all dictated by the sense of the highest order term in the polynomial. The distributions of available experimental along-fibre shear stress-strain data for a range of known materials [13,14] tend to show a trend in which as the stress and strain increase monotonically but the rate of increase in stress reduces, i.e. the first-order derivative being positive whilst the second-order derivation being negative. A straightforward polynomial is therefore not ideal in order to capture this trend. However, if the experimental data are viewed as the strain-stress relationship, as shown in Fig. 1, then the above-described features of polynomials are precisely what one would like to have, i.e. the rate of change in strain increases as stress and strain increase. Relatively, the rate of change in stress reduces. Therefore, the constitutive relationship will be given as strain versus stress.
(2) The order of the polynomial has been raised from 3rd as in [4] to the 5th. This is to obtain better fit to the experimental data as shown in Fig. 1 whilst backed with rational considerations. An underlying assumption is that the along-fibre shear behaviour is decoupled from the other stresses. There could be mechanisms to incorporate interactions between the along-fibre shear stresses and other stresses, such as one outlined in item 3 below. However, designated experiments will have to be conducted in order to generate data to characterise such behaviour. A fair assessment is that, as the state-of-the-art of the development, the composites community is perhaps not ready for that yet. (3) It can be noted that if one includes the term I 4 I 5 , it is possible to reflect the different shear stress-strain responses in the presence of transverse stresses of different senses. However, necessary experimental data are required in order to incorporate this effect. In absence of such data, this will not be pursued in this paper. (4) The experimental data in [13,14] was provided in the first quadrant only. When fitting them into a polynomial directly, it is likely to result in even terms in the polynomial which contradicts the odd nature of shear stressstrain relationship. The curve fitting employed was based on antisymmetrically extended experimental data, as shown in Fig. 1 to avoid undesirable even terms in the fitted function.

Obtaining stresses from strains
In the actual applications, it is often desirable to express stresses in terms of strains. Another rational point needs to be noted here if one intends to obtain an inverted relationship of (17). According to [8], the inverted relationship should not be obtained from the strain-stress relationship (17) directly with the stresses and strains having been linearised already. One has to recover their appropriate form according to their strict definition in nonlinear elasticity before the inversion takes place. Linearisation of stresses and strains should proceed afterwards because linearisation and inversion as two mathematical operations are not commutative, in general [8]. Practically, because it is impossible to invert (17) in general as a simultaneous set of high-order equations cannot be solved analytically, explicit inversion will not be pursued here. As an alternative, numerical solution will be sought as follows by employing the Newton-Raphson iterative approach. The shear stresses are updated through the iterations as follows: taken from the beginning of the current increment, assuming that an incremental and iterative procedure has been followed. The iteration continues until convergence is met which can be signified by ⎧ ⎨ ⎩ τ 12 where ε is a designated tolerance, taken as 10 −6 , for the examples to be shown later in this paper. Given (17), the Jacobean matrix involved in (22) where the values of stresses should be taken from those of the previous iteration in order to facilitate the current iteration. The iterations usually converge rapidly. It should be noted that the inverse of matrix (24) as involved in (22) also gives the tangential stiffness matrix for the along-fibre shear behaviour upon convergence of the iterations.

Verification, an example of application and discussion
Using the Newton-Raphson iteration scheme, the stress-strain relationship has been obtained effectively and given a set of along-fibre shear strains, their corresponding shear stresses can be obtained. The relationship applies in any coordinate system sharing the common axis-1 as appropriate for transversely isotropic materials with the same set of the material constants a, b and c.
The newly established strain-stress relationship (17), applicable to 3D stress states, obviously reproduces its 2D stress state counterpart given in (20) when the other along-fibre shear stress vanishes. It has also been verified that the implicit stress-strain relationship ⎠ as the inverse of (17) obtained through the iterative approach given in (22) reproduces the same of (20) under a plane stress condition. In other words, the 3D stress-strain relationship reproduces that in 2D stress state when the stress state is 2D. An interesting observation out of this exercise is presented as follows. When the material is under τ 12 alone, i.e. in 2D stress state, the tangential stiffness associated with this stress component reduces as the stress increases, as expected. As a result of the interactions between the stresses in the strain-stress relationship (17), the tangential stiffness associated with τ 13 also reduces in the absence of τ 13 , although at a different rate from the one associated with τ 12 . This has been shown in Fig. 2 where the stress-strain curve is the same as that in Fig. 1 in the 1st quadrant but with stress and strain swapped around. The local increase in both tangential stiffnesses at small values of strain (less than 0.01) is not genuine material behaviour. Over this strain range, the respective tangential stiffnesses should remain constant as the material behaviour is still within the linear elastic regime. The deviation from a straight horizontal line can be attributed to the nature of polynomial functions, in general, that tends to leave more pronounced errors in their derivatives. The variation of the tangential stiffness associated with τ 13 when the material is under τ 12 alone as shown in Fig. 2 can also be justified by examining (17), where the compliance associated with τ 13 when the material is under τ 12 alone apparently increases with τ 12 . The corresponding stiffness as the inverse of compliance should therefore reduce.
In order to employ the nonlinear shear stress-strain relationships established in the previous section in structural analyses, it has been implemented into a user defined material subroutine UMAT so that it can be incorporated in a commercial FE code, Abaqus [15]. The UMAT has been verified first using a single element under various along-fibre shear stress conditions, resulting in correct outcomes as expected.
As an example of applications, the torsion of a solid UD composite cylinder is studied with the results presented below. This example has been carefully chosen. It is simple enough in the sense that it could essentially be solved analytically with a little assistance of numerical integration. The dimensions of the cylinder were assumed to be D = 10 mm for diameter and L = 50 mm for length. The length is irrelevant as far as this analytical solution is concerned, because the deformation is described as angle of twist per unit length, but it is required for creating FE model as involved subsequently. To complete the material definition in the FE model, elastic properties E 1 = 168 GPa, E 2 = 8.4 GPa, ν 12 = 0.34 and ν 23 = 0.50 from [13,14] have been employed, in addition to a, b and c as given in (21). The mesh will be shown when the contour plots are presented later. The types of elements employed are C3D20R and C3D15R [15] which are both quadratic solid elements with reduced integrations.
The simplicity of the problem which allows an analytical solution comes from the fact that the stress state in a cylindrical coordinate system is plane stress, i.e. the only non-vanishing along-fibre shear stress is in the circumferential direction whilst the other along-fibre shear stress in the radial direction vanishes identically. In this case, (20) can be used perfectly correctly if direction-2 is understood to be the circumferential direction. The analytical solution is based on the mechanics of materials model, where plane section assumption has been adopted, and the shear strain is assumed to vary linearly with the radial coordinate, but the shear stress is obtained from (20) by solving it using the Newton iterations. The torque is obtained by integrating the shear stress as follows: where θ is the angle of twist per unit length of the cylinder and r is the radial coordinate which varies from 0 at the centre of the cylinder to D/2 on the surface of it. The numerical integration as indicated in (25) can be conducted with the radius being divided into n equal intervals with n being a sufficiently large integer and γ i = r i θ = i Dθ/2n with i = 1, 2, . . ., n. Expression τ (γ ) is the stress-strain relationship as an implicit function defined by (20). Given its implicit nature, integration can be conducted in a numerical manner and the results are shown in Fig. 3. This torsion case is also substantial and meaningful enough if it is presented in a rectangular coordinate system. In this case, the ratio between τ 12 and τ 13 varies from -∞ to +∞ depending on the location over the cross-section and is therefore an ideal case to verify the proposed model. This case will be too complicated to solve manually and FE analysis has been adopted. Eight sets of results altogether have been obtained, respectively, in order to draw informative comparisons: (1) Analytical results are readily available from mechanics of materials, e.g. the textbook by Gere and Timoshenko [16], based on linear elastic stress-strain relationship. They serve as a relevant representation of the nonlinear behaviour when the magnitude of deformation is sufficiently small. The result is shown in Fig. 3 and it is a straight line as expected.
(2) Results obtained with the analytical approach involving nonlinear shear as described above through (25) and (26). (3) FE analysis with linear stress-strain relationship in a rectangular coordinate system. (4) FE analysis with linear stress-strain relationship in a cylindrical coordinate system. (5) FE analysis with the nonlinear model as established in this paper using FE in rectangular coordinate system.
The contour plot of the circumferential along-fibre shear stress at the 0.012 rad/mm angle of twist per unit length has shown in Fig. 4a, where the along-fibre shear stresses have been transformed according to equation (5). An axisymmetric distribution of the stress can be seen clearly in Fig. 4a as it should be. The radial along-fibre shear stress has also been plotted as given in Fig. 5a at the same angle of twist per unit length. As the legend indicates, this stress is 6 orders of magnitude lower than its circumferential counterpart and therefore is negligible, again as expected. What has been shown in the contour plot is nothing but numerical noise of rounding errors.  Torque versus angle of twist per unit length (8 curves: referring to the legend, the 1st, 3rd and 4th are linear predictions for which curves coincide with each other; the 5th, 6th and 8th are nonlinear FE predictions and also coincide; the 2nd is also a nonlinear prediction but by analytical means, and the 7th is nonlinear FE prediction but lacks objectivity and hence standing out) (6) FE analysis with the nonlinear model as established in this paper using FE in a cylindrical coordinate system. The result is identical to that from the rectangular coordinate system as expected. This verifies the objective nature of the established model. (7) The result based on the following uncoupled strain-stress relationship: γ 12 = aτ 12 + bτ 3 12 + cτ 5 12 , γ 13 = aτ 13 + bτ 3 13 + cτ 5 13 .
It is not expected to be correct as it is not objective, but the purpose is to demonstrate the consequence of the lack of objectivity. The error in the torque obtained can be seen in Fig. 3 as a single case standing out of the cluster of other predictions. The contour plot of the circumferential along-fibre shear stress as shown in Fig. 4b is also at the 0.012 rad/mm angle of twist per unit length so that it can be compared with Fig. 4a directly. Apart from the increased stress level corresponding to the increase in torque as shown in Fig. 3, it is not difficult to notice that the axisymmetric appearance is no longer present, which is a clear sign of the dependence on the coordinate system selected. It is even more revealing if the radial along-fibre shear stress is plotted as shown in Fig. 5b where non-trivial values have been predicted, which should not be present in a correct solution based on common sense as well as strict formulation. The distribution of this stress appears to be systematic. The fact that its magnitude diminishes towards the surface of the cylinder is not because of the constitutive model becomes reasonable there, but because of the boundary conditions and the equilibrium requirements as imposed by the FE model. (8) The final one is to demonstrate the sign of lack of objectivity of (27) in a different way. If (27) is used along with the cylindrical coordinate system for the current torsion problem, it can in fact lead to correct result but only for the present problem where this is only one non-vanishing along-fibre shear stress in the cylindrical coordinate system. Changing the problem slightly, e.g. if the cross-section of the bar deviates from being circular, when both along-fibre shear stresses will show up, (27) will not be able to produce correct result anymore.
Through the results from the torsion problem shown above, clear conclusions can be drawn as will be presented in the next section.

Conclusions
The following conclusions have been established through the formulation, elaboration and verifications as presented in this paper, based on the assumptions of small deformation and monotonic loading.
(1) The intuitive strain-stress relationship, such as that given in (27), lacks objectivity and hence is wrong. It is not applicable to problems involving 3D stresses in general. (2) Strain-stress relationship (17) as established in this paper is consistent and objective backed by rational considerations, including nonlinear theory of elasticity and its appropriate linearisation, and stress invariants for transversely isotropic materials in terms of the irreducible integrity bases for the construction of the complementary strain energy. (3) Strain-stress relationship (17) is applicable to problems involving monotonic loading conditions in a static sense. Viscosity and plasticity may have to be considered otherwise. It remains a challenge to formulate the constitutive relationships under these conditions.
The applications of the established nonlinear along-fibre shear model in more practical problems, in particular, for fibre tows in modern textile composites, will be addressed along with other material and structural complications in a subsequent publication. Contour plots for the radial along-fibre shear stress based on a the objective strain-stress relationship as given in (17) and b the non-objective one as given in (27)