A two-way loose coupling procedure for investigating the buckling and damage behaviour of stiffened composite panels

Accepted Manuscript A two-way loose coupling procedure for investigating the buckling and damage behaviour of stiffened composite panels Sina Hühne, José Reinoso, Eelco Jansen, Raimund Rolfes PII: DOI: Reference:

S0263-8223(15)00916-2 http://dx.doi.org/10.1016/j.compstruct.2015.09.056 COST 6901

To appear in:

Composite Structures

Please cite this article as: Hühne, S., Reinoso, J., Jansen, E., Rolfes, R., A two-way loose coupling procedure for investigating the buckling and damage behaviour of stiffened composite panels, Composite Structures (2015), doi: http://dx.doi.org/10.1016/j.compstruct.2015.09.056

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A two-way loose coupling procedure for investigating the buckling and damage behaviour of stiffened composite panels Sina H¨ uhnea , Jos´e Reinosoa,b , Eelco Jansena,∗, Raimund Rolfesa b Group

a Institute of Structural Analysis, Leibniz Universit¨ at Hannover, Appelstr. 9A, 30167 Hannover, Germany of Elasticity and Strength of Materials, School of Engineering, University of Seville, Camino de los Descubrimientos s/n, 41092, Seville, Spain

Abstract This paper is concerned with the development of a novel FE-based two-way loose coupling approach for the analysis of composite stiffened panels. The aim of this numerical strategy is to investigate the global postbuckling behaviour as well as the local damage progression of composite structures using separated FE models with different levels of fidelity: (i) a relatively simple global model of the complete structure, and (ii) more complex local models of certain details that incorporate damage capabilities to simulate damage events. In the coupling process, information is exchanged between these diverse models to simulate the overall structural behaviour including geometrical as well as material nonlinearities. The two-way loose coupling character of the methodology allows, firstly, a direct interaction between the global and the local levels along the solution process and, secondly, a highly versatile adaption with regard to the definition of the local models. In addition, the separation of the models and analyses in the approach enables the use of standard FE software without involved implementations or modifications of the source code. The developed coupling procedure is assessed through two applications: (i) an academic composite stiffened panel, and (ii) a real stiffened panel taken from literature. The results of the proposed coupling approach are compared with the numerical and experimental reference data, exhibiting a satisfactory level of accuracy. Keywords: Loose coupling, Multiscale analysis, Composite structures, Stiffened panels, Damage mechanics, Buckling and postbuckling

1. Introduction The use of fibre-reinforced composites for lightweight constructions as aircraft structures or rotor blades of wind turbines has been gradually increased due to their excellent material properties, such as high specific stiffness and strength ratios. Adversely, the application of composites is related to relatively high production and material costs, which makes the design of such structures with a minimised structural weight a very important aspect. However, in the context of aviation and aerospace applications weight saving is in general a crucial topic, especially regarding cost reduction and fuel consumption. At present, most designs of composite structures are still rather conservative, since the overall structural response, including buckling and material damage of such structures, cannot be predicted in a reliable and efficient way because of the complex failure behaviour of composites, that involves different damage mechanisms. 1.1. Coupling methods With the aim of improving the level of reliability regarding the simulation of composite structures, the applicability of global/local and multiscale methods has been widely explored in the last two decades. These techniques connect several FE models of different size and/or fidelity throughout the computations. ∗ Corresponding

author. E-mail: [email protected], Tel.: +49 511 762 17425.

Preprint submitted to Composite Structures

October 6, 2015

From the computational perspective, the connection between the different models that compose the solution scheme can be carried out in various ways, leading to many different categories of global/local and multiscale methodologies. Possible distinctions among the existing methods regard one-way and two-way techniques as well as loose and tight coupling procedures. On the one hand, one-way compared to two-way coupling means that the information transfer between the different models is only conducted in one direction, i.e. from the global to the local model or from the local to the global model. On the other hand, two-way coupling techniques obviously transfer information in both directions. With regard to the interrelation between both levels of analysis, loose coupling implies the use of spatially separated models or at least of distinct simulations that exchange information but which are computed consecutively. Contrarily, the term tight coupling describes numerical techniques that solve the global and local systems of equations simultaneously. In this setting, from the authors’ perspective, the fundamental difference between the global/local and the multiscale denominations refers to the existence of scale separation. A scale separation is present if the characteristic lengths at the considered levels are significantly different, so that the strains and stresses of the local level can be regarded as constant on the global level. In the global/local strategies, usually no specific scale separation is assumed among the different models that compose the procedure, whereas in the multiscale approaches, a scale separation is compulsory and allows the use of the Representative Volume Element (RVE) concept. Specifically, multiscale techniques rely on the so-called classical homogenisation procedures, in which the mechanical properties at the large scale are determined through averaging the small scale behaviour using the aforementioned RVE. 1.2. Tight coupling Some typical two-way tight coupling approaches without scale separation deal with transition elements, multipoint constraints or other numerical features to insert different meshes into one combined model. Aside, most of the so-called coupled numerical homogenisation techniques involve individual models of separated scales and, therefore, belong to the group of multiscale approaches. Within this context established methods are the landmark contribution about the FE2 method of Feyel (2003), the second-order homogenisation by Kouznetsova et al. (2004) and other homogenisation-based approaches, e.g. by Terada and Kikuchi (2001), Fish et al. (1999), Ghosh et al. (2001), Miehe and Koch (2002), to quote a few of them. Hierarchical methods, as the variational multiscale method by Hughes et al. (1998), the bridging multiscale method by Kadowaki and Liu (2004) or the coupled volume approach by Gitman (2006) are instead suitable for cases without specific scale separation. In these latter approaches, the global solution is enriched with information from the local solution including local nonlinearities, via definition of different constitutive laws or interactions along with mesh and geometrical refinements, which are generally omitted in the global model. Since a tight connection always implies an interaction between the models involved in both directions, one-way tight coupling approaches do not exist. 1.3. Loose coupling The area of one-way loose coupling comprises many methods for the direction from the global to the local level, such as the global-local methods developed by Mote (1971), Noor (1986), Mao and Sun (1991), or more specific zooming techniques by Rank (1993) or submodelling approaches as by Reinoso et al. (2012), among other alternatives. For the direction from the local to the global level the uncoupled homogenisation by Zohdi et al. (1998), substructuring methods or the multiscale failure analysis by Ernst et al. (2010) are available among other employed strategies. Instead, the field of two-way loose coupling implies only few methods, so that this topic represents an important area of potential research. Mentionable contributions are the multiscale projection method by Loehnert and Belytschko (2007), the non-intrusive global-local technique by Gendre et al. (2009), the adaptive progressive damage modelling technique by Labeas et al. (2012) and the homogenisation-based iterative multiscale approach (HIMSA) by Chrupalla et al. (2011).

