Numerical modelling of the rigid polyurethane foam microstructure

The paper presents a numerical model of 2-D microstructure of rigid polyurethane foam. The selection of geometric and material parameters is presented. For a particular structure, its behavior has been studied for typical cases of external loads (or forced displacements). The characteristic phenomena have been identified and described. A parametric analysis was performed due to the dimension of the cross-section of the struts which form the cell edges. An analysis of the impact of support and loading conditions on the behavior of the cell structure was performed.


Introduction
Rigid cell structures of low density are very common in nature and have been produced by humans for many years. The relative density (the density of the cellular material divided by that of solid from which the cell walls are made) of typical rigid foam structures is in the range 0.025-0.05 and the Young modulus of the foam is about 1-5 MPa. Many believe that the importance of this class of structures is much greater than, for example, fiber composites. Despite this, the mechanics of rigid foam structures is often neglected.
The problem of deformation of cellular structures has been repeatedly taken up in the literature. The starting point is usually the study of two-dimensional cellular solids [1]. Describing phenomena taking place in these structures, we can distinguish in principle elastic behavior, elastic buckling and formation of plastic hinges. One of the fundamental determinants of cell structure and its behavior is the relative density of the foam [2]. Searching for relationships between the relative density factor and material properties was widely discussed in the literature [3,4], while the analysis of spatial structures is far more complex than 2-dimensional systems. The principal mechanism of linear-elastic deformation (cell wall bending) was correctly identified by Menges & Knipschild [5]. Much later it turned out that the density of the foam is not sufficient for a reliable characterization of the foam structure. The effect the manufacturing process on the foam microstructure and the properties of the matrix material is important but has received limited attention in the literature [6]. This paper presents a relatively simple model of a rigid cellular foam structure. The elementary problem of the influence of the structure stiffness and density on the macroscopic behavior is discussed. For typical schemes of external loads, the phenomena leading to the mechanism of cells damage are observed. Analyzing these issues (including numerical issues) is essential for further, in-depth studies on phenomena occurring in foam cores of sandwich panels.

Problem formulation
A two-dimensional, honeycomb structure consisting of regular hexes is considered (Fig. 1). Such structure has been repeatedly described [1,2,7] and recommended for preliminary analysis of the behavior of foam cores. Interestingly, it was mainly used to find the relationship between foam density and foam mechanical parameters. This structure corresponds roughly to the actual foam structure (Fig. 2), although of course the actual cells are of varying size and shape. Another issue is the presence of shells in closed cell structures (Fig. 2). Despite these shells, many authors emphasize that their effect on the strength of the cell structure is small because the shells rapidly lose stability. In addition, the edges of the cells are much stiffer than the shells are.
The structure consisting of 60×60 = 3600 open-cell hexagons was analyzed in the paper. This is a typical honeycomb structure. The length of the edge of the hex is 0.1 mm, so the whole structure is 10.48×9.05 mm. It is worth noting that this system is not symmetrical due to the small displacement of neighboring cells and the even number of cells. This was introduced consciously to avoid the symmetry of the numerical problem. The aim of the study is to analyze the behavior of the system under the influence of unidirectional compression, tension and shear. For this purpose, 3 groups of models were created in the Abaqus system environment. The appropriate boundary conditions were assumed. In addition, the impact of structural rigidity on the cell deformations and stress values was investigated.

