Numerical Analysis, Optimization, and Multi-Criteria Design of Vacuum Insulated Glass Composite Panels

The subject of this study is Vacuum Insulated Glass (VIG) panels, which consist of two glass panes with an evacuated space and evenly distributed micro-support pillars between them. The deflection of panes towards the centre of the structure caused by atmospheric pressure is a mechanical problem that occurs in this type of structure. The aim of this study was to extend previous research on the optimal arrangement of support pillars in terms of eigenfrequencies and dynamics to include aesthetic aspects. Using Abaqus/CAE v2017 software, a large number of numerical models were created and subjected to a comprehensive multi-criteria analysis. Fractal analysis was employed to automatically assess the aesthetics of the proposed solutions. The study presents theoretical solutions that could be implemented in industrial production. The presented study shows that it is possible to effectively extend the criteria for optimizing the arrangement of pillars with new design criteria. Most studies focus on pillar placement, amount, or shape in terms of panes thermal or mechanical properties. Due to the increasing number of VIG panels applications in places exposed to external vibrations, other design criteria for VIG panels are also required and are provided by the following study.


VIG Panels
Nowadays, energy management is an important problem in the world economy. Its price has risen sharply in recent years. It is caused by the need to move away from fossil fuels and the dwindling resources of the planet. The recent geopolitical problems are an additional factor deepening the upward trend in energy prices. Therefore, any initiative aimed at saving energy is a priority for the economy.
The subject of this work is VIG (Vacuum Insulating Glass) panels. They represent the latest technological achievement that is already implemented in many countries. The VIG technology consists of two high-strength glass panes with a vacuum between them ( Figure 1). In order to counteract the atmospheric pressure and the contact of the panels, they are connected to micro-support pillars-elements usually made of steel with a diameter of 0.3-0.5 mm [1]. Moreover, the panes are sealed around the edges. Consequently, this design provides excellent thermal properties to VIG panels. Their U-value can be as low as 0.2-0.4 W/m 2 K [2].

Fractals in Engineering
In recent years, scientists and engineers have been looking for inspiration to create new, original structures with better properties. Among the many interesting trends, the attempts to apply structures called fractals in various fields of technology are remarkable [3][4][5]. These are mathematical structures discovered and described by Mandelbrot in the 1970s. Their representation in nature is significant [6]. They are characterized by a sophisticated, self-similar structure, the similarity scale of which is described by a parameter called fractal dimension. Over the next years, fractals will be subject of the extensive research. The individual section of mathematics, such as fractal geometry, was created. The wide possibilities of fractal applications were also analysed, both to create new technical structures and to use the so-called fractal analysis to describe the features of existing structures [7]. Fractals application in architecture and civil engineering [8] deserves special attention-particularly in the dynamics of plate structures [9]. Fractals have been seen in architecture across many eras and cultures.
The richness and variety of fractal forms naturally occurs in nature but also in the architecture of many cultures. This leads to natural questions of the aesthetic value of fractals. The proportions related to the similarity of fractal structures provides the observer an aesthetic feeling which has already been the subject of scientific research [10]. The mathematical description of fractal structures gives the engineers possibility to quantify the aesthetic criterion that can be used for generative design of building objects [11].