2

1.4. Existing approaches However, the different multiscale and global/local approaches described in the literature can have different drawbacks and/or limitations. For tight coupling approaches with embedded models, the critical areas need to be known in advance, since the local models are inserted into the global model right from the start of the computation. Hence, the local models cannot be adapted to progressing damage that might take place at different locations within the structure. In contrast to this, the more flexible coupled numerical homogenisation techniques often have the disadvantage of high computational efforts due to the complex connections between the separated models and the involved iteration procedures. The main disadvantages of many loose coupling approaches basically concern the limitations to non-automatic and one-way procedures, which do not allow an information exchange between the different modelling levels, and hence result in ignoring the potential mutual interaction of global and local effects. The available approaches in the field of two-way loose coupling also have some room for improvement. The multiscale projection method by Loehnert and Belytschko (2007) as well as the homogenisation-based iterative multiscale approach by Chrupalla et al. (2011) cannot be applied instantly to any standard FE software as both approaches require a certain modification of the source code. Within the adaptive progressive damage modelling technique by Labeas et al. (2012) the degraded material properties for the global model are determined rather simplified by taking the mean values of the local engineering constants, which might lead to wrong values in some cases. 1.5. Developed approach The objective of the present work regards the development of a novel FE-based two-way loose coupling approach that enables the simulation of composite structures incorporating the effect of geometrical as well as material nonlinearities, and also attains the interaction between both computation levels. In this way, the existing load-bearing capacities can be fully exploited and the structural weight and at the same time the related costs can be reduced. For appropriate applications it is desired that the two-way loose coupling approach requires less computational time compared to a full 3D simulation, which is the current alternative to investigate detailed progressive material damage. It has to be noted that in line with the previous description of the state-of-the-art, the presented simulation strategy can be seen as a new type of global/local methodology without scale separation. The proposed procedure aims at overcoming some of the aforementioned limitations of precedent approaches through an automatic two-way loose coupling with separated global and local models. This implicates a great flexibility regarding the creation of local models, as their positions and sizes may be adjusted during the coupling procedure according to the propagating damage estimations. In case of local damage, the substitutional material parameter for the global model are determined by means of virtual local part tests that deliver reliable results. Moreover, an effortless implementation of the coupling method is envisaged to enable a broad application within many common FE programs, which would be advantageous against some elaborate tight coupling implementations. The applicability of the proposed methodology is assessed by means of two exemplary applications: (1) a reference 3D model of a composite stiffened panel, and (2) a real composite stiffened panel (from Orifici et al. (2008)), for which the corresponding numerical and experimental data were reported. Both applications are subjected to axial compressive loading up to collapse and include geometrical as well as material nonlinearities. 2. Two-way loose coupling approach: General description This section is devoted to the description of the two-way loose coupling technique herein proposed. This numerical methodology is particularly designed for the investigation of material damage progression within composite structures prone to undergo geometrical instabilities (buckling) prior collapse, e.g. stiffened panels. It is worth mentioning that the present procedure has been integrated into the general purpose FE code Abaqus by means of several user-defined routines and Python scripts, emphasizing the high versatility of the envisaged strategy. The implementation of the method into other FE codes is likewise possible. The computational procedure essentially consists of the computation of three separated modules with two different levels of fidelity: 3

1. Global analysis - a simple and coarsely meshed global model of the complete structure that aims to capture the global buckling behaviour (global level, Section 2.1) 2. Local analysis - the generation and computation of high fidelity local models to investigate the intralaminar damage progression within certain parts of the structure (local level, Section 2.2) 3. Local part tests - a set of local characterisation tests to perform the computation of the degraded elastic properties which are sent back to the global model (Section 2.3) Figure 1 shows schematically the main ideas of the procedure along with the information exchanged between both levels for a generic case. On the top, the global FE model is displayed with a central critical area and, on the bottom, the finer meshed local FE model. In the following, the local area belonging to one global element is denominated as local part.

Figure 1: Overview of coupling between global and local model

The transition from the global to the local model is performed using the global nodal displacements, which are applied onto the corresponding local boundary nodes as prescribed kinematic conditions. For the opposite direction, in case of damage progression, a set of equivalent elastic properties that account for the stiffness degradation due to intralaminar failure are transferred from the local model to the respective global elements. This information exchange between the models is performed in an iterative way as reported in the flowchart of Figure 2. As the simulations at global and local level are not executed in parallel but one after the other, the nonlinear analyses are not paused in the process. More precisely, the equilibrium in the global and local models is always fulfilled before the respective information is transferred to the other level. The concept of the developed two-way loose coupling approach is presented in Figure 3 for the case of a stiffened composite panel, which represents the main motivation of this research. In this illustration, the three main modules of the procedure (global analysis, local analysis and local part tests) are clearly specified along with the information exchange between them and the employed computational capabilities. The consecutive coupling steps that comprise the numerical strategy are thoroughly described in the forthcoming sections. 2.1. Module 1: Global analysis The initial step of the coupling procedure is the creation of the global model. The panel geometry is represented with two-dimensional shell elements, whereas the composite material is described by means of the user-defined subroutine UMAT. This capability allows the modification of the initial elastic properties due to the potential intralaminar damage events that are evaluated at the local level. The external loading 4