Numerical models
The structure shown in Fig. 1 (60×60 cells) has been modeled using two-dimensional B23 beam elements (two-node element, cubic shape functions). Each hexagon edge consists of 3 finite elements. All connections of elements are rigid, but of course the nodes can move and rotate. Based on the available literature, it is assumed that the structure is made of the elastic material characterized by E s = 1200 MPa and  s = 0.3. Of course, Poisson's ratio of the material does not affect the deformations observed macroscopically. Although only the elastic material was considered in the first models, a simple definition of plasticity was finally adopted. It was assumed that at the stress of 127 MPa the plastic strain is 0.0, and at the stress 140 MPa plastic strain is 0.1.
One of the challenges was determining the dimensions of the cross-section of the beam elements from which the structure is composed. From the macroscopic observations (stressstrain relation), the Plateau effect is initiated at a stress level of 90-100 kPa, which corresponds to a strain of 2-2.5%. Based on this information, and additionally taking into account that some deformations result from the shortening of the elements due to compression, a square cross-section of 0.015×0.015 mm was assumed. These dimensions are of a similar order as the values observed in reality (cf. Fig. 2). Finally, the structure consists of beams with a square cross-section of dimension t s = 0.015 mm, which form regular hexes connected to a system of 60×60 cells. Analyses were also carried out for different cross-sections, namely for square sections of dimension 0.0125 mm and 0.010 mm.
Due to boundary conditions, three basic cases were distinguished (compression, tension, shear). In certain areas, the appropriate displacement was forced, which quite well reflect the conditions of real experiments. On the bottom surface, displacements in both directions have been blocked (x, y). The upper surface had the following boundary conditions: a) compression: blocked horizontal displacement, forced vertical displacement downward, b) tension: blocked horizontal displacement, forced vertical displacement upward, c) shear: blocked vertical displacement, forced horizontal displacement to the right. In models 6, 7 and 8 (with a note: c -compression, t -tension, s -shear), these conditions were applicable to all nodes located on the respective surfaces of the model. Models 6, 7 and 8 correspond to the structures made of the beam elements with the cross-section dimension 0.015 mm, 0.0125 mm and 0.010 mm, respectively.
In order to analyze the other effects, 3 other models (9c_a, 9c_b, 9c_c) were also considered, wherein the starting point was model 6c. In the 9c_a model, instead of limiting the horizontal displacement (x-direction) of all the nodes on the two external surfaces, such a restriction was introduced only at one central point on the surfaces. As a result, the other nodes were free to move horizontally. In the 9c_b model, the horizontal displacement restriction was applied to all nodes located on two faces (as in the model 6c), but the vertical constraint was only applied to 6 points in the center of the top surface. Model 9c_c was different from model 6c in such a way that compression was forced in another direction (in x instead of y-direction). The displacements of the nodes on the left surface were blocked, and horizontal displacements (leftward) of all the extreme nodes on the right surface were forced.

