Nonlinear Homogenization in Masonry Structures

11th World Congress on Computational Mechanics (WCCM XI) 5th European Conference on Computational Mechanics (ECCM V) 6th European Conference on Computational Fluid Dynamics (ECFD VI) E. Oñate, J. Oliver and A. Huerta (Eds) NONLINEAR HOMOGENIZATION IN MASONRY STRUCTURES Georgios A. Drosopoulos∗,† , Maria E. Stavroulaki† , Konstantinos Giannis† , Leonidas Plymakis† , Georgios E. Stavroulakis† and Peter Wriggers∗ ∗ Leibniz University of Hannover, D-30167 Hannover, Germany e-mail: drosopoulos@ikm.uni-hannover.de, wriggers@ikm.uni-hannover.de web page: http://www.ikm.uni-hannover.de † Technical University of Crete, GR-73100 Chania, Greece e-mail: mstavr@mred.tuc.gr; gestavr@dpem.tuc.gr - web page: http://www.comeco.tuc.gr Key words: Multi-scale, Homogenization, Interpolation, Masonry, Contact, FE2 Abstract. Numerical homogenization is based on the usage of finite element analysis for the description of average properties of materials with heterogeneous microstructure. The practical steps of a computational homogenization approach and representative examples related to masonry structures and ceramic materials are presented in this article. The non-linear Representative Volume Elements (RVEs) of a masonry structure, including parts with elastoplastic material behaviour (mortar) and a ceramic material with a unilateral contact interface (crack), are created and solved. Parametric analysis has been chosen and used for the description of the strain loading. Results concerning the average stress and strain in the RVE domain are then calculated. In addition, the stiffness is estimated for each loading level. Finally, two databases for the stiffness and the stress-strain data are created, a metamodel based on MATLAB interpolation is used, and an overall non-linear homogenization procedure (FE2 ), is considered. The good comparison with direct heterogeneous macroscopic models created by commercial software shows that the proposed method can be used for the simulation of non-linear heterogeneous structures. 1 INTRODUCTION In this article, a multi-scale, computational homogenization method is presented, for the study of non-linear, heterogeneous materials. Several different approaches have been proposed in the past for the investigation of these materials. Analytical/mathematical methods, like asymptotic homogenization [1], can be more accurate in the description of the micro structure, for relatively simple microscopic patterns and constitutive laws. On the other hand, numerical methods may be used for the simulation of complex mi1 G.A. Drosopoulos, M.E. Stavroulaki, K. Giannis, L. Plymakis, G.E. Stavroulakis and P. Wriggers croscopic geometries, over a statistically defined representative amount of material [2]. These methods are indispensable for non-linear problems. Numerical/computational homogenization can be extended to cover several non-linear effects, like contact, debonding, damage and plasticity [3]. According to numerical homogenization, a unit cell is explicitly solved and the resulting average quantities are then used for the determination of the parameters of a macroscopic constitutive law [4]. For non-linear problems numerical homogenization is performed at each load level. From another point of view, multi-level computational homogenization incorporates a concurrent analysis of both the macro and the microstructure, in a nested multi-scale approach [5, 6, 7, 8, 9, 10]. Within this method, the macroscopic constitutive behaviour is determined during simulation, after solving the microscopic problem and transferring the necessary information on the macroscopic scale. This approach, which is generally called FE2 , offers the flexibility of simulating complex microstructural patterns, with every kind of non-linearity. In the present work, parametric analysis is used for the simulation of a non-linear Representative Volume Element (RVE), under several loading paths and loading levels. After solution of the microscopic structure, the average stress is estimated and a strain-stress database is created. In addition, stiffness information is obtained for each particular loading level and a second strain-stiffness database is obtained. Based on these databases, and on a MATLAB-based interpolation for the creation of a metamodel, an overall computational homogenization model, in a FE2 sense, is created for the simulation of a macroscopic structure. Comparison with direct heterogeneous macroscopic models shows that the adopted procedure leads to satisfactory results. The method can be applied to different RVEs, with different non-linear microscopic behaviour. In the present work two different microscopic models have been chosen: a non-linear masonry RVE and a ceramic material with a unilateral contact interface as a potential crack. 2 INTRODUCTION TO COMPUTATIONAL HOMOGENIZATION The approach adopted in this article is related to the concurrent analysis of the macroscopic and the microscopic structure, respectively. According to the classical formulation of the multi-scale computational homogenization [5, 6], two nested boundary value problems are concurrently solved. The initial heterogeneous macroscopic structure is equivalent with a homogeneous one, in each Gauss point of which a suitably defined RVE is correlated. This RVE includes every heterogeneity and non-linearity of the material. According to the Hill-Mandel condition or energy averaging theorem, the macroscopic volume average of the variation of work equals to the local work variation, on the RVE  [11]: 1 M M σ :ǫ = σ m : ǫm dVm (1) Vm Vm Three type of loading states, which satisfy the above condition, can be applied to the RVE: a) prescribed linear displacements, b) prescribed tractions, c) periodic boundary 2 G.A. Drosopoulos, M.E. Stavroulaki, K. Giannis, L. Plymakis, G.E. Stavroulakis and P. Wriggers conditions. With linear or periodic boundary conditions, a macroscopic strain is the loading of the RVE. After analysis and convergence of each RVE in every Gauss point, results concerning the average stress and the stiffness are given back to the macroscopic structure, Fig. 1. No assumption for the constitutive law of the macroscopic structure is a priori considered, thus the macroscopic constitutive behaviour is numerically obtained. This is a practical solution to the major question of homogenization, namely which are the properties of the homogeneous constitutive law. In this work the microscopic calculations of the RVE, which normally take place in each Gauss point and time step of the macroscopic model, are substituted with two databases, which carry the same information: average stress and stiffness of the macroscopic model. To obtain this information, several numerical simulations on the masonry RVE are explicitly considered. In the next paragraphs the steps of the proposed procedure are given in details. Figure 1: Schematic representation of the multi-scale, concurrent, computational homogenization 3 GENERAL FRAMEWORK OF THE PROPOSED APPROACH The key idea of the present work is to replace the microscopic simulation of the RVE, which is required within each iteration step of the computational homogenization method, with two databases containing information related to the stress and the stiffness of the macro model. This information is transferred back to the macroscopic structure. Thus, instead of solving the RVE in each Gauss point and iteration loading or time step, which is a time consuming procedure, a metamodel is used, or in other words an interpolation of the proper quantity from the databases is considered. This concept has the following steps: a) Creation of the RVE, for example within COMSOL Multiphysics. Two microscopic structures have been created, one representing a masonry RVE and another representing a ceramic material with a potential unilateral crack, as is explained in details in the following section. b) Linear displacement boundary conditions are applied to the RVE; a number of different loading strain combinations have been considered, by incorporating three parameters to the equations describing the linear displacement conditions on the boundaries of the 3 G.A. Drosopoulos, M.E. Stavroulaki, K. Giannis, L. Plymakis, G.E. Stavroulakis and P. Wriggers RVE. By using these parameters, the space of the strain vector is scanned, resulting in a significant number of possible strain loads. c) After analysis of each RVE, the average stress and strain are calculated. d) Steps b) and c) are repeated, but now with three test incremental loading strain vectors, which are applied to the boundaries of the RVE. Then, by incrementally solving the Hooke’s law, stiffness information is obtained for the particular loading conditions. e) After the previous steps, two databases have been created: one that corresponds strains to stresses and another that corresponds strains to stiffness information. These are incorporated in an overall FE2 computational homogenization scheme developed with MATLAB, for the simulation of the respective macroscopic structures. f) Comparison of the results with direct heterogeneous macroscopic models created in other commercial software packages is used to evaluate the whole procedure. 4 4.1 THE MICROSCOPIC MODELS The masonry RVE The masonry RVE consists of the brick parts and the mortar joints, thus the material that connects the bricks. The non-linearity of this model is concentrated on the mortar joints by using a perfect plasticity law, while the brick parts are considered to be linear. The dimensions and the mesh of the RVE are shown in Fig. 2. The mesh consists of rectangular plane stress elements, with out of plane thickness equal to 70mm. In addition, material properties have been taken from the literature: Eb = 4865N/mm2 , nb = 0.09 for the brick parts and Em = 1180N/mm2 , nm = 0.06 for the mortar joints. A tensile strength of 0.9N/mm2 has been also used for the mortar joints. (a) (b) Figure 2: (a) Geometry of the masonry RVE (mm) (b) Mesh of the RVE 4.2 The discrete ceramic RVE As a second implementation of the proposed method, a discrete ceramic microscopic structure has been chosen. The material is linear, but a unilateral contact interface representing a potential crack is found in the middle of the geometry, as is shown in Fig. 4 G.A. Drosopoulos, M.E. Stavroulaki, K. Giannis, L. Plymakis, G.E. Stavroulakis and P. Wriggers 3. The mesh consists of rectangular plane stress elements, with out of plane thickness equal to 10mm. In addition, material properties have been taken from the literature: E = 328389N/mm2 , n = 0.22968. (a) (b) Figure 3: (a) Geometry of the discrete, ceramic RVE (b) Mesh of the RVE 4.3 Averaging procedure The proposed numerical scheme has two basic parts: first the derivation of the average stress of the RVE, thus the creation of the strain-stress database; second the estimation of the effective constitutive tensor, thus the strain-stiffness database. In the overall computational homogenization scheme the average stress will be the macroscopic stress and the effective elasticity tensor will be used for the development of the tangent stiffness matrix of the macroscopic model. The averaging relations are given here. To obtain these quantities, the subdomain integration, postprocessing capability of COMSOL was used.   1 1 m < ǫ >Vm = ǫ dVm , < σ >Vm = σ m dVm (2) Vm Vm Vm Vm The procedure was repeated in order to obtain the stiffness information for the final macroscopic model. For this reason, for every load strain vector of the microscopic analysis, three test, incremental strain vectors were considered (equation 3a), and three incremental average stress vectors were calculated (equation 3b), respectively. Then, the effective elasticity tensor was calculated by using Hooke’s law, according to equation 3c. This effective elasticity tensor is used in the overall homogenization scheme, for the estimation of the macroscopic tangent stiffness matrix. For plane stress conditions, the following relations are given: [δǫM ] = [δǫM 1 δǫM 2 δǫM 3 ] (3a) [δσ M ] = [δσ M 1 δσ M 2 δσ M 3 ] (3b) −1 [δσ M ] = CM [δǫM ] ⇒ CM = [δσ M ][δǫM ] 5 (3c) G.A. Drosopoulos, M.E. Stavroulaki, K. Giannis, L. Plymakis, G.E. Stavroulakis and P. Wriggers 5 THE MULTI-SCALE COMPUTATIONAL HOMOGENIZATION The final step of the proposed approach includes the development of a non-linear, multi-scale, computational homogenization scheme, for the investigation of the macroscopic structure. The simulation of an RVE in each Gauss point and each time step of the macro model, has been now replaced with the usage of the strain-stress and strainstiffness databases, which were previously created. Thus, instead of solving a non-linear finite element model of the microscopic problem in each Gauss point and time step, the databases and some interpolation method are used in order to obtain the macro stress and the stiffness of the macroscopic model, Fig. 4. The whole numerical scheme has been implemented with MATLAB. Plane stress, first order, full integration finite elements have been used in the macroscopic model. The Newton-Raphson incremental iterative procedure has been chosen, to capture the nonlinear behaviour of the masonry. Analysis starts by introducing a test strain vector. This is used as an input, for receiving from the stiffness database the information for the consistent stiffness of the macroscopic model. Then, for any current value of the macroscopic strain vector, a stress vector and an elasticity tensor are chosen from the databases. An interpolation method must be used, to obtain these quantities from the databases. In this work the MATLAB function ”TriScatteredInterp” is used, however other possible solutions for the creation of the metamodel (interpolation) can be used, for instance Neural Networks. By considering stress interpolation, each strain vector (3x1) corresponds to one average stress value. For the elasticity tensor (3x3), each strain vector corresponds to one value of the tensor. This tensor is the consistent stiffness of the Newton-Raphson incremental iterative procedure, used in the macroscopic model. Figure 4: The final multi-scale homogenization scheme 6 G.A. Drosopoulos, M.E. Stavroulaki, K. Giannis, L. Plymakis, G.E. Stavroulakis and P. Wriggers 6 6.1 RESULTS AND DISCUSSION: THE MASONRY STRUCTURE The masonry RVE In this section some results for the average stress - average strain relation and the failure of the masonry RVE will be presented. According to the diagrams of Fig. 5, the RVE simulation results in non-linear stress-strain behaviour. (a) (b) Figure 5: Average stress-strain diagrams obtained from the RVE analysis (a) direction x, (b) direction y. Loading along linear paths in the strain space defined by angles a and b The failure mode of some RVEs is shown in Fig. 6. According to this Figure, plastic strains are developed only in the mortar joints. Moreover, as the value of the parameters used to describe the gradually increased displacement loading are increased, the effective plastic strains given by COMSOL are also increased, from zero to a maximum value. (a) (b) (c) (d) Figure 6: Effective plastic strain from the RVE analysis, for gradually increased load, from (a) to (d) 6.2 Overall multi-scale homogenization The last step of the proposed approach is related to the development of an overall multiscale, computational homogenization scheme, in a FE2 sense, for the study of macroscopic masonry structures. To compare the results obtained from the proposed model, a second heterogeneous, macroscopic model is developed, for each masonry structure. MARC software has been used for this simulation of heterogeneous masonry structures, directly at 7 G.A. Drosopoulos, M.E. Stavroulaki, K. Giannis, L. Plymakis, G.E. Stavroulakis and P. Wriggers the macroscopic scale (Direct Numerical Simulation, DNS models). Thus, plane stress, first order, full integration elements have been used, with material properties and yield behaviour of the constitutive materials equal to the ones considered in the RVE analysis (Section 4). The first model which is presented here is a rectangular masonry wall, with dimensions equal to 1.82mx1.69m. Loading of this wall is a distributed vertical displacement equal to 5mm at the right vertical edge of the model, while fixed boundary conditions are applied to the left vertical edge of it. Fig. 7 shows the degradation of the strength of the structure. In particular, the plastic strain distribution is shown in Fig. 7b, for the direct heterogeneous macro model. The distribution of the trace of the elasticity tensor has been chosen as a qualitative only measurement of the degradation of strength, for the proposed multi-scale model, Fig. 7a. The dark blue colour shows bigger values of trace, while the light blue colour which gradually becomes red, smaller values, respectively. According to this Figure, the degradation of the strength obtained from the two models, has the same distribution in the domain. (a) (b) Figure 7: Degradation of the strength of a macroscopic masonry wall (1.82x1.69m) in the end of the analysis (a) proposed multi-scale homogenization (b) direct macroscopic simulation (MARC) Two more diagrams obtained from the two different methods are shown in Fig. 8. These diagrams represent the distribution of the vertical displacements along the top side of the masonry (Fig. 8a) and the distribution of the horizontal displacements, along the right side of it (Fig. 8b), respectively. The comparison between the two models leads to similar results, indicating that the proposed approach can be used for the simulation of non-linear, heterogeneous structures. In addition, in Fig. 8b the direct heterogeneous simulation presents a fluctuation, which is attributed to the alteration between brick and mortar material, in the model. In the multi-scale homogenization approach no such fluctuation is obtained, as the material is homogeneous. Another example of a big masonry wall with two openings (windows) is presented in the following lines. A distributed displacement loading of 5mm is applied to the top edge of the structure, while the bottom of it is fixed, Fig. 9a. Figs. 9b and 9c show that the image 8 G.A. Drosopoulos, M.E. Stavroulaki, K. Giannis, L. Plymakis, G.E. Stavroulakis and P. Wriggers (a) (b) Figure 8: Distribution of (a) vertical displacements along the top side (b) horizontal displacements along the right side, in the final time step of the degradation of the strength is the same for both the multi-scale homogenization and the direct heterogeneous model. In addition, a concentration of plastic strains appears in the corner of the windows, Fig. 9c. Finally, Fig. 9.d,e show the stress distribution obtained from the proposed homogenization approach and the direct macroscopic analysis. The comparison is satisfactory, regarding both the distribution of stresses and the limit values of them. 7 RESULTS AND DISCUSSION: THE DISCRETE CERAMIC MATERIAL Some results regarding the second RVE considered in this study, are presented here. Fig. 10 shows the discrete ceramic RVE and the two possible states of the unilateral contact interface: contact in Fig. 10(a) and no contact in Fig. 10(b). Calculation of the RVE includes the effect of contact nonlinearity, as it can be shown from the curve of Fig. 10 (c). A macroscopic rectangular geometry with dimensions equal to 0.5mx0.5m, tensile loading in the top edge and fixed displacement on the bottom edge have been initially chosen, for the overall implementation of the multi-scale analysis. Fig. 11 shows that similar displacement distribution is obtained from both the FE2 approach and the direct heterogeneous macroscopic model. Convergence of classical FE2 is certainly influenced by the nonsmooth unilateral contact nonlinearity and in some cases of compressive loading can not be achieved. This difficulty reminds us similar convergence behaviour in bimodulus elasticity and will be investigated in more details in the future. 8 CONCLUSIONS A multi-scale computational homogenization numerical scheme is presented in this article, for the study of non-linear masonry structures and discrete ceramic materials. The main idea is to replace the simulation of an RVE in the microscopic scale, with the usage of two databases carrying the information for the stress and the stiffness of the macroscopic structure. Parametric analysis of the RVE is initially conducted, in order to 9 G.A. Drosopoulos, M.E. Stavroulaki, K. Giannis, L. Plymakis, G.E. Stavroulakis and P. Wriggers (a) (c) (b) (d) (e) Figure 9: (a) The masonry wall with two openings (b) Degradation of the strength proposed multi-scale homogenization (c) Degradation of the strength - direct macroscopic simulation (MARC). Stresses Syy obtained (d) from the direct macroscopic model, (e) from the proposed homogenization approach (last time step) obtain this information. An overall homogenization scheme is then developed, in a FE2 sense. Comparison with a direct heterogeneous macroscopic model for the same masonry structures, shows that the proposed approach works well. According to the steps which have been followed, a general purpose structural analysis code is needed for the RVE analysis; a programming language is also necessary for the implementation of the overall homogenization and an interpolation method should be used for picking stress and stiffness values from the databases. The method is general, since it can be used in other masonry structures (arches), or other materials (composites). More complicated constitutive laws in the microscopic scale can be also considered. Finally, interpolation with Neural Networks can be used for the creation of the databases. All these investigations are left open for future research. 9 ACKNOWLEDGEMENTS The research project is implemented within the framework of the Action Supporting Postdoctoral Researchers of the Operational Program ”Education and Lifelong Learning” (Actions Beneficiary: General Secretariat for Research and Technology), and is co- 10 G.A. Drosopoulos, M.E. Stavroulaki, K. Giannis, L. Plymakis, G.E. Stavroulakis and P. Wriggers (a) (b) (c) Figure 10: The discrete ceramic RVE (a) contact (b) no contact conditions, (c) stressstrain curve of the ceramic RVE indicating the effect of contact nonlinearity (a) (b) Figure 11: Vertical displacements: (a) FE2 (b) DNS model financed by the European Social Fund (ESF) and the Greek State. REFERENCES [1] Sanchez-Palencia, E. Non-homogeneous media and vibration theory. Springer, Lecture notes in physics, (1980). [2] Zohdi, T.I. and Wriggers, P. An introduction to computational micromechanics. Springer, (2008). [3] Nguyen, V.P., Stroeven, M. and Sluys, L.J. Multiscale continuous and discontinuous modeling of heterogeneous materials: A review on recent developments. J. Multisc. Model (2011) 3:1–42. 11 G.A. Drosopoulos, M.E. Stavroulaki, K. Giannis, L. Plymakis, G.E. Stavroulakis and P. Wriggers [4] Dascalu, C., Bilbie, G. and Agiasofitou, E.K. Damage and size effects in elastic solids: A homogenization approach. Int. J. Solids Struct. (2008) 45:409–430. [5] Smit, R.J.M., Brekelmans, W.A.M. and Meijer, H.E.H. Prediction of the mechanical behaviour of non-linear heterogeneous systems by multi-level finite element modeling. Comput. Methods Appl. Mech. Engrg. (1998) 155:181–192. [6] Kouznetsova, V.G. Computational homogenization for the multi-scale analysis of multi-phase materials. Technical University Eindhoven, PhD thesis, (2002). [7] Drosopoulos, G.A., Wriggers, P. and Stavroulakis, G.E. Contact analysis in multiscale computational homogenization. Proceedings of 3rd international conference on computational modeling of fracture and failure of materials and structures (CFRAC), Prague, Czech Republic, (2013). [8] Drosopoulos, G.A., Wriggers, P. and Stavroulakis, G.E. Incorporation of Contact Mechanics in Multi-level Computational Homogenization for the Study of Composite Materials. Proceedings of 3rd international conference on computational contact mechanics (ICCCM), Lecce, Italy, (2013). [9] Miehe, C. and Koch, A. Computational micro-to-macro transitions of discretized microstructures undergoing small strains. Arch. Appl. Mech. (2002) 72:300–317. [10] Massart, T.J., Peerlings, R.H.J. and Geers, M.G.D. Structural damage analysis of masonry walls using computational homogenization. Int. J. Damage Mech. (2007) 16:199–226. [11] Hill, R. Elastic properties of reinforced solids: Some theoretical principles. J. Mech. Phys. Solids (1963) 11:357–372. 12