Derivation of Ply Specific Stiffness Parameters of Fiber Reinforced Polymer Laminates via Inverse Solution of Classical Laminate Theory

Period. Polytech. Mech. Eng. L. Kovács, G. Romhány Abstract The realistic estimation of the ply stiffness parameters of polymer composite laminates is a big challenge nowadays in industrial practice. In this paper a new, innovative concept is introduced that is based on the backward use of Classical Laminate Theory (CLT). The innovation in this new concept is (amongst others): possibility to infer the stiffness constants from the simple mechanical tests of specimens with multidirectional ply stack-up identical to the part to design. In addition the new method is manifested in a form of a compact equation that surely returns the measured deformation of the tested specimen on laminate level. The mathematical background of this concept is slightly more complex than what the conventional techniques offer, however its explicit form allows to code it in any automatic systems (e.g. user script) that can be run in Finite Element environment or as part of the software of a mechanical testing frame.


Introduction
The (mainly Finite Element Method -FEM) based design and analysis of composite laminates consisting of fiber reinforced polymer (FRP) plies requires to know the stiffness characteristics of each individual ply.Using these the deformation and thus, the stress state of the part can be estimated that helps to assess it from strength point of view.
In engineering practice the ply specific stiffness is approximated by a bespoke orientation dependent material model (e.g.orthotropic, transversely isotropic or isotropic) and from the ply behavior the stiffness and deformation of the whole laminate can be forecast by using an appropriate plate theory (e.g.Kirchhoff plate theory [1,2], Reissner-Mindlin theory [1] etc.).The hierarchy between the individual building blocks of the laminate level material model is shown in Fig. 1.

Fig 1 Estimation process of composite laminate stiffness
As it follows from Fig. 1 the literature and industry knows numerous ways to estimate the stiffness of a composite laminate.There are techniques predominantly used in the scientific world that attempt to handle the problem in micromechanical as well as molecular scale.The molecular dynamics (MD) methods provide promising results in researching the mechanical behavior of non-oriented nanocomposites.However, these approaches can mainly be exploited when the effect of nano-sized particles (e.g.carbon nanotubes -CNT [3][4][5] or graphene nanosheets GNS [5]) on the base matrix material is quantified, therefore these are useful when relative trends need to be understood.
These techniques have only relevance from composite standpoint if the micromechanical or molecular models of the individual constituents are linked with an appropriate micromechanical FEM or analytical rule of mixture.In this case from the characteristics of the independent constituents the behavior of the ply can be inferred.
When micromechanical FEM is being used, a micro scale Representative Volume Element (RVE) of the ply is built such a way that using the mechanical properties of the individual particles those are approximated by appropriate finite elements and the interactions between them are defined by realistic constrain equations [3,6].A somewhat more simple approach is to estimate the stiffness of the composite ply from the constituents' properties with the aid of the rule of mixture (ROM).Such widely used models are for instance the Halpin-Tsai and the Voigt-Reuss [7] models.The research works described in [3,6] and [8] are good examples for this.
The approaches mentioned above due to their complexity need validation via mechanical testing of appropriate specimens.This is partly the reason why in the industry the stiffness is directly derived from ply mechanical testing instead of implementing any sophisticated micromechanical approaches.These measurements are usually well described and easily implementable processes covered by a wide range of international standards that are hard-coded in the control software of testing frames.The most common mechanical testing types are collected in Table 1 along with the relevant standards.
In addition, there is another mechanical measurement technique the so called Dynamic Mechanical Analysis (DMA) which also provides a complete view on the static and dynamic stiffness (and damping) of a single ply or a whole laminate even as a function of temperature.About the practical use of DMA [9] yields a comprehensive review.
Also in scientific world the aforementioned mechanical tests are widely utilized to investigate inhomogeneous composite ply stackups [10,11], short fiber reinforced composites [12] or even nanocomposites [13].