Results
For each considered model, a numerical calculation was carried out until a loss of convergence of the solution. The loss of convergence occurred with a specified displacement. In the case of compression and tension, these were displacements in the ydirection (u 2 ), whereas in the case of shear the displacements were in the x-direction (u 1 ). The characteristic results, corresponding to the moment when the solution convergence was lost, are shown in Tab. 1.
Analyzing the results presented in Tab. 1, it should first be noted that the numerical problems lost the convergence at similar values of the compressive normal stress (of the order of 140 MPa). The stress level is equivalent to the yield stress of the material, but a very similar result was obtained for a material defined as ideal elastic. This means that the loss of convergence is rather related to the magnitude of the deformation of the elements. This is illustrated in Fig. 3, where the deformations of the entire specimens during compression (model 6c), shear (6s) and tensile (6t) are presented. The symbol S11 denotes normal stress in the beam element. It is worth noting that the extreme deformations occur in farthest cells located at the bottom and top of the specimens.  The extreme values of support reactions RF2 (Tab. 1) are very different and obviously decrease considerably as the beam section size t s decreases. The mean values of the RF2 (reaction in y-direction) are much more interesting. These values are practically equivalent to the mean value of normal force in vertical pillars of the cellular structure. In the case of model 6c, the normal force in the pillar is 0.0015 N, which at the cross-section area of 2.25×10 -4 mm 2 means the normal stress of 6.67 MPa. This means that the bending of the structural elements is decisive for the level of normal stress and deformation of the system. It is easy to see this by comparing the strain in the pillar caused by the axial compression: with the deformation of the structure observed macroscopically: Assuming a certain thickness of the structure, e.g. b s = 0.1 mm, we can also calculate the macroscopically observed stress level and (effective) modulus of elasticity: If the density of the material from which the structure is made is assumed in accordance with literature [2,3] as  s = 1200 kg/m 3 , the above values are obtained for a foam structure with a density of 30.96 kg/m 3 . When we reduce the cross-sectional dimension of the beam component of the structure t s (models 7c i 8c), the density, structural rigidity and effective E-modulus decrease. These results provoke some reflections. First of all, the rigidity of the structural elements is too low. Certainly the E s modulus can be taken at the limit of maximum values known from the literature (1600 MPa). Moreover, the cross-sectional dimension of the beam element can be slightly increased to obtain a specimen density of 40 kg/m 3 , which is consistent with typical values for structural applications. It is much more difficult to achieve the correct (geometrically and physically) nonlinear behavior of the structure and to obtain well-known mechanisms of destruction.
Analyzing deformations of a structure (including cellular) it is difficult not to refer to the issue of Poisson's ratio. In the case of cellular structure, a question immediately arises on how to determine and interpret the ratio. This is not the value set for one cell. Rather it is a value averaged over a certain area. Looking at the deformations of the entire structure in model 6c, it can be seen, however, that it is difficult to define the area for which we should determine Poisson's ratio. Boundary conditions (limitation of displacements u 1 ) significantly disturb deformations of the structure. If you rely on the ratio of the maximum expansion of the specimen to its shortening, then we will get the value  = 0.985. To test how boundary conditions affect deformations of the system, in model 9c_a the compression was induced, but the horizontal displacements were limited only in one point of the lower and upper external surface (Fig. 4). The system deforms more regularly, however, measurements of the displacements of the entire specimen give  = 0.834.
The regular hexagonal cell systems have different properties if we apply load (or forced displacement) in another direction. The anisotropy of the structure highly influences the structural behavior of mechanical systems [8]. In the 9c_c model, the boundary conditions have been changed. These are analogous to those in model 6c, except that the limitation of displacements u 1 and u 2 refers to the left surface (instead of the lower one), while on the right-hand side, displacements u 2 are blocked and the compression is induced horizontally. The deformed system was slightly less oval than in the case of the model 6c, but the Poisson ratio determined using the macroscopic extreme displacements of the specimen was equal to  = 0.802. It is worth noting that all the obtained Poisson's ratios are higher than 0.5, which is greater than the theoretical limit for natural materials. Moreover, it is much higher than the values given in the literature for cellular structures. An interesting issue is also the impact of concentrated load on deformations of cell structures. Model 9c_b was developed based on model 6c, except that the vertical displacement (compression) was imposed not on the entire upper surface, but only on the six nodes located near the center of the upper surface. Despite the considerable rigidity of the system, surprisingly rapid vertical displacements were achieved at the place of extortion. These displacements occurred despite the fact that the deformations of individual cells were not too large (Fig. 5). In addition, the concentrated extortion was very quickly redistributed to a large area of the structure.

Conclusions
The paper presents a simple 2-dimensional cellular structure. This structure was modeled using beam elements. Linear elasticity and plasticity have been assumed for the material from which the structure is constructed. The issues of the influence of stiffness of structural elements and boundary conditions (support and load) on structural behavior were discussed. The problem of Poisson's ratio determination for the structure was discussed.
The results of the conducted analyses revealed surprisingly large local deformations of cells. These cells are degraded, leading to a loss of convergence of the numerical solution. Based on the deformations achieved for different structures, it can also be stated that the stiffness of the elements from which the structure is made is too low. For rigid foams used in civil engineering (e.g. used as cores in sandwich panels), a higher stiffness of beam elements should be assumed. On the other hand, the structural elements should lose stability much earlier, which is in some contradiction to the earlier conclusion. The probable solution to this problem is to use a model with a different geometry, different edge geometry, or a spatial model. This latter solution is appropriate after a slightly deeper understanding of the problems encountered in 2-D structures.
The intriguing issue is the determination of the Poisson ratio for the foam structures. The results obtained for the honeycomb structure are very high and do not confirm the observations made for rigid polyurethane foams. The redistribution of local loads also deserves attention. This issue is important from an engineering point of view [9]. The results presented in the paper indicate rapid redistribution of forces and surprisingly small deformations of individual cells.