Influence of Boundary Conditions on the Simulation of a Diamond-Type Lattice Structure: A Preliminary Study
Emergent additive manufacturing processes allow the use of metallic porous structures in various industrial applications. Because these structures comprise a large number of ordered unit cells, their design using conventional modeling approaches, such as finite elements, becomes a real challenge. A homogenization technique, in which the lattice structure is simulated as a fully dense volume having equivalent material properties, can then be employed. To determine these equivalent material properties, numerical simulations can be performed on a single unit cell of the lattice structure. However, a critical aspect to consider is the boundary conditions applied to the external faces of the unit cell. In the literature, different types of boundary conditions are used, but a comparative study is definitely lacking. In this publication, a diamond-type unit cell is studied in compression by applying different boundary conditions. If the porous structure’s boundaries are free to deform, then the periodic boundary condition is found to be the most representative, but constraint equations must be introduced in the model. If, instead, the porous structure is inserted in a rigid enclosure, it is then better to use frictionless boundary conditions. These preliminary results remain to be validated for other types of unit cells loaded beyond the yield limit of the material.
Porous metals, also called metallic foams, are increasingly used in many different applications, especially in the aerospace and medical fields [1, 2]. The main advantage presented by these materials is the possibility of having their properties (stiffness, permeability, thermal resistance, etc.) adjusted by controlling the porosity level and the morphology (spatial distribution, size, shape, connectivity, etc.) of the pores. Simulations based on numerical techniques such as finite elements are generally used in the design process to determine the optimal porosity and morphology of the foam. Unfortunately, numerical simulations using the complete and detailed geometry of the entire porous structure become impracticable in terms of computational cost as soon as the studied geometry comprises a large number of pores. A multiscale modeling approach, otherwise known as a homogenization technique, must then be used , and the behavior of a porous volume is assumed to be similar to that of a fully dense homogeneous volume having fictive properties (Young’s modulus, yield stress, thermal resistance, etc.). This strategy is frequently used for heterogeneous materials such as composites or foams. For a foam material having a given morphology, these fictive properties are related to the porosity level through the use of scaling relations , which are generally obtained by performing finite element simulations on a relatively small volume of porous material. If the foam has a random architecture, a volume called a “representative volume element” is used in the simulations, and it must be sufficiently large to be representative of the foam. If that is not the case, the foam will have ordered architecture, and only the “unit cell” can be analyzed. The unit cell is the smallest geometry that is copied in several directions to obtain the ordered porous structure.
Performing a simulation with a unit cell implies that boundary conditions must be defined at its bounding surfaces. Several types of boundary conditions can be selected, and the behavior of the unit cell will be greatly influenced by this choice . It is thus important to choose boundary conditions that are truly representative of the behavior of the real cell surrounded by many other identical cells in order to build reliable scaling relations. In this work, four types of boundary conditions are considered: free, frictionless, planar, and periodic. They are schematized in Figure 1 in a two-dimensional representation. In all the cases, the initially undeformed unit cell is compressed by imposing a uniform displacement of the upper bound, while the lower bound is restrained only in the direction of this displacement. The difference between the studied cases lies in the restrictions that are imposed on the other bounds, namely, the two other vertical faces in Figure 1 (note that, in 3D, there are four vertical faces).
The free boundary condition is the simplest one to implement, and, in it, the vertical faces can deform freely without any restriction. This boundary condition does not take the presence of neighboring cells into consideration. Nevertheless, Öchsner and Mishuris  did use it to simulate the behavior of a cell located on the periphery of a sample consisting of a large number of cells. The frictionless boundary condition is also very simple to implement because it is assumed that vertical faces are forced to remain on their initial plane, while the displacements on this plane occur without friction . The planar boundary condition is based on the assumption that vertical faces always remain flat and parallel during the deformation of the structure [4, 6, 8–11]. Finally, the periodic boundary condition forces two opposite faces to deform in exactly the same manner [12–17]. The implementation of the periodic boundary condition generally requires that constraint equations be defined in order to link the displacement of a node to that of the corresponding node on the opposite face [18, 19]. Unfortunately, comparative analyses studying the effect of these four types of boundary conditions on a unit cell are lacking. Öchsner et al.  carried out such a comparative study, but only two different boundary conditions were considered.
The objective of this paper is to study the effect of the boundary conditions on the validity of a simulation performed on a unit cell. From a series of numerical simulations, we will determine the boundary condition that gives the best results, depending on the type of symmetries (reflectional or translational). Firstly, the unit cell used in this work is described. Secondly, a mesh convergence analysis is carried out to determine the optimal size of the elements that do not influence the calculated results. Then, a structural convergence analysis is presented in order to determine the reference behavior of a unit cell located in the middle of a sample comprising many other similar cells. Finally, the results obtained with different boundary conditions applied on the unit cell are compared to those obtained with the reference cell.
2. Materials and Methods
2.1. Presentation of the Diamond Unit Cell
A diamond-type unit cell consisting of 12 struts of equal length is employed in this work. The red lines in Figure 2 represent the wireframe schematization of the struts. The dimensions of the rectangular box that encompasses the unit cell are . Dimensions and are dependant, according to (1). The rectangular bounding box has 6 faces called front, back, right, left, top, and bottom. This unit cell is chosen in order to have two types of symmetries (reflectional and translational) that should be considered differently . Indeed, opposite pairs of top-bottom and front-back faces are reflectional symmetries. This means that the adjacent cell can be obtained by reflecting the unit cell with respect to a plane normal to the - and -axes (see Figure 2 for the -, -, and -axes). The third pair of faces (right-left) is a translational symmetry and the adjacent cell is obtained by translating the unit cell along the -axis.The porosity level of the unit cell is controlled by varying the strut diameter. In Figure 3, two unit cells having the same strut length , but with different diameters , are shown. On the left, the unit cell has a ratio and a porosity of 91.5%, while, on the right, a ratio leads to a porosity of 57.7%. In this work, the strut length and strut diameter are kept constant at mm and mm, which gives a ratio equal to 1.25 and a porosity of 55.1%.
The unit cell could be further simplified by using only one-quarter of the bounding box shown in Figure 1. Four struts would then be required to represent the unit cell, and it would be impossible to apply periodic boundary conditions due to the absence of corresponding nodes on opposite faces. Therefore, the geometry shown in Figure 3 is indeed the smallest diamond-type unit cell having three identical pairs of faces (top-bottom, front-back, and right-left).
2.2. Mesh Convergence
If the struts shown in Figure 2 were very slender, a geometry having a very high porosity would be obtained, and it could be meshed with beam elements. To be valid, it is recommended that elements based on Timoshenko beam theory should have a slenderness ratio greater than 30 . The slenderness ratio is defined as where(i) is the shear modulus of the material (),(ii) is the area of the cross-section (),(iii) is the length of the member (),(iv) is Young’s modulus of the material,(v) is the moment of inertia of the cross-section ().For the studied geometry, the slenderness ratio is approximately 9.6. The formulation of beam elements is no longer valid, and the use of volume elements becomes unavoidable. The aim of the mesh convergence analysis presented in this section is to determine the optimal size of the elements below which the results become independent of this size.
Meshing of the unit cell is carefully executed in order to have an identical tessellation (i.e., the same number and position of nodes) on each pair of opposing faces (front-back, top-bottom, and right-left). For example, the front face is first meshed with shell elements, and then the mesh is copied to the back face. Afterwards, the unit cell is meshed with volume elements and the software uses the shell element to build the volume elements. Finally, shell elements are removed, and only volume elements remain.
The unit cell is first meshed with a fairly large element size , and a displacement is imposed on the nodes of the front and back faces to simulate a 0.1% compression of the unit cell along the -axis (see Figure 2 for the representation of the coordinate system). A free boundary condition is applied to the other four faces (top, bottom, right, and left). Modeling, meshing, solving, and postprocessing are carried out with ANSYS Mechanical APDL 15.0 finite element software. The geometry is meshed with SOLID187 10-node tetrahedral elements, and a linear elastic material model is chosen ( GPa, ). The simulation is then repeated with an element size twice as small as that of the previous iteration. As shown in Figure 4, reducing the element size increases the ratio as well as the number of nodes in the mesh. This optimal element size is determined when the results at a given iteration show variations that are less than 1% of those of the previous iteration.
2.3. Structural Convergence
The purpose of the structural convergence analysis is to determine the behavior of a reference cell surrounded by other identical cells. When a sample comprising a large number of cells is compressed, all cells behave similarly, except for those that are located at the periphery of the sample (the edge effect). The idea is to successively increase the number of cells in the sample and extract the results for the middle one, which requires samples comprising an odd number of cells in each direction (, , , etc.). When the behavior of the central cell no longer changes, it can be assumed that the reference behavior is obtained and that it is not influenced by the boundaries of the sample. Figure 5 shows the samples that are used for structural convergence.
Each sample is compressed similarly as in the previous mesh convergence analysis. For example, the displacement along the -axis (front and back faces) for the sample is five times larger than that for the sample. For the other boundary faces (right, left, top, and bottom), two different boundary conditions, called “boundary condition of the structure” or “sBC,” are used: free or frictionless. The free sBC simulates a sample that can deform unrestrictedly at its boundaries, while the frictionless sBC simulates a sample inserted in a completely rigid enclosure.
2.4. Effect of the Boundary Condition Applied to the Unit Cell
Knowing the behavior of a reference unit cell surrounded by other identical cells for two different sBC, the behavior of a single unit cell is then analyzed, but using different boundary conditions (free, frictionless, planar, or periodic). The choice of this boundary condition is referred to as the “boundary condition applied to the unit cell” or “ucBC.” The objective is to determine which one most accurately represents the behavior of the reference cell. For the back and front faces, the unit cell is compressed by fixing a uniform displacement along the -axis of all the nodes resting on these faces, and the following assumptions are made, depending on the choice of the ucBC:(i)Free, frictionless, or planar: all the nodes on the back and front faces can move freely along the two other directions (-axis and -axis).(ii)Periodic: constraint equations are added to guarantee that the displacements in the - and -directions are the same for each pair of corresponding nodes on the back and front faces.For the other four faces (top, bottom, right, and left), the following actions are taken, depending on the choice of the ucBC: (i)Free: all the nodes on these four faces can move freely (no constraint), except for a few fixed displacements, to prevent rigid body motion.(ii)Frictionless: the -displacement of the nodes lying on the top and bottom faces, as well as the -displacement of the nodes lying on the right and left faces, are set to zero. No friction is considered, and, as a result, these nodes can move freely in the other directions.(iii)Planar: constraint equations are defined such that the -displacement of the nodes resting on the top face is identical. The top face then remains planar when the structure is deformed. Similarly, constraint equations are defined for the -displacement of the nodes on the bottom face and for the -displacement of the nodes on the left and right faces. No friction is considered, and, as a result, these nodes can move freely in the other directions.(iv)Periodic: constraint equations are defined for each pair of corresponding nodes of opposite faces (top-bottom or right-left) such that the displacement is the same in the three directions. The number of constraint equations is significant. Also, precautions must be taken for the nodes resting simultaneously on more than one face in order to avoid conflictual definitions of constraint where a degree-of-freedom is simultaneously a master and a slave.The resulting displacement of each node of the unit cell is compared to that of the same node in the reference cell. A global deviation is then calculated according to (2), where in the total number of nodes in the unit cell; , , and are the components of the displacement of the th node of the unit cell; and , , and are, respectively, the same values for the same nodes in the reference unit cell. If a unit cell is deformed exactly in the same way as the reference cell, the deviation would then be equal to zero. The less representative the boundary condition is, the higher the deviation is.Two other deviations are also calculated. Instead of using all the nodes of the unit cell, only those resting on the top and bottom faces are considered. The deviation calculated with only these nodes gives an indication of the representability of the different boundary conditions for faces involved in a reflectional symmetry. Similarly, the deviation is calculated with only the nodes resting on the right and left faces for the representability of a translational symmetry.
3. Results and Discussion
3.1. Mesh Convergence
To determine the optimal size of the tetrahedral elements, the unit cell is meshed with different element sizes, ranging between 1.6 mm ( nodes) and 0.05 mm ( nodes). The maximum displacement, the maximum Von Mises stress, and the total strain energy are calculated as reported in Table 1 and graphically plotted in Figure 6. The percentages given in parentheses in Table 1 are the variations between successive iterations. For example, for an element size of 0.8 mm, a variation of 18.0% is obtained for the maximum stress when compared to the same metric calculated with an element size of 1.6 mm. The Von Mises stress distribution is plotted in Figure 7, and it is obvious that stress raisers are generated at the sharp edges at the junction of struts.
Convergence is deemed to be reached when a reduction of the element size by half causes a variation in the results of less than 1% compared to the previous iteration. Therefore, from the results of Table 1, convergence in terms of maximum displacement and strain energy is achieved for an element size of 0.2 mm (). For the maximum stress, it is obvious that there are singularities at the sharp edges, and it is normal that an increasing stress is calculated for smaller elements since this stress is theoretically infinite at these locations; that is why the maximum stress should not be used in a convergence analysis where local stress raisers are present, and another global measurement such as strain energy should rather be employed . For the following analyses, an element size of 0.2 mm will be used.
3.2. Structural Convergence
The exact same mesh comprising 36,995 elements of 0.2 mm is used to mesh each cell. Thus, for the sample consisting of 125 cells, a mesh comprising more than 4.4 million nodes is required. Since the finite element mesh is identical for every cell, the element size is no longer an issue for comparison purposes, especially near stress raisers. Therefore, maximum stress and maximum displacement are the two metrics that are used to study the structural convergence. Results for the free and frictionless sBC are, respectively, given in Tables 2 and 3.
It is obvious that the results are already converged for the sample. Indeed, when using a sample, variations of less than 1% are obtained when compared to the sample. Therefore, the behavior of the cell located in the middle of the sample could be used as the reference unit cell; however, since the behavior of the sample is already calculated with the same number of nodes, it is the cell at the middle of the sample that is considered as the reference unit cell.
3.3. Effect of the Boundary Condition Applied to the Unit Cell
Figures 8 and 9 report the three deviations calculated for the four ucBC, when the behavior is compared to reference cell obtained with the free and the frictionless sBC, respectively.
On the one hand, when a porous ordered structure is compressed with its boundaries unrestrained (free sBC, Figure 8), it is obvious that the periodic ucBC is the one that is the most representative, and this is from a global point of view, as well as for both types of symmetries (reflectional and translational). The frictionless ucBC is the least representative with deviations that are 70 to 80 times larger than those of the periodic ucBC. In between those two boundary conditions, the planar ucBC and free ucBC offer similar deviations that are approximately one order of magnitude larger than those of the periodic ucBC. One exception concerns the use of a planar boundary condition on a face representing a reflectional symmetry. Also, for all the boundary conditions, deviations for faces representing translational symmetries were systematically larger than those for faces representing reflectional symmetries. On the other hand, if a porous ordered structure is compressed inside a rigid enclosure (frictionless sBC, Figure 9), the conclusions are completely different and the frictionless ucBC becomes the most representative.
The objective of this paper was to determine the type of boundary conditions to be applied on the bounding faces of a unit cell in order to obtain a behavior that is representative of a cell surrounded by many other similar cells. Mesh and structural convergence analyses were first carried out after which four types of boundary conditions were tested. By comparing the behavior of the unit cell with that of the reference cell, it was found that the most representative boundary condition that should be applied on the unit cell depends on how the boundaries of the entire structure are restrained. If they can move freely, then a periodic boundary condition is the most representative, but if, instead, they are fully restrained, then a frictionless boundary condition is more accurate.
If the boundaries of a structure can freely deform, but its geometry is not conducive to the use of periodic boundary conditions, then planar boundary conditions should be employed instead of free boundary conditions, even if both choices give similar deviations (see Figure 8). This affirmation is based on geometrical compatibility considerations. If the deformed geometry of a unit cell with free boundary condition is patterned, geometrical discontinuities can be observed. Indeed, if the deformed geometry with the free boundary condition shown in Figure 1 is horizontally copied, it is obvious that incompatibilities will be observed at the boundaries. With a planar boundary condition, compatibility is guaranteed since the faces always remain flat.
The conclusions stated above are valid for a 55.1% porosity diamond-type unit cell behaving according to a linear material law and subjected to small strains. In future works, additional numerical investigations will be carried out to verify if these conclusions remain valid for other types of unit cells showing reflectional and translational symmetries, for other porosities, or when nonlinear material behavior is considered, including plasticity and failure.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
This work has been funded by National Science and Engineering Research Council of Canada (NSERC) through the Discovery Grant of Professor Terriault aimed at the modeling of porous materials.
E. Del Olmo, E. Grande, C. R. Samartin et al., “Lattice structures for aerospace applications,” in Proceedings of the 12th European Conference on Spacecraft Structures, Materials and Environmental Testing, Noordwijk, Netherlands, 2012.View at: Google Scholar
G. W. Rodgers, E. E. Van Houten, R. J. Bianco, R. Besset, and T. B. F. Woodfield, “Topology Optimization of Porous Lattice Structures for Orthopaedic Implants,” in Proceedings of the 19th World Congress of the International Federation of Automatic Control, Cape Town, South Africa, 2014.View at: Google Scholar
C. Simoneau, P. Terriault, J. Rivard, and V. Brailovski, “Modeling of metallic foam morphology using the representative volume element approach: development and experimental validation,” International Journal of Solids and Structures, vol. 51, no. 21-22, pp. 3633–3641, 2014.View at: Publisher Site | Google Scholar
G. Maîtrejean, P. Terriault, and V. Brailovski, “Density dependence of the superelastic behavior of porous shape memory alloys: representative volume element and scaling relation approaches,” Computational Materials Science, vol. 77, pp. 93–101, 2013.View at: Publisher Site | Google Scholar
R. Glüge, “Generalized boundary conditions on representative volume elements and their use in determining the effective material properties,” Computational Materials Science, vol. 79, pp. 408–416, 2013.View at: Publisher Site | Google Scholar
A. Öchsner and G. Mishuris, “Modelling of the multiaxial elasto-plastic behaviour of porous metals with internal gas pressure,” Finite Elements in Analysis and Design, vol. 45, no. 2, pp. 104–112, 2009.View at: Publisher Site | Google Scholar
J. Ramírez, M. Cardona, J. Velez et al., “Numerical modeling and simulation of uniaxial compression of aluminum foams using FEM and 3D-CT images,” Procedia Materials Science, vol. 4, pp. 227–231, 2014.View at: Publisher Site | Google Scholar
M. Panico and L. C. Brinson, “Computational modeling of porous shape memory alloys,” International Journal of Solids and Structures, vol. 45, no. 21, pp. 5613–5626, 2008.View at: Publisher Site | Google Scholar
B. F. Oliveira, L. A. B. Da Cunda, and G. J. Creus, “Modeling of the mechanical behavior of metallic foams: damage effects at finite strains,” Mechanics of Advanced Materials and Structures, vol. 16, no. 2, pp. 110–119, 2009.View at: Publisher Site | Google Scholar
G. Maîtrejean, P. Terriault, D. Devís Capilla, and V. Brailovski, “Unit cell analysis of the superelastic behavior of open-cell tetrakaidecahedral shape memory alloy foam under quasi-static loading,” Smart Materials Research, vol. 2014, Article ID 870649, 11 pages, 2014.View at: Publisher Site | Google Scholar
A. Öchsner and K. Lamprecht, “On the uniaxial compression behavior of regular shaped cellular metals,” Mechanics Research Communications, vol. 30, no. 6, pp. 573–579, 2003.View at: Publisher Site | Google Scholar
T. Kanit, S. Forest, I. Galliet, V. Mounoury, and D. Jeulin, “Determination of the size of the representative volume element for random composites: statistical and numerical approach,” International Journal of Solids and Structures, vol. 40, no. 13-14, pp. 3647–3679, 2003.View at: Publisher Site | Google Scholar
Z. Xia, C. Zhou, Q. Yong, and X. Wang, “On selection of repeated unit cell model and application of unified periodic boundary conditions in micro-mechanical analysis of composites,” International Journal of Solids and Structures, vol. 43, no. 2, pp. 266–278, 2006.View at: Publisher Site | Google Scholar | MathSciNet
L. T. Harper, C. Qian, T. A. Turner, S. Li, and N. A. Warrior, “Representative volume elements for discontinuous carbon fibre composites—part 1: boundary conditions,” Composites Science and Technology, vol. 72, no. 2, pp. 225–234, 2012.View at: Publisher Site | Google Scholar
V.-D. Nguyen, E. Béchet, C. Geuzaine, and L. Noels, “Imposing periodic boundary condition on arbitrary meshes by polynomial interpolation,” Computational Materials Science, vol. 55, pp. 390–406, 2012.View at: Publisher Site | Google Scholar
S. M. Mirkhalaf, F. M. Andrade Pires, and R. Simoes, “Determination of the size of the representative volume element (RVE) for the simulation of heterogeneous polymers at finite strains,” Finite Elements in Analysis and Design, vol. 119, pp. 30–44, 2016.View at: Publisher Site | Google Scholar
S. E. Stapleton, L. Appel, J.-W. Simon, and S. Reese, “Representative volume element for parallel fiber bundles: model and size convergence,” Composites Part A: Applied Science and Manufacturing, vol. 87, pp. 170–185, 2016.View at: Publisher Site | Google Scholar
S. Das, A. Maroli, and N. Neithalath, “Finite element-based micromechanical modeling of the influence of phase properties on the elastic response of cementitious mortars,” Construction and Building Materials, vol. 127, pp. 153–166, 2016.View at: Publisher Site | Google Scholar
R. Wang, L. Zhang, D. Hu et al., “A novel approach to impose periodic boundary condition on braided composite RVE model based on RPIM,” Composite Structures, vol. 163, pp. 77–88, 2017.View at: Publisher Site | Google Scholar
A. Öchsner, W. Winter, and G. Kuhn, “On an elastic-plastic transition zone in cellular metals,” Archive of Applied Mechanics, vol. 73, no. 3-4, pp. 261–269, 2003.View at: Publisher Site | Google Scholar
S. Li, “Boundary conditions for unit cells from periodic microstructures and their implications,” Composites Science and Technology, vol. 68, no. 9, pp. 1962–1974, 2008.View at: Publisher Site | Google Scholar
ANSYS Documentation—Release 15.0, "Beam189 Element Description", ANSYS Inc., 2013.
W. Frei, “How to identify and resolve singularities in the model when meshing,” Comsol Blog, 2013, http://www.comsol.com/blogs/how-identify-resolve-singularities-model-meshing/.View at: Google Scholar