Motivation
Thermal properties of VIG panels [12,13], strength properties of support-pillars and their static properties [14,15] were subject of numerous studies. However, in recent years VIG technology has been analysed also in terms of its acoustic properties [16].
An interesting trend in the recent research on VIG technology are studies on its dynamic properties. In [17], the authors, using an advanced FEM model, analysed a number of plates with different geometries and compared their dynamic properties. By determining the value of the natural frequency of these structures, it is possible to determine the possibility of using VIG plates in buildings exposed to external dynamic factors, such as traffic or industrial vibrations. Another work worth mentioning is [18], in which the authors made an attempt to optimize pillar distribution in terms of the dynamic properties of the entire structure. They used an analytical model based on the analysis of classical equations of motion of thin plates elastically connected with connectors. Then, using genetic algorithms and applying dynamic optimization criteria, they obtained non-trivial topological solutions for the structure of VIG plates.
The motivation to undertake this work was the willingness to summarize and also continue the research contained in [18]. It is worth emphasizing that authors use a com-pletely different theoretical model of VIG plates. Combining, comparing and analysing advanced FEM models with a simple analytical approach allows the selection of models of this type of structure that is easier to apply in the future.
Another motivation to start this work were the conclusions from the work [18]. The authors, after analysing the dynamic properties of many VIG-type structures, concluded that the optimization used criteria require further refinement and extension, as they slightly affect the distribution of pillars. Hence, the need to extend the optimization criterion with new aspects, which is undoubtedly the fractal analysis. By determining the correlation coefficient of a given structure with the mathematical structure of the fractal, it is possible to evaluate the aesthetic properties of the arrangement of pillars in the VIG panels. This would enable a more comprehensive approach to the optimization of these types of structures.

Purpose of the Work
The first aim of this work is to compare the accuracy of various numerical models of VIG plates in terms of dynamics. Demonstrated models are simple and easier to apply. They provide correct and accurate results in relation to advanced three-dimensional structure models. They allow for faster computations, which is important in the case of genetic calculations and optimization.
The second important goal of the work is the use of multi-criteria optimization to design the arrangement of connectors in VIG panels. By extending the objective function with a parameter that determines the degree of similarity of a given structure to a fractal, a simple quantitative parameter characterizing the aesthetic properties of panels was obtained. The relationship between the aesthetic properties of the given structure and its fractal properties has been investigated by other authors for various problems [10]. The distribution of pillars, which in the case of structural elements intended for use in exposed places such as window partitions, is of key importance for the comfort of users. To determine a parameter characterizing the degree of similarity in the arrangement of pillars on the VIG plate, and thus to investigate its aesthetic properties, the mass-radius method was used [19].

Geometry and Support Conditions
The main issue in selecting the geometry and support conditions was the need to enable compatibility between all analysed modelling methods, which differ significantly in their approach. Only by providing the same support conditions obtained results were comparable to each other. Ultimately, simple support conditions were assumed on all edges of the plate (i.e., all displacements were blocked and all rotations were released at all nodes- Figure 2). In the case of the three-dimensional FEM model, the displacements were blocked only in the internal panes of both plates (Figure 3). In the case of blocking the displacements along the entire thickness of the glass elements, restraint would have appeared on a given edge [17].  In the case of the three-dimensional FEM model, the displacements were blocked only in the internal panes of both plates ( Figure 3). In the case of blocking the displacements along the entire thickness of the glass elements, restraint would have appeared on a given edge [17].
In the case of the three-dimensional FEM model, the displacements were blocked only in the internal panes of both plates ( Figure 3). In the case of blocking the displacements along the entire thickness of the glass elements, restraint would have appeared on a given edge [17]. Support of the middle panes of both glass plates was not considered as this would significantly increase the stiffness of the entire model. The thickness of the glass plates is much greater (6-8 mm) than the vacuum thickness (0.3 mm). This makes the arm of horizontal support forces much (about 20 times) greater than the arm of corresponding forces in the simplified model.

Simple FEM Model
A simple FEM model was made using the commercial engineering software Autodesk Robot Structural Analysis Professional v2022. A shell-member model was built, in which the glass plates were represented as flat shell elements simply supported in their middle panes. Pillars, represented as bar elements with six degrees of freedom in each node, were placed between plates ( Figure 4). The division of the regular FEM mesh of both plates was adjusted so the mesh stands out at the point of the pillar's connection. Support of the middle panes of both glass plates was not considered as this would significantly increase the stiffness of the entire model. The thickness of the glass plates is much greater (6-8 mm) than the vacuum thickness (0.3 mm). This makes the arm of horizontal support forces much (about 20 times) greater than the arm of corresponding forces in the simplified model.