Figure 2: Flowchart of the two-way loose coupling procedure

conditions at the global computation are imposed via prescribed displacements at the load introduction areas of the specimen. In general, the global model should be kept relatively simple but should preserve an adequate level of accuracy, so that quick and reliable results for the overall structural behaviour can be obtained from the global analysis. During the postprocessing stage corresponding to each load increment of the global simulation, computed lamina stress and strain components are compared within the chosen damage initiation criterion in order to assess if there exist any critical elements in the global model which can undergo intralaminar failure. These critical regions are usually denominated as “hot-spot areas”. In the current procedure, without loss of generality, the damage criterion proposed in Linde et al. (2004) is applied, that distinguishes between the two intralaminar failure modes fibre and matrix failure. Damage is predicted to occur when the index fm exceeds YT for the case of matrix failure (MF), and the index ff exceeds XT for the case of fibre failure (FF) as shown in Equation (1) and (2). These failure conditions read s    2 YT2 YT YT 2 fm = (22 ) + YT − 22 + (12 )2 > YT (1) YC YC SA s   X2 XT ff = (11 )2 + XT − T 11 > XT (2) XC XC In Equation (1), the strain components ij are referred to the local material setting at the lamina level. As usual the index 1 identifies the fibre direction, and 2 and 3 stand for transverse in-plane and out-of-plane to the fibre directions. The material strengths are the longitudinal tensile strength XT , the longitudinal compressive strength XC , the transverse in-plane tensile strength YT , the transverse in-plane compressive strength YC and the axial shear strength SA . It is worth mentioning that other failure criteria such as LaRC5

Figure 3: Two-way loose coupling procedure: application to a composite stiffened panel

type (D´ avila et al. (2005)) or Puck (Puck and Mannigel (2007)) could also be used in this computational strategy without any remarkable limitations. Once a certain area in the global model is predicted to fail, the location and size of the set of local models to be subsequently analysed can be accordingly determined. In particular, the determination of the size of the local models is defined by ensuring a regularised stress state at the area of interest. In this way, undesirable boundary effects due to imposition of the boundary conditions from the global level are prevented (Reinoso et al. (2012)). The adaptable character of the developed computational strategy with respect to the definition of the local models should be emphasized. This issue represents one of the main advantages with respect to other coupling procedures, due to the fact that in the present approach the identification of different local areas can be redefined along the loading process. However, it has to be noted also that this adaptability can limit the complete automation of the process, since the local model definitions are not fixed before starting the computations. 2.2. Module 2: Local analysis For the critical global regions, determined on the basis of the global simulation results, the local models are generated. Differing from the global shell model, the local models consist of three-dimensional solid elements with much denser meshes, and include a detailed material degradation model to investigate the intralaminar damage progression. Consistent with the global damage initiation criterion, the threedimensional version of the Linde degradation model (Linde et al. (2004)) is employed, which is available as a user subroutine in Abaqus Documentation (2013). The local damage analyses are performed through the application of the displacements from the global analysis to the boundary nodes of the local models lying on the delimiting surfaces by means of the shellto-solid submodelling functionality (see Abaqus Documentation (2013)). The different mesh densities at global and local level result in some nonconforming nodes in the global and local models. Depending on the position of the local nodes, the global displacements can be transferred directly to the local nodes (coincident local nodes) or they need to be interpolated (intermediate local nodes). 6

From the given strain state, in case the damage criterion is met, the matrix (dm ) and the fibre (df ) damage variables are computed as follows: dm = 1 −

YT [−C22 YT (fm −YT )Lc /Gm ] e fm

(3)

df = 1 −

XT [−C11 XT (ff −XT )Lc /Gf ] e ff

(4)

where Lc is the characteristic element length used to mitigate the mesh dependency relying on the crack band concept (Bazant and Oh (1983)), and Cij identifies the ij-component (i,j = 1,6) of the undamaged elasticity tensor. The effective elasticity tensor Cd defined in Linde et al. (2004) renders,  (1 − df )C11    Cd =    

(1 − df )(1 − dm )C12 (1 − dm )C22

(1 − df )C13 (1 − dm )C23 C33

symmetric

0 0 0 (1 − df )(1 − dm )C44

0 0 0 0 C55

 0 0   0  . 0   0  C66

(5)

Accordingly, the stress tensor σ and the algorithmic consistent tangent moduli C can be computed as follows: σ = Cd : ; (6) C=

X ∂di ∂σ = Cd + ⊗ σ0 , ∂ ∂

(7)

i=d,m

where, as customary, σ0 = C0 :  is the effective stress tensor and C0 stands for the intact elasticity tensor that contains the original mechanical properties (Simo and Ju (1987)). Based on this damage model, the values of the predicted degraded material properties are stored in order to be employed in the subsequent local part tests (see Section 2.3). 2.3. Module 3: Local part tests To return the information concerning the local damage extent estimation from the local level to the global level, an additional “homogenisation-like” step is necessary as the local mesh is much finer as compared to the global mesh and thus, the updated material properties cannot directly be given to the global model. As a consequence of the plane stress state that is defined at the global level, three linear characterisation tests are performed through the use of the so-called local parts to determine equivalent engineering constants for the global shell model. Figure 4 shows the main principle of the two tensile tests (leftmost and central sketch) and the shear test (rightmost sketch). It has to be noted that the local parts are actually three-dimensional models, but for the sake of convenience simplified two-dimensional drawings are displayed. For the sake of clarity, 3D FE models with the corresponding displacements of the three tests are depicted in Figure 5. Basically, one local part test comprises the local elements that are associated with one particular laminate layer of one global element. In addition, each of these local part tests is characterised by a linear material definition and contains the final degraded material properties from the local analysis. Consequently, for any damaged global element and laminate layer, a set of updated global engineering constants can be computed and returned by means of the associated local part. Due to transversely isotropic UD layers definition at global level, the four independent material parameter E11 , E22 , G12 and ν12 can be determined by means of the following tests: • Test 1: Tension in global 1-direction for the calculation of E11 = σ11 /11 and ν12 = −22 /11 7

