Application of Substructure Techniques to Syntactic Metal Foams in a Finite Element Environment

The presented work focuses on the development of a novel method that can numerically describe the properties of metal matrix syntactic foam (MMSF) with low memory requirements and short computational times without losing the properties of the interior structure. In this paper, we propose a novel method that avoids using the homogenization technique and instead rearranges stiffness matrices and constructs specific substructures to perform the overall construction. The two-dimensional cases are discussed in order to focus on the methodology itself. First, the reductions and structural design with solid mesh structures were performed


Introduction
One of the main problems with the modelling of metallic foams lies in their stochastic structure. The solutions for this problem are: • to investigate representative elements of the cell structure, or • to model the whole random foam structure.
The first case is a strong simplification that may lead to a poor agreement with the measured properties, especially in the case of more complex structures like metal matrix syntactic foams (MMSFs), in which the porosity is ensured by embedded porous second-phase particles. In the second case, the calculation time can reach days even on high-performance computers. In this paper a new, so-called substructure technique is proposed, to solve the task and problems mentioned before. The main goal of this work is to create an accurate, reduced-order model. This way the running time can be drastically decreased, which is a common problem in modelling of metal foams. This paper is based on mechanical models, but it is also important to briefly describe models based on other principles (e.g., numerical calculations) to get a more complex picture of the modelling methods in the field of MMSFs.
Besides the low density of metal foams, a very favorable property is their high energy absorption capacity. It is difficult to formulate generalized properties to describe the behavior of the metal foams [1]. Their elastic modulus is a property that allows the foams to be modelled relatively well, within certain limits. Several parametric modelling techniques have been developed that can give good estimates of compressive curves and elastic modulus. From previous studies, it can be seen that some papers process and compare several different modelling techniques (e.g., numerical, finite element, analytical) [2][3][4][5]. The most commonly used modelling process is the homogenization technique, which is often combined with the representative volume elements (RVE) [5][6][7][8][9]. This method has advantages but disadvantages too. In many cases, the models do not even consider the interior foam structure but assign foam properties to compact models. This method is commonly found in commercial finite element software. Neglecting or simplifying the marrow structure often causes differences, to a lesser or greater extent, between calculated values and experimental results [10]. Some researchers have found that preserving the interior defects and accurately representing the disordered interior structure lead to significantly improved results. [11][12][13]. In any cases, CT-based models are prepared for finite-element calculations, which often show exact agreement with experimental results, but their extremely time-consuming preparation makes calculations with larger-scale models nearly impossible [8,13]. However, CT scanning is often used on its own as an examination of the response to external loads [9,[14][15][16]. Occasionally special, more manageable interior structures are created to try to replicate the interior assembly behavior. However, these usually only give good results for certain load cases. The capabilities of the different methods are often presented on two-dimensional structures and in the case of satisfying results, the methods are extended to spatial models [17,18]. A further step forward from the use of conventional metal foams is the use of MMSFs, where the foam structure is provided by filling with hollow spheres. With MMSFs better properties can be achieved than with conventional closed-cell foams, but their accurate modelling is a more difficult task (due to the second phase filler) [19]. In this case, modelling the chemical interaction at the interface between the hollow spheres and the matrix becomes an important factor. The calculations will be greatly improved if the interfacial separation can be simulated [11]. Another important factor is the distribution of hollow spheres in the matrix. It has been shown that the Gaussian distribution provides the closest result to the experiments [20]. In the last couple of years, the porosity in the matrix has been getting attention, regarding how it can be a negligible factor or a significant influence under external loads [21].
To improve the modelling of MMSFs, the application of the substructure technique is proposed. The basic idea of the substructure technique was conceived a long time ago for the reason that large models (e.g., airplanes) could be treated as separate units and then reunited along the joining perimeter [22]. This paper builds on this basic idea but applies it to a very different set of problems. This work is based on the creation of so-called super finite elements. In general, model reduction procedures based on the substructure technique are commonly used for studies involving model updating (e.g., modal analysis) [23,24]. In addition, if the program is parallelized, the runtime can be reduced even more drastically. Fortunately, this method is easy to parallelize [25]. In this way, the substructure technique can reduce computation time by up to three orders of magnitude [25,26]. In the modal analysis of large structures (such as bridges), it has been shown that if many identical substructures can be formed, then the calculations can be run faster. This finding will be an important reference later in this paper [27]. Regardless of whether the substructure technique is used for dynamic or quasi-static cases, the stiffness matrix is always decomposed in a similar way. The nodes are clustered in the connection points of the substructures and other points. From this separation, a reduced matrix of stiffnesses will be derived that combines the stiffness of all points into connection nodes. This will give significantly lower degrees of freedom (DoF) system [28,29]. Often, large structures are divided into completely periodic, repeating parts to form substructures. In some cases, 2D models are used to test the efficiency of the methods, as their effectiveness is already becoming apparent [30]. Several studies compare different numerical and analytical techniques, but it is clear from all of them that no universal formula exists, and that each method must be optimized for each problem [31,32].
The main goal of this work is to use the substructure technique to provide a reduced DoF model that considers the disorder of the interior structure of MMSFs while keeping the model treatable at runtime. The presented program is based on the displacement method within finite element calculations. To develop the methodology, a two-dimensional model was used to make the presentation of the subsequent steps even easier to follow. The method was mainly developed for small deformations and static loads.

Introduction of the substructure technique
Describing the real geometry of metal foams and solving its mechanical behavior is a challenging task even for a high-performance computer. To improve the effectiveness of the finite element method (FEM) various techniques can be applied such as reduced ordered model (ROM) or sub-structuring. The inhomogeneous interior structure of metal foams is usually treated as homogeneous with generalized material properties, but much interior, strength-relevant information is lost. The foundation of these techniques, assuming that dynamic effects are disregarded, is the discretized equation of equilibrium (Eq. (1)): where K is the stiffness matrix of the structure, u is the nodal displacement vector, and f is the load vector.
Here, only small displacements and stains are taken into account. An example of a plane strain problem is shown in Fig. 1, which in this case shows a rectangular geometry for simplicity. But it can be any other two-dimensional geometry. Fig. 1 shows a quadrilateral mesh, but these can be replaced by triangular elements too. Displacements are fixed in the x-y directions at the left side of the meshed region and the upper right side is loaded. The primary problem with this approach is that for models with a rising number of DoF the memory usage and computation time also increases significantly. The basic idea of the substructure technique is to decompose the mechanical model into parts and incorporate their interior stiffness values into to the boundary points. While preserving the original geometry this method results in a mechanical model that behaves under external loads in the same way as the originally meshed model, but with a significantly reduced DoF.
To demonstrate the sub-structuring consider a bit more complex model ( Fig. 2) with FE mesh than the previous one. The substructures are marked with a big yellow number which can be called super finite elements. The black dots represent the inner nodes and the green dots the edge nodes. The boundary nodes lie on both common and non-common edges of the substructures. Interior nodes are eliminated, and their effect is taken into account through a modified stiffness matrix.
The method is illustrated on a 2 × 2 substructure with mesh and then extended to more complex structures. The stiffness matrices of the adjacent substructures, grouped into the interior and external nodes, look like Eq. (2): where u b is the boundary and u i is the interior column vector of the displacements in the directions of x and y, K bb is the stiffness matrix concerning the boundary points, K ib and K bi describe the connection between the boundary and the interior points, and K ii is the stiffness matrix concerning the interior points. , , and are equal to zero, f b including the load given for all the perimeter points. Because the stiffness matrix K bb needs to be split into four separate parts, each substructure has its own boundary point, several of which coincide. Then, the load f b and displacement u b vectors must also be decomposed into as many pieces as many substructures exist. Each of them must be considered separately with its traction and displacement on the surface. The matrix equation is:  Performing the matrix multiplication and expressing the first and last two lines in terms of the displacements of the interior nodes results in the following algebraic system of Eqs. (4) to (7): After explaining the remaining third row, Eq. (8) Here it can be seen that Eqs. (4), (5), (6) and (7) result in the back-transforming matrices: where n = 1, 2, 3, 4. It will be indispensable later.
The new, generalized formula of Eq. (11) is the following for n substructures and m substructure loads: where n denotes the number of substructures and m the number of load vectors. Based on Eq. (11), it is easy to see that if the given part is divided into many substructures, the stiffness values of the same substructures do not need to be calculated from scratch but can be stored in memory and simply added to the stiffness values in the next step. The effectiveness of the substructure technique depends largely on the ratio of the edge points to the interior points. If this ratio is as small as possible, the running time can be significantly reduced.

Application of substructure technique to metal foams
Modelling conventional metal foams is a challenging task due to their internal random porosity. For MMSFs, this is even more complicated. The problem with CT-based mechanical models or other techniques describing the internal structure is that the created models have too high DoF, which is extremely time-consuming to solve. For this reason, the substructure technique was used, but in a way, it had never been applied before. The aim was not to develop a new homogenization technique, but to maintain the internal inhomogeneous structure. The first step was to create a few differentiating substructures. An important criterion is that the boundary points of all substructures (super finite elements) must be the same, as in the case of conventional finite elements. In the following, it will be weighty to take into account the fact that the more identical substructures are in the whole structure, the more efficient the memory usage will be. In addition, if the points on the edges are equivalent, then one super finite element can be replaced by another. Compared to the solid structure, the only modification required about the foams is to remove the corresponding finite elements from the model with their appropriate stiffness values. Exactly how many finite elements should be left out of a given substructure is determined by what percentage of the volume of the foam to be modelled as a cavity. The hollow parts are filled with gas in the real case, but the existence of gas is mechanically negligible.
In the case of MMSFs, this is complemented by the need to substitute a different material constant matrix around the cavity in a given layer when calculating the stiffness matrix. Fig. 3 shows a substructure where the filler hollow sphere is shown in green. 63% of the internal part of the substructure is removed, because this is the theoretical maximum gas percentage of metal MMSFs [7]. Prefabrication of substructures with different internal spherical geometries is the key to this method. The efficiency of the method is that the whole structure is built by producing and assembling the predefined finite number of pieces of different substructures (e.g., Fig. 3). The stiffness matrices of these small substructures can be easily Fig. 3 Representative MMSF unit for substructure technique generated. Thus, the external geometry and stiffness of the whole model are the same as that of the conventional model, but with a much smaller number of DoF. Finally, only summarizing stiffnesses the substructures is needed to produce the stiffness matrix of the reduced structure. The connectivity matrix of the super finite elements will provide information on the summing of stiffnesses. The geometry of the resulting model and the displacement values calculated on the external surface are completely consistent with the conventional method. The advantage will be the faster solvable algebraic equation system due to the smaller number of unknowns.
An example of this method can be seen in Fig. 4, where nine substructures with different interiors have been predefined. In such a case, only the stiffness of the nine differentiating super finite elements is calculated and separated into boundary and interior points. Then the matrix operations behind the first sum of Eq. (11) must be performed for the nine cases. Afterward, the matrix operations behind the first sum of Eq. (11) must be performed for the nine cases. An important criterion is that the substructures must be identical at the edges, as the nodes must overlap for continuity when joining. As an analogy, this is understood in the same way as in the rules of dominoes that only pairs of identical elements can be joined. A normal distribution is the best way to approach real inhomogeneity. Thus, it is only needed to insert the calculated stiffness values of the nine cases into the corresponding positions. In this way, it is possible to construct inhomogeneous random internal structures in a memory-efficient way using a few different sub-structures.
The irregularity of the internal structure of the model can be further increased by varying the diameters of the hollow spheres. This method is applied in the same way to the three-dimensional case, but with three-dimensional substructures.

Demonstrating the effectiveness of the substructure technique on different MMSF structures
Consider a practical implementation of the two-dimensional foam structure. The model built from substructures is subjected to loading. The 10 × 10 mm model substructure partitioning and specifying boundary conditions are shown in Fig. 5. Each substructure includes a mesh spacing of 0.0125 mm edge length. The internal structure of the super finite elements includes a 63% space-filling of the spherical cavities, but this is no longer visible on the reduced structure. Equation (11) was used to calculate the displacements of the boundary points in the full substructure model.
The displacements of the boundary points were interpolated to produce the reduced displacement field shown in Fig. 6. For better visibility of displacements and line widths the scale factor is ten.
The resulting plots will usually show the largest displacements, since the accurate displacements of the full contours of the models will always be known. However, this method does not yet provide information on the  internal points at this stage. Equation (9) gives the possibility to determine the complete displacement field of the selected substructure. The entire displacement field of such a super finite element is shown in Fig. 7.
The key to the effectiveness of the method lies in the minimized use of memory. Previously, it was mentioned that the geometry of the model is described by predefined and reduced-size substructures, but the memory efficiency of this was not discussed. These predefined super finite elements just need to be stored in the memory to build up the whole model with them. Relative efficiency increases with the larger model size and when the number of types of substructures is fixed. Consider a bar diagram (Fig. 8.) where a predefined substructure is used to describe models of increasing size. Two pre-defined substructures were used to determine the values of the diagram. The vertical axis shows the ratio of memory allocation of the stiffness matrix between the predefined substructures and the total reduced structure.
Equation (12) is used to calculate the memory utilization rate (MUR), where m denotes the pre-defined substructures and M denotes the total reduced memory allocation. It can be seen, that the preallocated memory requirement decreases exponentially as the number of required substructures increases (see Table 1 and Fig. 8).
The proportion of predefined substructures that allocates memory becomes smaller compared to the total structure. Since a real model will consist of many substructures, it is easy to see that the initial memory allocation will converge to zero.
Increasing the internal mesh density of the substructures will not affect the rate of run-off, as it will increase the size of the overall reduced structure. The only parameter affecting the ratio reduction is the number of predefined substructures. Increasing this number reduces the convergence efficiency. Since the vertical axis is a ratio of the memory allocations, increasing the number of pre-defined substructures would increase the values of the ratios in direct proportion as shown in Fig. 9. Each color represents an increasing number of predefined substructures.
One method to achieve more accurate results in FE calculations is to increase the density of the mesh. But this increases the size of the system of equations to solve.   Regardless of the way the system of equations is solved, this significantly increases the computation time.
If the number of substructures in a model is fixed, but the basic mesh is continuously densified, the effectiveness of the method becomes apparent which can be clearly seen in Table 2. This difference can be illustrated in Fig. 10.
In the case where the number of substructures is not fixed the interpretation of running times becomes more complex. In such cases, the shortest running time can be found by varying the reduction ratio of the boundary points (Eq. (12)) and the number of substructures.
If the calculation is performed on models with a fixed number of DoF at multiple points, the colored curves can be seen in Fig. 11 are obtained. The reason for including these two data on the horizontal axes is that their mutual interaction will clearly influence the running time. An important aspect of Fig. 11 is that since data on horizontal axes interact with each other, the value of one axis must always be fixed when searching for the optimum. The interaction of the horizontal axes can be easily seen from the fact that if the number of substructures is increased, the reduction ratio will obviously change and vice versa. Side views 1 and 2 in Fig. 11 show how the running time increases dramatically with the reduction ratio and slightly with the number of substructures. From the top view, it can be seen that the curves with a constant number of DoF take on a hyperbolic shape.
Each curve has a minimum point associated with the minimum running time. It can be found that the minimum  points of curves with increasing DoF belong always to higher reduction ratios (Eq. (12)). The reduction ratios for the minimum running time are summarized in Table 3. This optimal ratio shows a slowly increasing trend, as shown in Fig. 12. The implication for models with larger memory allocation sizes or denser meshes is that more substructures can be used to achieve minimum runtimes. However, this trend is not directly proportional to the number of DoF, but only slightly increasing.

Conclusions
The substructure technique has been presented for calculating small deformations of MMSFs. The conventional finite element computation method for the plane deformation problem has been compared with the substructure technique in 2D. The substructure technique can significantly reduce the running time, which is becoming more pronounced as the size of the model increases. This work has shown that homogenization techniques or CT-based meshing are not the only methods to describe metal foam structures. The substructure technique can combine the beneficial features of both methods and neglect their disadvantages. An important finding is that with pre-defined substructures it is possible to build the internal foam structure in a memory-efficient way without increasing the runtime as in CT-based methods. The other major result is that as the number of DoF is increased, the reduction ratio associated with the minimum running time shows a slight increase. This slightly increasing reduction trend shows the number of substructures that, depending on the size of the model, can achieve the minimum running time.  . 12 The trend of the minimum running time points with increasing DoF