Simple FEM Model
A simple FEM model was made using the commercial engineering software Autodesk Robot Structural Analysis Professional v2022. A shell-member model was built, in which the glass plates were represented as flat shell elements simply supported in their middle panes. Pillars, represented as bar elements with six degrees of freedom in each node, were placed between plates ( Figure 4). The division of the regular FEM mesh of both plates was adjusted so the mesh stands out at the point of the pillar's connection.

3D FEM Model
The model implemented in ABAQUS based on three-dimensional C3D8R elements: first order, reduced integration element with hourglass control activated was selected as the most advanced 3D model using FEM [17] (Figure 5).

3D FEM Model
The model implemented in ABAQUS based on three-dimensional C3D8R elements: first order, reduced integration element with hourglass control activated was selected as the most advanced 3D model using FEM [17]

3D FEM Model
The model implemented in ABAQUS based on three-dimensional C3D8R element first order, reduced integration element with hourglass control activated was selected a the most advanced 3D model using FEM [17] (Figure 5).  (3 where is volume of the element and I = 1, 2, 3. In the centroidal strain formulation th gradient matrix is simply given: Considering the above, it may be seen that centroidal strain formulation reduces th amount of effort required to compute the gradient matrix. In ABAQUS the artificia stiffness method and the artificial damping method showed in [21] are used to control th The interpolation function for C3D8R element shown in Figure 5 is given by the equation: where I denote the node of the element. The last four vectors Γ I α , α = 1, 2, 3, 4 are the hourglass base vectors. The gradient matrix B I is defined by integrating over the element: The shape functions are the same as for the C3D8 element and can be found in [20] (2): where V el is volume of the element and I = 1, 2, 3. In the centroidal strain formulation the gradient matrix is simply given: Considering the above, it may be seen that centroidal strain formulation reduces the amount of effort required to compute the gradient matrix. In ABAQUS the artificial stiffness method and the artificial damping method showed in [21] are used to control the hourglass modes in these elements. In the work [22], the effectiveness of C3D8R elements with hourglass control was compared to other types of elements.

Analytical Model
The analytical model included in this study was based on the assumptions made in [18]. Let Ox 1 x 2 x 3 be an orthogonal Cartesian coordinate system. Two plates of thicknesses h i , i = 1, 2 made of the same material (density ρ, Young's modulus E and Poisson number ν) were considered. There are n springs, each of stiffness k j , located in points with coordinates x 1,j , x 2,j , j = 1, . . . , n are located between those plates ( Figure 6). hourglass modes in these elements. In the work [22], the effectiveness of C3D8R elements with hourglass control was compared to other types of elements.

Analytical Model
The analytical model included in this study was based on the assumptions made in [18]. Let Ox1x2x3 be an orthogonal Cartesian coordinate system. Two plates of thicknesses hi, i = 1, 2 made of the same material (density ρ, Young's modulus E and Poisson number ν) were considered. There are n springs, each of stiffness kj, located in points with coordinates ( 1, , 2, ), = 1, … , are located between those plates ( Figure 6). Due to the small thickness of the glass plates and their slight deformations, the classical theory of thin plates, known as the Kirchhoff-Love theory, was applied. According to the studies conducted by the authors in [15], the results obtained using the linear and classical theory for frequency values are fully realistic. According to the theory of thin Due to the small thickness of the glass plates and their slight deformations, the classical theory of thin plates, known as the Kirchhoff-Love theory, was applied. According to the studies conducted by the authors in [15], the results obtained using the linear and classical theory for frequency values are fully realistic. According to the theory of thin plates, the displacement takes the form of: strain-displacement relations are presented as: The action functional of the problem is given by: where the strain energy is equal to: and kinetic energy is given by: Let w i = w i (x 1 , x 2 , t) be the transverse deflection; overdot stands for the derivative with respect to time t = t 0 , t 1 . The coefficients D i αβγδ represent the plate stiffness tensor elements and are given by the following formula: where the plate flexural rigidity is equal to: The derivatives are given as follows: Coefficient ∂ αβ stands for Kronecker delta. The indices α, β, γ, ω are related to the Ox 1 x 2 coordinate system and run through 1, 2.

Natural and Forced Vibrations
In this study, the dynamic properties of VIG plates were investigated in terms of natural vibrations in which various approaches to numerical modelling were compared and in terms of forced vibrations. The eigenproblem can be expressed in a standard form: where K, M and u stand for matrix of stiffness, matrix of mass and vector of displacement, respectively. The displacement vector is equal to: Elements of matrices K and M are, respectively, given by following formulas: Considering forced damped vibrations the system of equations is assumed as: where K, C, M stand for stiffness, damping and mass matrices, respectively, Q is load, u is displacement vector Ω is forcing frequency and t is time. The damping matrix has been defined so that the damping could be considered as slight: where α is assumed equal to 10 −2 [18]. The formula describing a uniformly distributed load is assumed as: where q j are known coefficients dependent on the loading case. Application of the multimodal approach leads to the following solution: Coefficients a C and a S are defined as: 2.6. Simple Genetic Algorithm (SGA) SGA (simple genetic algorithm) was selected as the tool for topology optimization of the plate. A single individual (VIG plate) was encoded by genotype in the form of an array with values 1 or 0 representing the existence or absence of a pilaster at a given position in the plate [18]. An example of coded pillars arrangement is shown in Figure 7.

Simple Genetic Algorithm (SGA)
SGA (simple genetic algorithm) was selected as the tool for topology optimization of the plate. A single individual (VIG plate) was encoded by genotype in the form of an array with values 1 or 0 representing the existence or absence of a pilaster at a given position in the plate [18]. An example of coded pillars arrangement is shown in Figure 7. Other SGA parameters, such as the selection of an appropriate method of interbreeding of individuals taking into account the characteristics of the problem, are described in detail in [9].

Optimization Criteria and Fractal Analysis
In this paper, several optimization criteria of various types have been defined. The first criterion used was the dynamic one, defined analogously to the work [18], i.e., as the average of the ratio of the maximum deflection of the lower and upper plate in the frequency range from 220 Hz to 440 Hz: The function (23) for a solid plate (i.e., a plate with evenly spaced pillars) is not less than 0.985. However, as shown in [18], limiting ourselves to this criterion gives only a small possibility of topological optimization of the joint arrangement.
Fractal analysis was used as another criterion for advanced optimization of pillar distribution. It determines the degree of similarity of a given structure to the theoretical fractal structure. There are many methods of determining the so-called fractal dimension of the structure. Due to the point nature of the pillars distribution decision to use an algorithm for the objective function calculation based on the simplified mass-radius method described in [19] was made. This algorithm was implemented in two steps: 1. The so-called spinning radius R (Figure 8) varying from 0 to the maximum value depending on the size of the given plate, and for each R, the number of pilasters in total M(R) was counted in circles with centres in all pillars of a given spinning radius R; 2. Because for an ideal fractal structure there is a relationship [19]: Other SGA parameters, such as the selection of an appropriate method of interbreeding of individuals taking into account the characteristics of the problem, are described in detail in [9].

Optimization Criteria and Fractal Analysis
In this paper, several optimization criteria of various types have been defined. The first criterion used was the dynamic one, defined analogously to the work [18], i.e., as the average of the ratio of the maximum deflection of the lower and upper plate in the frequency range from 220 Hz to 440 Hz: The function (23) for a solid plate (i.e., a plate with evenly spaced pillars) is not less than 0.985. However, as shown in [18], limiting ourselves to this criterion gives only a small possibility of topological optimization of the joint arrangement.
Fractal analysis was used as another criterion for advanced optimization of pillar distribution. It determines the degree of similarity of a given structure to the theoretical fractal structure. There are many methods of determining the so-called fractal dimension of the structure. Due to the point nature of the pillars distribution decision to use an algorithm for the objective function calculation based on the simplified mass-radius method described in [19] was made. This algorithm was implemented in two steps:

1.
The so-called spinning radius R (Figure 8) varying from 0 to the maximum value depending on the size of the given plate, and for each R, the number of pilasters in total M(R) was counted in circles with centres in all pillars of a given spinning radius R; 2.
Because for an ideal fractal structure there is a relationship [19]: where D is the so-called fractal dimension. In order to investigate the degree of similarity between VIG plate structure and the ideal fractal a logarithmic approximation of M(R) data series was used. Then, r coefficient was determined which defines the degree of adaptation the logarithmic curve to the data described in the previous step. Classical formulas for the least squares approximation were used: where x i = log(R) and y i = log[M(R)]. The correlation coefficient r obtained in (25) was used directly as an objective function. Its values are in the range from 0 to 1. Equal weights and linear combination of the above objective functions were used for the multi-criteria analysis: The results obtained solely on the basis of the criterion defined by the Formula (21) were shown and analysed in [18]. This paper focuses on the results obtained for the criteria defined by the Formulas (23) and (24).

where
= log( ) and = log [ ( )]. The correlation coefficient r obtained in (25) was used directly as an objec tion. Its values are in the range from 0 to 1.
Equal weights and linear combination of the above objective functions wer the multi-criteria analysis: The results obtained solely on the basis of the criterion defined by the Fo were shown and analysed in [18]. This paper focuses on the results obtained f teria defined by the Formulas (23) and (24).

Comparison of Various Numerical Models
Comparative numerical calculations were performed for two pillars system cases, a plate with a glass pane thickness of h = 7.5 mm, vacuum width hv = 0. plate dimensions 2.0 m × 1.0 m was used. Glass with Young's modulus E Poisson's ratio ν = 0.33 and density = 2500 kg/m 3 was assumed. Support-pill mm height and diameter. The material for pillars is steel, with Young modulu GPa, the mass density ρs = 7850 kg/m 3 , and the Poisson's ratio νs = 0.3. The foll lars systems were used ( Figure 9):  For the analysis using FEM, in addition to the given material data, the following approaches were adopted:

Comparison of Various Numerical Models
I.
Simple FEM model-flat plate quadrilateral finite elements with a relatively sparse division of 50 mm. II.
The 3D FEM model-three-dimensional finite elements were used. The glass panel were divided with a mesh size of 2.5 mm (0.5 mm near the support pillars), while the support pillars themselves were divided more accurately. For the analysis using FEM, in addition to the given material data, the following approaches were adopted: I. Simple FEM model-flat plate quadrilateral finite elements with a relatively sparse division of 50 mm. II. The 3D FEM model-three-dimensional finite elements were used. The glass panels were divided with a mesh size of 2.5 mm (0.5 mm near the support pillars), while the support pillars themselves were divided more accurately.
In the case 0375, the following results of the first several natural frequencies were obtained (Table 1). The 0375X was analysed as the second case, which differed in more pillars in the centre of the slab. In the case of the 0375X, the following results of the first sixteen natural frequencies were obtained (Table 2). Table 2. The first sixteen eigenfrequencies calculated in three different numerical models for the case number 0375X. The frequencies that are most influenced by the problem of mapping with the method of supporting the 3D model are marked in green. Vibrations in counter-phase are marked in blue. When comparing the values of individual frequencies of natural vibrations, the modes of natural vibrations obtained in all models were analysed. They can be divided into two groups: vibrations in phase and vibrations in counter-phase. In the former case, both plates swing in the same direction at the same time (Figure 10a-d), in the latter case, plates swing in the opposite direction (Figure 10e-h). When comparing the values of individual frequencies of natural vibrations, the modes of natural vibrations obtained in all models were analysed. They can be divided into two groups: vibrations in phase and vibrations in counter-phase. In the former case, both plates swing in the same direction at the same time (Figure 10a-d), in the latter case, plates swing in the opposite direction (Figure 10e-h). Figure 10. Selected z-normalized wane vibration modes for plate 0375: (a-d) vibration in phase, (e-h) vibration in counter phase; (a) 1st mode shape, (b) 2nd mode shape, (c) 4th mode shape, (d) 11th mode shape, (e) 3rd mode shape, (f) 6th mode shape, (g) 10th mode shape, (h) 14th mode shape.

Shape Number
It is worth noting that in the case of vibrations in the phase, the pilasters have no significant effect on the shape and the result of the natural vibrations values because both panes vibrate identically-the pillars do not undergo major deformations. However, in the case of vibrations in the counter-phase, the pillars undergo significant deformations, and their parameters and modelling method are of key importance for the final result. In Figure 10. Selected z-normalized wane vibration modes for plate 0375: (a-d) vibration in phase, (e-h) vibration in counter phase; (a) 1st mode shape, (b) 2nd mode shape, (c) 4th mode shape, (d) 11th mode shape, (e) 3rd mode shape, (f) 6th mode shape, (g) 10th mode shape, (h) 14th mode shape.
It is worth noting that in the case of vibrations in the phase, the pilasters have no significant effect on the shape and the result of the natural vibrations values because both panes vibrate identically-the pillars do not undergo major deformations. However, in the case of vibrations in the counter-phase, the pillars undergo significant deformations, and their parameters and modelling method are of key importance for the final result. In general, the obtained results are in line with expectations. In most cases, the results obtained from all three models are consistent with each other, with differences not exceeding 5-6%.
However, in both cases, both 0375 and 0375X, the results obtained for models I and II (analytical model and FEM model) in the case of the lowest vibration frequencies significantly differ from the values obtained for model III (FEM 3D)-even by 24% (marked green in Tables 1 and 2). The results for the 3D model are significantly higher. This is due to a completely different method of supporting these models, which is of greatest importance precisely for the first forms and frequency of vibrations. As the thicknesses of the glass plates is much greater than the thickness of the vacuum, a realistic 3D model is much stiffer than the corresponding 2D model. It was observed after analysing the fragment of the first mode of natural vibrations located in the support zone ( Figure 11).
The form of deformation, in the case of 3D FEM model (Figure 10a), resembles a fixed scheme more than a simply supported scheme (Figure 10b).
For some cases of counter-phase vibrations large discrepancies in the results between the models were obtained, i.e., 3 ( Figure 10e) and 10 ( Figure 10g). Method of modelling and meshing of pillars may have had an influence on them. This problem was analysed in [17].
The results from the 0375 and 0375X models were also compared. It can be clearly seen that, in the case of the 0375X, the use of pillars in the centre of the slab completely eliminated the 3rd form of vibrations in the 0375 model (Figure 10e, Tables 1 and 2). green in Tables 1 and 2). The results for the 3D model are significantly higher. to a completely different method of supporting these models, which is of g portance precisely for the first forms and frequency of vibrations. As the thic the glass plates is much greater than the thickness of the vacuum, a realistic 3 much stiffer than the corresponding 2D model. It was observed after ana fragment of the first mode of natural vibrations located in the support zone (F The form of deformation, in the case of 3D FEM model (Figure 10a), r fixed scheme more than a simply supported scheme (Figure 10b).
For some cases of counter-phase vibrations large discrepancies in the tween the models were obtained, i.e., 3 ( Figure 10e) and 10 ( Figure 10g). modelling and meshing of pillars may have had an influence on them. This pr analysed in [17].
The results from the 0375 and 0375X models were also compared. It can seen that, in the case of the 0375X, the use of pillars in the centre of the slab eliminated the 3rd form of vibrations in the 0375 model (Figure 10e, Tables 1 a

Multi-Criteria Optimization
The rest of the numerical calculations for a number of cases which var mensions of the VIG plate were performed (length L1, width L2 and the thick glass plate h). Other geometrical and material parameters were assumed as vious analysis. The list of case numbers along with the corresponding comb dimensions are summarized in Table 3.

Multi-Criteria Optimization
The rest of the numerical calculations for a number of cases which varied the dimensions of the VIG plate were performed (length L 1 , width L 2 and the thickness of the glass plate h). Other geometrical and material parameters were assumed as in the previous analysis. The list of case numbers along with the corresponding combinations of dimensions are summarized in Table 3. The analytical model of VIG plate and optimization using previously described genetic algorithms were used for the calculations. The following pillars distributions inside the VIG plate were obtained ( Figure 12).
For square plates, the optimal forms are point-symmetrical with rotational symmetry of 90 • angle (cases 0125 and 0175). As the disproportion between plate length and width increases, the optimal forms develop toward axisymmetric systems. An interesting relationship for the cases with a thicker glass plate (0175 and 0275) may be noticed. Due to the greater stiffness of the individual glass panes and the consequent different degree of energy transfer between the plates correlated with the objective function (21), the optimal forms aim at the ones shown in Figure 13.

0350
2.00 1.00 5.00 0375 2.00 1.00 7.50 The analytical model of VIG plate and optimization using previously described genetic algorithms were used for the calculations. The following pillars distributions inside the VIG plate were obtained ( Figure 12). For square plates, the optimal forms are point-symmetrical with rotational symmetry of 90° angle (cases 0125 and 0175). As the disproportion between plate length and width increases, the optimal forms develop toward axisymmetric systems. An interesting relationship for the cases with a thicker glass plate (0175 and 0275) may be noticed. Due to the greater stiffness of the individual glass panes and the consequent different degree of energy transfer between the plates correlated with the objective function (21), the optimal forms aim at the ones shown in Figure 13. In the case of using the fractal objective function defined by the Formula (23), the optimal solutions become more complicated. It is clear that trivial solutions are rejected. In the case of a square plate (0125), it is characterized by quasi-rotational symmetry with 90° angle. However, in the case of rectangular plates (0225 and 0325) the evolution of forms in relation to those shown in Figure 13 towards more differentiated solutions, but with the preservation of some areas of a linear arrangement of pillars.  For square plates, the optimal forms are point-symmetrical with rotational symmetry of 90° angle (cases 0125 and 0175). As the disproportion between plate length and width increases, the optimal forms develop toward axisymmetric systems. An interesting relationship for the cases with a thicker glass plate (0175 and 0275) may be noticed. Due to the greater stiffness of the individual glass panes and the consequent different degree of energy transfer between the plates correlated with the objective function (21), the optimal forms aim at the ones shown in Figure 13. In the case of using the fractal objective function defined by the Formula (23), the optimal solutions become more complicated. It is clear that trivial solutions are rejected. In the case of a square plate (0125), it is characterized by quasi-rotational symmetry with 90° angle. However, in the case of rectangular plates (0225 and 0325) the evolution of forms in relation to those shown in Figure 13 towards more differentiated solutions, but with the preservation of some areas of a linear arrangement of pillars. In the case of using the fractal objective function defined by the Formula (23), the optimal solutions become more complicated. It is clear that trivial solutions are rejected. In the case of a square plate (0125), it is characterized by quasi-rotational symmetry with 90 • angle. However, in the case of rectangular plates (0225 and 0325) the evolution of forms in relation to those shown in Figure 13 towards more differentiated solutions, but with the preservation of some areas of a linear arrangement of pillars.

Summary and Conclusions
As a result of the numerical calculations and analyses the following conclusions were formed: • Natural frequencies analysis of VIG panels requires 3D model application, the results obtained from simplified models are significantly underestimated; • VIG panels are characterized by two types of vibration: in phase (both glass panes bend in the same direction) and in counter-phase (both glass panes bend in opposite directions), pillars, their geometry and modelling method are of key importance for vibrations in counter-phase; • It is possible to effectively extend the criteria for optimizing the arrangement of connectors with new design criteria; • Fractal analysis can be a tool for VIG panels design.
In further research, the authors intend to focus on refining the programming genetic algorithm and determining the fractal size of the structure, so that it is possible to significantly tighten the pilaster mesh, as well perform a detailed analysis of pillars parameters on the values and forms of natural vibrations.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.