• Test 2: Tension in global 2-direction for the calculation of E22 = σ22 /22 • Test 3: Shear in global 1-2 plane for the calculation of G12 = τ12 /γ12

Figure 4: Local part tests (Left: Test 1, Middle: Test 2, Right: Test 3)

Figure 5: Displacements of local part tests (Left: u1 of Test 1, Middle: u2 of Test 2, Right: u2 of Test 3)

2.4. Transformation process Since the global model requires the material properties in local directions, i.e. parallel and perpendicular to the fibre orientation, a transformation from the global to the local directions is needed to be accomplished for fibre orientations different from 0◦ . This transformation process is performed through a set of equations, that are usually employed for the determination of the global stiffness matrix from the known engineering constants of the single composite plies. Here, the customary process is performed in the reverse way, starting from the partly known global stiffness matrix Q and ending up with the local stiffness matrix Q0 and the respective local engineering constants. For the global stiffness matrix Q the values for Q11 , Q22 , Q33 , Q12 and Q21 can be determined from the results of the local part tests, i.e. the engineering constants related to the global setting, whereas the entries Q13 , Q23 , Q31 and Q32 are initially unknown, see Equation (8). The entries of the local stiffness matrix Q0 are completely unknown, with the exception of those equal to zero that result from the transversely isotropic definition of the UD layer, see Equation (9). These two matrices (Q and Q0 ) are related through the standard transformation matrix T , that depends on the fibre orientation and whose definition is given in Equation (10). The determination of the previous set of nine unknowns is carried out by solving the system of equations expressed in Equation (11), that is automatized using a Matlab script. Finally, the desired engineering constants, that are used as entry data in the global UMAT, can be computed via the local flexibility matrix, see Equation (12).   Q11 Q12 Q13 Q = Q21 Q22 Q23  (8) Q31 Q32 Q33 8

 0 Q11 Q0 = Q021 0 cos2 (α) 2  sin (α) T = −2sin2 (α) 

Q012 Q022 0

 0 0  Q033

(9)

 sin2 (α) sin(α)cos(α) cos2 (α) −sin(α)cos(α)  2sin(α)cos(α) cos2 (α) − sin2 (α)

(10)

Q = T T · Q0 · T

(11)

S = Q0

(−1)

(12)

2.5. Iteration process The new material properties, that are obtained from the local part tests, are then incorporated into the associated global elements using the global user subroutine UMAT that reads them from an external auxiliary file. In this regard, in line with Labeas et al. (2012), the current solution of the global model does not represent an equilibrium state any more due to the modification of the elastic properties accomplished at the local level, and furthermore, an update of the global displacements is required. This implies that the global simulation should be recomputed within a global/local iteration loop from which the updated global displacements are newly obtained. Subsequently, the new global displacements are imposed as updated boundary conditions to the local models that were defined in the precedent iteration. However, such recomputed global displacements might also provoke the development of further intralaminar damage at the local level, so that the degraded global material properties need to be recomputed as well. The converged state can be obtained after carrying out this iteration loop until the material damage state does not change within certain limits between two consecutive iterations, compare Figure 2. Once, this convergence condition is satisfied, the iteration process can be regarded as finished and the final global state is reached. Based on experience, usually two or three iteration runs are sufficient to reach a stable damage state. Subsequently, the next coupling loop can be started, i.e. the loading of the global model can be increased until the global damage initiation criterion indicates the development of further damage within the specimen. Figure 6 shows an exemplary load-displacement curve for the global model of the two-way loose coupling procedure. It can be seen that for the initial global analysis the original material properties are used up to the global displacement load u1 resulting in the initial reaction force F1,initial (black curve). In case that the global criterion indicates damage initiation in the composite material (red point) a progressive damage simulation is performed for the critical area by means of a local model. Thereby the used displacement, here denoted as u1 , is slightly higher than the critical displacement u1,crit , to capture not only the damage initiation in the local model, but already a progressed material degradation. After the local degradation analysis and the virtual local part tests have been completed the updated global material properties are used for the next global iteration loop run. In doing so, the global simulation is repeated, at first with the initial undamaged material properties and then, from the critical loading u1,crit on, using the new degraded values (red curve). This process is iterated until the first coupling loop run is converged to the final reaction force F1 (black square). In the diagram only the final curves are depicted, the different iteration runs are not shown. The next coupling loop run is conducted accordingly. The load is raised to the displacement u2 that is shortly higher than the second critical loading u2,crit , at which the global criterion shows a damage propagation either because the existing critical global region enlarges or another critical area arises (blue point). For the following global iterations the original material properties (black curve), the degraded and converged properties from the first coupling loop run (red curve) and finally the current properties from the second coupling loop run (blue curve) are successively applied. An overall equation for the two-way loose coupling process can be formulated as follows: global F k+1,i+1 = K global · (uk+1,crit − uk,crit ) + K global original · u1,crit + ... + K k k+1,i+1 · (uk+1 − uk+1,crit )

9

(13)

Figure 6: Exemplary load-displacement curve of the two-way loose coupling procedure