Comparison of conventional and new concept
The outcome of standard mechanical tests is recorded load vs. strain curves (in predefined global directions).Processing these curves leads to the stiffness parameters sought.Thus, these material constants can be calculated with relatively simple equations laid down in the relevant standards.Table 1 provides a non exhaustive list of typical testing types and the in-plane stiffness constants that can be evaluated from them.
The constants derived from the experiments listed above represent the behavior of fully uniform UD laminates with homogenous fiber orientation (or in some cases laminates with [0/90] symmetric and balanced layup).In practice, these parameters are then inputted to the calculations estimating the deformation and stresses of a composite part consisting of a multidirectional ply stackup.These calculations use FEM with the aid of an appropriately selected plate theory (e.g.CLT) as it is shown in Fig 1 .This procedure is called as the "conventional way".
The simplicity of the conventional way is counteracted by a number of drawbacks such as: • the formulas included in the standards are not applicable to infer the stiffness parameters for specimens having special multidirectional ply stackup.However, in practice the structures and components do have such ply structure.different test types usually mean different loading manner and hardware with dissimilar accuracy and quality, therefore the scatter (distribution) of the inferred constants varies and this does not necessarily come from the natural "built-in" uncertainty of the material.Wrongly coupled stiffness constants can lead to wrong deformation predictions for the part or component to design.This concern is especially high when the test data exhibit a significant scatter.The best example for this is the experimental evaluation of in-plane shear modulus (G 12 ).The different loading manners listed in Table 1 can provide significant differences in the calculated G 12 even if the specimens are extracted from exactly the same base laminate.
To remedy the issues collected above, as part of an R&D project called NVKP-16-1-2016-0046 [14] launched within the scope of the National Competitiveness and Excellence Programme (NCEP) and supported by the National Research, Development and Innovation Fund (NKFIA) a new, innovative concept was elaborated that approaches the problem of inferring material constants from experiment from a different standpoint.This new concept is more efficient and can provide more accurate solution than the conventional way for the following reasons: • it is capable of back-calculating the in-plane stiffness constants from specimen testing with almost any ply stackup sequence, • when the ply stackup is selected wisely it can return all four parameters in one go from one single test.• If all specimen tests are performed on the same frame in the same manner, the statistical concern about the conformity of the parameters mentioned earlier does not exist and therefore the scatter (distribution) of the derived constants will much more likely cover the uncertainty purely induced by the material itself.

• If the specimens to test have the same ply stackup
sequence as what the finished product has, the new concept ensures that the inferred stiffness constants will more accurately estimate the deformation of the product in operation than the conventional way.This is because the stiffness parameters do not need to be extrapolated to the laminate since the new method actually starts from the experienced deformation at that level.• It follows from the statements above that the new concept provides more accurate results than the conventional one even if the stackup sequences between the test pieces and the finished product differ.
The importance of the new concept is even bigger when it comes to computing the stiffness of multidirectional laminates reinforced by non-crimped fabric (NCF).The biggest advantage of NCF composites is that the ply stackup sequence of those can be flexibly set in the design phase of the actual product in order to ensure the desired deformation of the structure in operation.This application gradually becomes more important as the "morphing" composites get more widespread (e.g.controlled operational twisting of an airplane wing)."Morphing" composites have necessarily asymmetric and multidirectional layup.For such layups the use of the conventional method is rather questionable.
The new concept shown here utilizes the classical laminate theory (CLT) by solving it backwards in order to find the basic ply specific in-plane stiffness constants.The input data to it are the deformations and external loads taken from the mechanical tests conducted on plates with arbitrary layup.Since this technique is in fact an inverse solution of CLT it was briefly named as CLT -1 .

The CLT -1 method
The basic assumptions that need to be adhered to in order to make the concept work are: • all plies of the laminate have the same fiber-resin structure (e.g.all UDs or all woven fabric reinforced), • the composite material can be reasonably approximated with an orthotropic deformation model on ply level, • the engineering constants describing one layer are constants through the thickness of one ply, • the kinematics of the entire composite plate is approximated with the Kirchhoff plate theory [2], • the stiffness of the entire composite plate is not influenced by external restrains (i.e. it can deform freely due to the external forces), or the reaction forces and moments acting at the restrains are known.
The method is initiated from the basic system of equations of CLT (including thermal expansion) [1]: Notations to the equations above: t k -thickness of the ply with index "k" [mm], ( (1) z k -coordinate of the middle plane in thickness direction of the ply with index "k" measured from the laminate middle plane as zero level [mm].The structure of the laminate along with the notation of the geometrical parameters measured from the middle plane of the laminate are highlighted in Fig. 2.  4) equations into it leads to the following system of equations: Similarly, splitting the matrix equation in (2) into individual equations as well as substituting the (4)-( 5) equations into it leads to the following system of equations: In Eqs. ( 6)- (11) the formulation for the introduced extra variables is as follows: The extra variables detailed in Eqs. ( 12)-( 17) need to be evaluated for each ply in the laminate once the ply specific thermal strains have been evaluated.The thermal strains in a given direction are to be calculated as the product of the thermal expansion coefficient (CTE) in the given global direction and the temperature difference.The CTE values in the local ply specific material coordinate system can be evaluated from designated experiments.Finally, the CTE vector of [α 11 , α 22 , α 12 ] in the local coordinate system has to be modified with a rotation transformation to represent the CTEs in global coordinate system.
The components of the ply specific stiffness matrices considered in the global coordinate system (denoted as Q ij k     in Eqs. ( 6)-( 11)) can be expressed with the ply stiffness parameters representing the actual local material coordinate system.To do this, the rotation transformation of the local stiffness matrix is required that is expressed by the following matrix equation [15]: where [T σ (Θ k )] and [T ε (Θ k )] are the rotation transformation matrices of the stress and strain vectors of the given ply respectively and Θ k is the orientation of ply with index "k" compared to the global x axis.
[Q] k is the stiffness matrix of the ply with index "k".The structure of the transformation matrices is detailed in [15].After expressing the Q ij k     components as a function of the stiffness parameters in the local ply specific material coordinate system (Q ij ) with the aid of Eq. ( 18) and substituting those expressions into Eqs.( 6)-( 11) leads to a system of equations that can be transformed into matrix form: In Eq. ( 19) the matrix denoted as [T] contains all information about the layup structure and measured deformation of the laminate in question, therefore it has named as Stackup Specific Deformation (SSD) matrix.The 4 element vector denoted as [Q] in this equation collects all in-plane stiffness parameters of a ply in the local material coordinate system that are supposed to be the same as per the initial assumptions.Without giving the details of the derivation the elements of SSD matrix are to be calculated as follows: In the equations above the following ply orientation dependent variables were introduced: ) ( ) The linear system of equations in (19) can then be solved in one step by inverting the SSD matrix.The outcome of this will be the [Q] vector that in fact contains all four independent stiffness parameters of the ply.
Lastly the individual stiffness parameters can be evaluated from [Q] vector like this: ).