3. Application This section describes the application of the two-way loose coupling procedure to two examples: • The first example examines the applicability of the proposed methodology through its application to an academic composite stiffened panel. The results obtained using the coupling procedure are compared with those corresponding to a reference solid model. The reference model and the local models that are used for the coupling procedure have the same mesh density, along with the same material definition to inspect intralaminar failure based on the continuous damage model proposed by Linde et al. (2004). • The second example concerns the application to a real composite stiffened panel. The results of the coupling computations are correlated with the experimental and numerical data reported in Orifici et al. (2008). 3.1. Panel 1: Academic application The first application under analysis is a flat stiffened panel with one centered T-stringer, see Figure 7. The system comprises the following entities: (1) the primary skin, (2) the stringer and (3) the adhesive layer of 0.2 mm in thickness that is located between the previous parts. The geometrical characteristics are listed in Table 1. The laminate of the skin and the stiffener has a symmetric composite layup with [0◦ , 90◦ ]S , where the unidirectional layers have a nominal thickness of 0.25 mm. The properties of the composite material and the adhesive layer are reported in Table 2, where the longitudinal direction of the panel equals the reference material orientation. The boundary conditions applied to the panel are the following: (1) one transverse edge of the panel is fully clamped, (2) the opposite transverse edge is partly fixed with a free deformation along the longitudinal direction (imposed displacement), and (3) free conditions at both longitudinal edges. Simulations are conducted until buckling and damage are estimated to occur. 3.1.1. Global model: Geometrically nonlinear postbuckling simulation The global model of the panel is created with standard linear shell elements (S4R of Abaqus), using an approximate element-side length of around 5 mm (resulting in 280 elements in total). To facilitate the transition between the prebuckling and the postbuckling evolution ranges, the model is geometrically perturbed using the kinematic field corresponding to the first eigenmode of a precedent linear buckling simulation. The composite damage criterion by Linde et al. (2004) is used as a postprocessing utility to determine the critical areas of the structure, which correspond to the locations for the detailed local models. 10

Description Panel length Panel width Stringer width Stringer height Laminate thickness Adhesive thickness

Symbol l w wstringer h t tglue

Value 100 mm 40 mm 20 mm 8 mm 1 mm 0.2 mm

Table 1: Geometry of stiffened panel

Figure 7: Section geometry of stiffened panel [mm]

3.1.2. Local models: Material nonlinear simulations The local models corresponding to each of the critical areas consist of three-dimensional hexahedral elements (C3D8 of Abaqus), using an approximated in-plane and out-of-plane side-lengths equal to 1 mm and 0.25 mm, respectively, so that in thickness direction one element represents one laminate layer. The size of the local models is chosen according to the size of the critical areas that are identified in the global simulation. The local models are created automatically using a Matlab preprocessor that requires the desired size of the local model and its position within the global model as manual input. At this level of analysis, to preserve the consistency, the material nonlinearity is incorporated by means of the user subroutine UMAT with the degradation model proposed in Linde et al. (2004). 3.1.3. Reference model: Geometrically and material nonlinear simulation The reference model is a full 3D model with 32800 elements (C3D8 of Abaqus) with the same mesh density as the local models. This model also comprises the same material model by Linde et al. (2004) to investigate intralaminar damage events. Similar to the global shell model, the first eigenmode of the stiffened panel is applied on the reference model as an initial geometrical imperfection. 3.1.4. Coupling results The coupling procedure is performed with four coupling loops, each of them with several iteration loops. For each coupling loop run, the displacement loading is increased in such a way that the critical global areas expand, and accordingly the size of local models are enlarged. Figure 8 shows overlay plots of the separated global shell and the local volume models for the four successive coupling loops. Along the global simulation, the damage initiation criterion detects the first potential material damage at a compressive displacement equal to 0.56 mm (Figure 9.a). At this load level, damage estimations correspond to matrix cracking at the longitudinal borders of the stringer flange. For these regions two local models are created, and the first coupling loop is repeated until convergence is reached. Both local damage analyses confirm the occurrence of matrix failure (Figure 9.b) and the related global material properties are updated according to the local degradation procedure. The second coupling loop conducts a displacement of 0.60 mm and implies a progression of the matrix damage at both locations previously identified. After the second loop has been completed, the global loading is further increased to a displacement of 0.63 mm until an additional critical area related to fibre failure is predicted at the middle of the stiffener, see Figure 10.a and 10.b. It should be mentioned that up to this 11

Stiffness properties Young’s Modulus in 1-direction E11 Young’s Modulus in 2-direction E22 Shear Modulus in 12-plane G12 Poisson’s ratio ν12 Young’s Modulus of adhesive Eglue Poisson’s ratio of adhesive νglue -

Value 146.5 GPa 9.7 GPa 5.1 GPa 0.28 3.0 GPa 0.4 -

Strength and fracture properties Tensile strength in 1-direction XT Compressive strength in 1-direction XC Tensile strength in 2-direction YT Compressive strength in 2-direction YC Shear strength in 12-plane SA Fracture energy of fibre Gf Fracture energy of matrix Gm

Value 2.583 GPa 1.483 GPa 0.092 GPa 0.270 GPa 0.106 GPa 12.5 N/mm 1.0 N/mm

Table 2: Material data for composite and adhesive

Loop 1: u = 0.56 mm

Loop 2: u = 0.60 mm

Loop 3: u = 0.63 mm

Loop 4: u = 0.67 mm

Figure 8: Overlay plots of global and local models from coupling loops 1 to 4

Figure 9: Coupling loop 1. (a) First critical global areas related to matrix failure, (b) Matrix damage initiation for the Regions 1 and 2

third coupling loop the local material damage is relatively small, and therefore these damage events do not influence the global structural stiffness distinctly. Subsequently, in the fourth coupling loop the panel is further shortened up to 0.67 mm. This load increase leads to a further propagation of all three damage locations, so that a combination to one large local model for the complete critical area is suitable. Along the process, the local model is damaged considerably, which 12

results in a notable stiffness reduction of the global structure as well. Consequently, the total collapse of the panel the coupling procedure can be regarded as finished. The final matrix and fibre damage of the combined local model is shown in Figure 11.a and Figure 11.b, respectively.

Figure 10: Coupling loop 3. (a) First critical global area related to fibre failure, (b) Fibre damage initiation of middle local model

Figure 11: Coupling loop 4. (a) Final matrix damage of the local model, (b) Final fibre damage of the local model

The results of the four consecutive coupling loops are shown as a combined load-displacement curve in Figure 12. The coupling curve is in good agreement with the full 3D reference curve that is displayed as well. The structural stiffness agrees in both analyses at the beginning of the loading as well as after buckling, which can be seen from the coincident slopes of the curves. Buckling occurs at a displacement of around 0.14 mm and 0.15 mm in the coupling and in the reference simulation, respectively, and leads to a small stiffness reduction of the structure. This is visible as a slight kink in the load-displacement curve. Some important results of the coupling and the reference simulation are summarised in Table 3 by means of a list with the respective loading displacements of buckling, matrix and fibre damage initiation along with the total collapse load. The initiation of material damage is found at similar displacements and locations in the two analyses. The onset of degradation happens at a longitudinal compression of 0.54 mm in the coupling analysis, whereas it is predicted to occur at 0.55 mm in the reference solid analysis. At this point, matrix cracking starts in two areas in the stringer flange close to the stringer web in the upper laminate layer with a fibre orientation of 0◦ . As was previously stated, this first interfibre failure has no noticeable influence on the total stiffness of the panel. The matrix damage increases mainly in the stringer flange area until fibre failure initiates at a displacement of 0.62 mm and 0.63 mm for the reference and coupling computations, respectively. This latter damage event is predicted to occur at the stringer flange edge in the middle of the panel again in the upper 0◦ -ply as can be seen in Figure 11.b. Fibre damage propagates rapidly in a path throughout the panel width and finally results in the total damage of the panel at a shortening of 0.63 mm, that corresponds to a maximum load of 22 kN in the coupling analysis and 22.5 kN in the reference analysis. 13

Figure 12: Load-displacement curves of the coupling and reference analysis

In both analyses, a strong reduction of the panel stiffness leading to a distinct drop in the load-displacement curves constitutes the final failure as depicted in Figure 12. Based on the previous results, it can be seen that the coupling procedure is capable to approximate the results of the full 3D reference model with a high level of accuracy. Simulation Coupling Reference

Buckling 0.14 mm 0.15 mm

Matrix damage 0.54 mm 0.55 mm

Fibre damage 0.63 mm 0.62 mm

Total failure 0.63 mm 0.63 mm

Table 3: Comparison between coupling and reference results: Loading displacements of initiation points

3.2. Panel 2: Real composite stiffened panel In this section, the application of the two-way loose-coupling procedure to a real stiffened composite panels taken from the literature is examined. The chosen specimen have been investigated within the COCOMAT project in Orifici et al. (2008), which corresponds to the typology denoted as D1 in that investigation. Figure 13 shows a sketch of the panel configuration, whereas the geometry data and the material disposal for each region are reported in Tables 4 and 5, respectively. The skin and the stringer of the panel are manufactured by the unidirectional CFRC material IM7/8552 UD, whose mechanical properties are reported in Table 6. The reference material orientation accords with the longitudinal axis of the panel. Description Total length Free length Panel width Stringer width Stringer height

Symbol L Lf b w h

Value 400 mm 300 mm 64 mm 32 mm 14 mm

Table 4: Geometry data of the panel D1 (from Orifici et al. (2008))

14

Figure 13: Geometry of stiffened panel designs D1 (from Orifici et al. (2008))

Region Skin layup Stringer flange layup Stringer web layup

Laminate [90, ±45, 0]S [06 , (±45)3 ]S [(±45)3 , 06 ]S

Symbol tskin tflange tweb

Thickness 1 mm 1.5 mm 1.5 mm

Table 5: Laminate definition of the panel D1 (from Orifici et al. (2008))

3.2.1. Reference models To evaluate the results of the coupling procedure two different reference models are used: (i) one solid model, and (ii) one shell model. Both models are equipped with the same material degradation model (Linde et al. (2004)). After a preliminary mesh convergence analysis with respect to the lowest buckling load of the panel, 2700 shell elements (S4R of Abaqus), with a side-length equal to 4 mm, are selected for the reference shell model, whereas 680800 volume elements (C3D8 of Abaqus), with an in-plane side-length equal to 1 mm and one element per layer over the thickness, are employed for the reference solid model. For each of these reference models, the first eigenmode is used as an initial geometric imperfection to perturb the nominal geometry of the panel in the corresponding postbuckling simulations. A description of the applied boundary conditions to both models is shown in Figure 13, which are: (1) one transverse edge fully restrained (clamped end in Figure 13), (2) at the opposite transverse edge, all the kinematic degrees of freedom are clamped except the longitudinal one (loaded end in Figure 13), and (3) both longitudinal edges are completely free. 3.2.2. Global model: Geometrically nonlinear postbuckling simulation The global coupling models have in principle the same structure as the reference shell model. The only distinction with respect to this reference model regards the material definition, that, in the global model, is defined through a user-defined subroutine with a linear elastic material definition (global shell UMAT). 3.2.3. Local models: Material nonlinear simulations The local coupling models are created with fully integrated 8-node volume elements (C3D8 of Abaqus). The selected mesh densities are in accordance with the FE meshes of the reference solid model. The total size of the local models is dependent on the critical areas, that are determined within the global analysis. 3.2.4. Coupling results As a consequence of the particular mechanical performance of this panel, the coupling procedure is performed with one single coupling loop. According to the used damage criterion (Linde et al. (2004)), the 15