The SSD matrix
In classical sense from mathematical point of view inverting the SSD matrix is not possible since it is not a symmetric nXn type matrix (by having 4 columns and 6 lines).However, there is a technique called "pseudo-inversion" that can return a matrix which provides an nXn dimension unity matrix if the original gets multiplied with it from the left or an mXm dimension unity matrix if from the right.
Only one of the aforementioned two cases can work at a time when it comes to pseudo-inverting a non-symmetric matrix.When the matrix to invert is a full column rank matrix (which is the case here for a 6X4 dimension SSD matrix with no full zero rows or columns) the formulation to find the pseudo-inverse is the following: where [T] T denotes transposing.The conditions for a matrix to be pseudo-invertible are defined by the Moore-Penrose law.When pseudo-invertibility is possible, building SSD matrix and the vector of external loads can be easily automated when the deformations of the laminate are known.From that, as it is described above, obtaining the in-plane stiffness data is one additional step only (see (50)).If a large number of tests are performed, the method described here needs to be repeated that gives a data pool for all derived constants.This enables to do rigorous statistical assessment on the parameters as well that provides reliable information on the variability of the parameters in question.
sequence of the specimen to test wisely it is possible to obtain all four in-plane stiffness parameters from one single experiment.This new method due to its compactness can deliver more reliable results than the UD test based conventional methods since it is able to find the constants for test specimens having a ply stackup identical to the product to design and it does not require multiple different type tests with different fidelity to derive all four constants.With this method the questionably accurate in-plane shear tests can be replaced by simple tensile or flexural experiment of test plates with appropriately chosen stackup sequence.The new procedure largely increases the fidelity of the inferred stiffness data for NCF composites that will likely be the base material for morphing structures of the future.
A key element of this new method is the so called Stackup Specific Deformation (SSD) matrix.The analysis of SSD matrix and in light of that the assessment of the practical usability of this concept along with the validation with real experiments will take place in 2018 as part of the NVKP-16-1-2016-0046 [14] project.
ij -element of the laminate stiffness matrix for tension, B ij -element of the tension-flexure coupling matrix, D ij -element of the laminate stiffness matrix for bending.N -number of plies building the laminate, Q ij k     -element of the stiffness matrix of the ply with index "k" in the global coordinate system [MPa], ε xx tot-0 , ε yy tot-0 , ε xy tot-0 -membrane strain components in the global coordinate system acting in the middle plane of the laminate, κ x , κ y , κ xy -curvature components of the middle plane of the laminate in the global coordinate system [1/mm], [ε xx T ] k , [ε yy T ] k , [ε xy T ] k , -thermal strain components of the ply with index "k" in the global coordinate system, N x , N y , N xy -external force components per unit length (in the width direction of the laminate) in the global coordinate system [N/mm], M x , M y , M xy -external moment components per unit length (in the width direction of the laminate) in the global coordinate system [Nmm/mm].

Fig 2
Fig 2 Parameters describing the geometry of a composite laminate

Table 1
Typical standard mechanical tests to find in-plane stiffness parameters ν 12 and ν 21 are the major and minor Poisson ratios respectively