Stiffness properties Young’s Modulus in 1-direction E11 Young’s Modulus in 2-direction E22 Shear Modulus in 12-plane G12 Shear Modulus in 13-plane G13 Shear Modulus in 23-plane G23 Poisson’s ratio ν12 -

Value 147 GPa 11.8 GPa 6 GPa 6 GPa 4 GPa 0.28 -

Strength and fracture properties Tensile strength in 1-direction XT Compressive strength in 1-direction XC Tensile strength in 2-direction YT Compressive strength in 2-direction YC Shear strength in 12-plane SA Fracture energy of fibre Gf Fracture energy of matrix Gm

Value 2.583 GPa 1.483 GPa 0.092 GPa 0.270 GPa 0.106 GPa 12.5 N/mm 1.0 N/mm

Table 6: Material data of the unidirectional CFRC material IM7/8552 UD

first two critical areas are related to fibre damage starting at an axial compressive displacement of around 1.89 mm at the upper region of the stringer blade. Respectively, two local models are created for each of these regions, see Figure 14. Along the coupling procedure, the evolution of the fibre and matrix damage events in those two local models during the four iteration loops can be seen in Figures 15–18. In these plots, notable fibre and matrix damage events are developed in both local models. Moreover, it is visible that there is no substantial change between the depictions corresponding to the iteration loops 3 and 4. Hence, already after the third iteration loop, the local damage state can be regarded as constant and the iteration procedure as concluded.

Figure 14: Critical global areas (left) and overlay plot of respective local models at global model location (right) of panel D1

Figure 15: Fibre damage of local model 1 of panel D1 (Iteration loops 1 to 4)

Figure 16: Matrix damage of local model 1 of panel D1 (Iteration loops 1 to 4)

16

Figure 17: Fibre damage of local model 2 of panel D1 (Iteration loops 1 to 4)

Figure 18: Matrix damage of local model 2 of panel D1 (Iteration loops 1 to 4)

Figure 19.a displays the load-displacement curves of the four iteration loops obtained from the coupling procedure. The four coupling loops show different drop-offs of the evolution curve after damage initiation, marking the total failure of the panel. By reaching this initiation point, the coupling procedure can be considered as finished, so that no further coupling loop including a load increase is necessary. In Figure 19.b, the resulting load-displacement curve of the coupling method is compared with the reference shell and solid models and with the results reported in Orifici et al. (2008). The coupling results satisfactory agree with the shell reference results, where the difference between the maximum load levels is less than 1%. With respect to the experimental and numerical results given in Orifici et al. (2008), it can be observed that the coupling results agree very well with the experimental data, whereas the other numerical simulations slightly overestimate the maximum load achieved during the test. Particularly, the maximum load from the Orifici simulation deviates about 15% from the maximum load of the coupling results. This difference might be presumably caused by the use of different damage criteria and degradation models in those computations along with the role of interlaminar damage events. In comparison with the solid reference analysis, the coupling results are mildly deviated with a difference of around 10% between the predicted maximum load levels. It is worth mentioning that this difference is also observable between the reference shell and solid models. Potential causes of this difference might be attributed to the role of the interlaminar stress components (τ13 , τ23 and σ33 ) that are taken into account in solid models but neglected in shell simulations. To support this argument, Figure 20 depicts the respective normal and shear stresses of the solid reference model at the onset of degradation. It is visible that the stresses related to the thickness direction have the same magnitude as some other stresses, so that they should not be neglected and certainly have an influence on the damage behaviour. This explains that in the present case the solid reference model reports material damage at a different loading as the shell reference model. Nevertheless, the load-displacement evolution curves corresponding to the reference shell and the coupling procedure are below the reference solid curve, and therefore they do not overestimate the ultimate load and remain on the safe side. A more comprehensive analysis of this feature is beyond the scope of this work. A summary of the deviations between the maximum load of the coupling simulation with respect to the numerical and the experimental data reported in Orifici et al. (2008) is given in Table 7. Finally, regarding the type and location of the material damages, it is worth mentioning that the present simulations are in close agreement with the experimental and numerical results from the literature. 4. Conclusions In this paper, a novel two-way loose coupling approach for an efficient and reliable simulation of composite structures, which allows the combined incorporation of material and geometrical nonlinearities, has been 17

Figure 19: (a) Load-displacement curves of coupling iterations for panel D1, (b) Load-displacement curves of coupling analysis for panel D1 in comparison with reference results and results by Orifici et al. (2008)

Figure 20: Details of stresses at damage initiation in solid reference model for panel D1

developed. This strategy is based on the analysis of separated FE models with different levels of fidelity: (i) a coarse mesh global model, and (ii) a set of detailed local models which include numerical techniques to represent intralaminar damage events. The two-way loose coupling character of this approach permits the interrelation between the global and local levels throughout the numerical analysis, so that the size and locations of the local models can be suitably adapted. Single composite stiffened panels have been considered as benchmark applications of the developed technique. In particular, the suitability of the method has been examined by means of two composite applications: (1) an academic composite stiffened panel, and (2) a real T-stringer specimen previously investigated in Orifici et al. (2008). Regarding the first panel, results from the coupling computations have been contrasted with the simulations of a full 3D reference model. The load-displacement curves of the solid reference model and the coupling simulation show a satisfactory agreement as well as coinciding predictions with regard to the location and the sequence of the intralaminar material damages. With respect to the second panel, the developed numerical strategy accords to the reference results in terms of the peak load level achieved during the reference test and the location and sequence of intralaminar damage events. This agreement of the coupling and the respective reference results has demonstrated the effectiveness 18

Solid reference Shell reference Orifici simulation Orifici experiments

Other results max. Load [kN] 52.05 47.46 55.00 50.00

Coupling result max. Load [kN] 47.52 47.52 47.52 47.52

Rel. deviation [%] 9.53 0.13 15.73 5.21

Table 7: Relative deviations between maximum reaction forces of coupling result and other results for panel D1

and applicability of the proposed novel approach in simulating postbuckling behaviour of composite stiffened panels up to collapse. In future developments this methodology might be enhanced by introducing the capability to simulate the initiation and evolution of interlaminar failures. 5. Acknowledgements Parts of the research leading to these results has received funding from the European Community’s Seventh Framework Programme (FP7/2007-2013) under grant agreement number 234147. JR acknowledges the support of the Spanish Ministry of Economy and Competitiveness (DPI2012-37187) and the Andalusian Government (Project of Excellence No. TEP-7093). References References F. Feyel, A multilevel finite element method (FE2 ) to describe the response of highly non-linear structures using generalized continua, Computer Methods in Applied Mechanics and Engineering 192 (28–30) (2003) 3233–3244. V. Kouznetsova, M. Geers, W. Brekelmans, Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy, Computer Methods in Applied Mechanics and Engineering 193 (48–51) (2004) 5525–5550. K. Terada, N. Kikuchi, A class of general algorithms for multi-scale analyses of heterogeneous media, Computer Methods in Applied Mechanics and Engineering 190 (40–41) (2001) 5427–5464, ISSN 0045-7825. J. Fish, Q. Yu, K. Shek, Computational damage mechanics for composite materials based on mathematical homogenization, International Journal for Numerical Methods in Engineering 45 (11) (1999) 1657–1679. S. Ghosh, K. Lee, P. Raghavan, A multi-level computational model for multi-scale damage analysis in composite and porous materials, International Journal of Solids and Structures 38 (14) (2001) 2335–2385. C. Miehe, A. Koch, Computational micro-to-macro transitions of discretized microstructures undergoing small strains, Archive of Applied Mechanics 72 (4-5) (2002) 300–317. T. J. Hughes, G. R. Feij´ oo, L. Mazzei, J.-B. Quincy, The variational multiscale method - a paradigm for computational mechanics, Computer Methods in Applied Mechanics and Engineering 166 (1–2) (1998) 3–24. H. Kadowaki, W. K. Liu, Bridging multi-scale method for localization problems, Computer Methods in Applied Mechanics and Engineering 193 (30–32) (2004) 3267–3302. I. M. Gitman, Representative volumes and multi-scale modelling of quasi-brittle materials, Ph.D. thesis, TU Delft, Delft University of Technology, 2006. C. D. Mote, Global-local finite element, International Journal for Numerical Methods in Engineering 3 (4) (1971) 565–574. A. K. Noor, Global-local methodologies and their application to nonlinear analysis, Finite Elements in Analysis and Design 2 (4) (1986) 333–346. K. M. Mao, C. T. Sun, A refined global-local finite element analysis method, International Journal for Numerical Methods in Engineering 32 (1) (1991) 29–43. E. Rank, A Zooming – Technique Using a Hierarchical hp–Version of the Finite Element Method, in: The Mathematics of Finite Elements and Applications - Highlights, John Wiley & Sons, Uxbridge, 1–9, 1993. J. Reinoso, A. Bl´ aquez, A. Estefani, J. Ca˜ nas, E. Ar´ evalo, F. Cruz, Experimental and three-dimensional global-local finite element analysis of a composite component including degradation process at the interfaces, Composites: Part B (43) (2012) 1929–1942. T. Zohdi, M. Feucht, D. Gross, P. Wriggers, A description of macroscopic damage through microstructural relaxation, International Journal for Numerical Methods in Engineering 43 (3) (1998) 493–506. G. Ernst, M. Vogler, C. H¨ uhne, R. Rolfes, Multiscale progressive failure analysis of textile composites, Composites Science and Technology 70 (1) (2010) 61–72. S. Loehnert, T. Belytschko, A multiscale projection method for macro/microcrack simulations, International Journal for Numerical Methods in Engineering 71 (12) (2007) 1466–1482.

19

L. Gendre, O. Allix, P. Gosselet, F. Comte, Non-intrusive and exact global/local techniques for structural problems with local plasticity, Computational Mechanics 44 (2) (2009) 233–245. G. Labeas, S. Belesis, I. Diamantakos, K. Tserpes, Adaptative progressive damage modeling for large-scale composite structures, International Journal of Damage Mechanics 21 (3) (2012) 441–462. D. Chrupalla, S. Berg, L. K¨ arger, M. Doreille, T. Ludwig, E. L. Jansen, R. Rolfes, A. Kling, A Homogenization-Based TwoWay Multiscale Approach for Composite Structures, in: R. Rolfes, E. L. Jansen (Eds.), Proceedings of the 3rd ECCOMAS Thematic Conference on the Mechanical Response of Composites, Hannover, Germany, 263–270, 2011. A. Orifici, I. Ortiz de Zarate Alberdi, R. Thomson, J. Bayandor, Compression and post-buckling damage growth and collapse analysis of flat composite stiffened panels, Composites Science and Technology (68) (2008) 3150–3160. P. Linde, J. Pleitner, H. de Boer, C. Carmone, Modelling and simulation of fibre metal laminates, in: ABAQUS Users’ conference, 421–439, 2004. C. G. D´ avila, P. P. Camanho, C. A. Rose, Failure criteria for FRP laminates, Journal of Composites Materials (39) (2005) 323–345. A. Puck, M. Mannigel, Physically based non-linear stressstrain relations for the inter-fibre fracture analysis of FRP laminates, Composites Science and Technology (67) (2007) 1955–1964. Abaqus Documentation, Abaqus 6.13 Documentation, Dassault systemes, 6.13 edn., 2013. Z. Bazant, B. Oh, Crack band theory for fracture of concrete, Materials and structures (16) (1983) 155–177. J. C. Simo, J. W. Ju, Strain and stress-based continuum damage models-I. Formulation, International Journal of Solids and Structures 23 (1987) 981–840.

20