Dynamic Scaling of a Wing Structure Model Using Topology Optimization

Abstract

In this paper, a dynamic scaling methodology is introduced to devise reduced scaled models of aircraft with the objectives of minimizing the development cost and exploring the design space. A promising way to accomplish this is using Topology Optimization (TO) for Additive Manufacturing (AM). Here, TO is employed to design a reduce scale model by matching its natural frequencies and mode shapes to those of a full scale model. Different TO strategies based on density approach are tested with the goal of achieving a dynamically scaled structure that can be manufactured. To achieve this goal, the TO solution should be free from intermediate densities, which is observed for some TO strategies but not all. When no penalization factor is applied: (i) the relative difference between natural frequencies is less than 1% and (ii) the estimated Modal Assurance Criteria (MAC) metric to evaluate the correlation between mode shapes is close to the ideal identity matrix. These results demonstrate the effectiveness of the dynamic scaling methodology. However, when using a penalization factor to avoid intermediate densities, the dynamic behavior correlation between full and scaled models degrades. This trend is more visible in the MAC metric, where off-diagonal terms above 20% and diagonal terms below 90% appear.

1. Introduction

Designing any complex new product requires several working hours, including expensive experimental tests, that usually implies a high Research and Development (R&D) cost. A good example is an aircraft design project, where the development cost represents a significant part of the final product cost, which will depend on selling several units to dilute its price on each sale. Furthermore, the aviation industry is studying more environmentally sustainable aircraft to reduce its carbon footprint. Some of the proposed solutions are unconventional such as the Blended Wing Body (BWB) [1] or the Truss-Braced Wing (TBW) [2]. These designs are likely to have a different dynamic behavior than the conventional tube and wing configuration. Therefore, it is important to study and develop them in an integrated way using not only Multidisciplinary Design Optimization (MDO) strategies [3], but also experimental test beds such as wind tunnel [4] or flight test [5] models. In this paper, the focus is on devising scaled test models that are able to mick the dynamic behavior of its full size counterpart. Producing meaningful scaled models to test unconventional features at early design stages might lead to the early detection of problems and avoid expensive redesigns after performing tests with a full-size model.

A dynamic scaling methodology is proposed to design meaningful reduced scale models to be experimentally tested. This methodology aims at designing a scaled model by means of Topology Optimization (TO) for Additive Manufacturing (AM) that presents the same dynamic behavior of the full size model. To obtain dynamic similitude [6] the mode shapes and associated natural frequencies should be adequately scaled down and matched. Dynamic scaling methodologies have been applied to devise both wind tunnel [7,8] and flight [9,10,11,12] test models for evaluating aeroelastic behavior of unconventional designs. As seen in previous publications [13,14], an adequate dynamic scaling strategy when combined with a suitable selection of aerodynamic parameters enables aeroelastic similarity for some flow conditions, even though Mach and Reynolds number similarities are difficult to ensure [15,16]. Thus, achieving dynamic similitude is a key step towards aeroelastic similitude. Most of these methodologies [9,11,13,14] recur to optimization routines for replicating the dynamic behavior of the full size model in the scaled one by means of sizing a prescribed structural layout [17] with an exact scaled down external geometry. Alternative strategies such as changing the wing planform [12] by means of MDO or using TO [18] has been recently proposed.

TO is seen as a promising numerical method to be included in a dynamic scaling methodology given the higher design freedom and the possible synergies with AM. Topologically optimizing a structure [19] has proved to be an efficient tool in minimizing weight for several engineering fields, including aeronautics [20]. Associated to the advantage of a higher design freedom is the problem dimension, which is also higher and often leads to complex designs to manufacture. However, its interest has exponentially grown in the past years with the recent progress in AM [21], including aeronautical applications such as Unmanned Aerial Systems (UAS) [22,23]. Even though substantial progress have been made to allow for manufacturing complex topologically optimized parts; issues such as overhanging structures, poor surface quality, porous infill layouts and intermediate densities [21,24,25,26,27] should be addressed when designing a part. These issues if handled in a post-processing phase for a dynamic scaling problem, are likely to modify both mode shapes and associated frequencies.

In this work, the main objective is to propose a dynamic scaling methodology that recurs to TO, instead of the commonly used sizing optimization, to design scaled models that represent the dynamic behavior of a full size model. To the authors best knowledge there are very few papers in the literature employing TO in the process of dynamically scaling structures which might lead to synergies with additive manufacturing. The effect of a penalization factor to avoid intermediate densities on the ability to replicate this dynamic behavior is also assessed.

The remaining of this paper is organized as follows. In Section 2, the proposed dynamic scaling methodology is introduced. Then its application to several test cases is discussed in Section 3 and the concluding remarks are presented in Section 4.

2. Dynamic Scaling Methodology

A dynamic scaling methodology comprising 3 main steps is proposed. First the full scaled model is defined and its modal parameters are determined (step-1). Based on this model and a predefined set of scaling factors, the scaled model is first built (step-2) and its structure is then topologically optimized to ensure that dynamic similarity between models (step-it). The flowchart of this methodology is depicted in Figure 1, where the reader should follow a counter-clockwise direction from the external green lines (step-1) passing through the blue lines (step-2) and converging to the red lines in the yellow box, representing the iterative topology optimization process (step-it). A detailed explanation of the proposed methodology applied to a wing model is given in the subsequent sub-sections.

2.1. Full Scale Model

For the full scale model the wing geometry proposed by [28] is scaled up to the dimensions depicted in Figure 2. Regarding the structure two parts are defined: an internal solid part (Part-1 in Figure 2b); and an external part modeled as a shell (Part-2 in Figure 2c). The external shell represents the wing skin and will be kept unaltered when scaled down to maintain the aerodynamic shape for a future aeroelastic scaling. The material properties of the skin are set to be “weaker” than those of the internal solid part. Two “fictitious” sub-parts with different material properties (Material-A and Material-B), as depicted in Figure 2b, are defined for the internal structure (this is a purely academic problem since no boundary condition between the different materials is applied.) This is conducted to enable the easy adjust of the contrast between different sub-parts to control the problem’s difficulty in reaching dynamic similarity. For instance, if two considerably different metals (e.g., Al and Ti) are chosen for the two sub-parts, it will be more difficult to reach the desired dynamic similarity than if the same material is used for both sub-parts.

A Finite Element Method (FEM) model of the above described wing is built in Abaqus^®. The skin and solid parts are modeled using Shell S4R and Solid C3D8R, respectively. Regarding the boundary conditions, these two parts are tied between each other, as illustrated in Figure 3a, and all degrees of freedom (DoFs) at the wing root are clamped, as depicted in he Figure 3b. To compute the natural frequencies (eigenvalues) and corresponding mode shapes (eigenvectors) necessary for the dynamic scaling methodology a modal analysis problem is solved. The mode shapes are extracted from several points distributed on the wing (marked in Figure 3c), such that those points are able to reproduce any mode shape in the range of interest. Here, only the first five mode shapes are considered for the cases analyzed.

A mesh size of 58.9k elements was defined for the FEM model as a compromise between computational cost and accuracy for the modal analysis problem, after performing a mesh convergence study. With this mesh the simulation takes around 4 minutes (Workstation with 24 cores Intel(R) Xeon(R) E5-2620 v2 @2.10 GHz totaling @50.4 GHz with 32 GB of RAM (Random Access Memory).) and the errors relative to the finest mesh considered are lower than 2%, as one can see from Figure 4.

After a mesh is chosen, this process is carried out automatically in the dynamic scaling methodology (Figure 1) be means of a python script (inpFile1.py) that is built in Matlab^® and runs Abaqus^®. The resulting modal information is written in a file (outFile1.txt) that will be read to establish the target parameters for the scaled model.

2.2. Scaled Model

There are several strategies that can be applied to derive the required set of scaling factors ($\lambda$) to define a dynamically scaled model. To form this set, the fundamental quantities of mass, length and time should be included in the three primary quantities chosen to derive the scaled model. One possibility is to define length (l), density ($\rho$) and pressure (p) in the scaled model as follows:

$${\lambda }_l=\frac{l_s}{l_f}{\lambda }_{\rho }=\frac{{\rho }_s}{{\rho }_f}{\lambda }_p=\frac{p_s}{p_f},$$

(1)

where the subscripts s and f denote scaled and full size models, respectively. The natural frequency scaling factor can be derived from the above factors, resulting in:

$${\lambda }_{\omega }=\sqrt{\frac{{\lambda }_p}{{\lambda }_{\rho }{\lambda }_l^2}}.$$

(2)

Here, for simplicity the density and pressure scaling factors are unitary and the scaled model is 10 times smaller than the full scale in size (i.e., ${\lambda }_l=0.1$). With these assumptions, the scaled model frequencies are targeted to be 10 times higher than those of the full scale (i.e., ${\lambda }_{\omega }=10$). To ensure an exact dynamic similarity also the damping factors and mode shapes should be matched. Since structural damping is difficult to model numerically it is often neglected. The mode shapes of full size and scaled models can be matched using the Modal Assurance Criteria (MAC) if a gradient free optimizer (e.g., a genetic algorithm) is used, which will be time consuming given the number of design variables involved in a topology optimization problem. For these reasons, neither damping nor mode shapes are directly matched in the dynamic scaling methodology here proposed. Nevertheless, the MAC is used to track the correct mode shapes for matching the natural frequencies, in the proposed methodology. MAC is a statistical indicator that is widely used to provide a quantitative comparison between two mode shapes [29] and can be expressed as follows:

$$\mathrm{MAC}(A,B)=\frac{{\left|{\Psi }_A^T\times {\Psi }_B\right|}^2}{\left({\Psi }_A^T\times {\Psi }_A\right)\times \left({\Psi }_B^T\times {\Psi }_B\right)},$$

(3)

where ${\Psi }_A$ and ${\Psi }_B$ denote the eigenvector A and B, respectively. A $n\times n$ square matrix can be constructed by comparing two sets of n eigenvectors, which ideally should be an identity matrix to ensure a perfect correlation between the two sets of mode shapes, for instance of the full-size and scaled models.

Based on the previous considerations, a scaled model 10 times smaller in size than the full scale one is built in Abaqus^®. Despite the geometry is merely a scaled down of the previous model, now only one material is used and the internal geometry needs to be prepared for optimization. To this end, a proper grid is needed to define the design variables using datum planes. For practical reasons, the internal mesh is much coarser than the one used for the full scale model since each of the internal mesh elements will be set as a design variable. Two different grids, depicted in Figure 5, are considered, (18, 24, 1) and (54, 72, 1), where (X, Y, Z) indicate the number of grid cells in x, y and z planes, respectively. It is worth to note that only the grid cells belonging to the wing geometry are set as design variables. When designing this grid, attention should be taken in regard to: (i) the points chosen to extract the mode shape information for comparison with the full scale model; and (ii) the mesh quality in the most incline boundary.

Again, python scripts, inpFile2.py and inpF_it.py, are created to generate and run the FEM model automatically to initialize the scaling process and during the optimization, respectively. In what concerns the output files: outFile2.txt provides the same modal information as before; while outF_it.txt besides storing the mode shapes and natural frequencies also stores information regarding strain energy for the topology optimization.

2.3. Topology Optimization

Topology optimization is used to determine an internal topological layout for the scaled model that minimizes the relative differences between the first five natural frequencies of full and scaled models. The mathematical statement of this optimization problem is the following:

$$\begin{array}{cccc}\hfill & \mathrm{minimize}\hfill & \hfill & f\left(\mathbf{x}\right)=\sum _{i=1}^m{\left(1-\frac{{\omega }_{s_i}^2\left(\mathbf{x}\right)}{{\omega }_{t_i}^2}\right)}^2\hfill \\ \hfill & \mathrm{w}.\mathrm{r}.\mathrm{t}\hfill & \hfill & 0⩽\mathbf{x}⩽1,\hfill \end{array}$$

(4)

where $\mathbf{x}$ is the vector of design variables, $\omega$ denotes natural frequency, subscripts s and t refer to scaled and target models, respectively, and m represents the number of mode shapes considered, which in this case is five. The vector of design variables ($\mathbf{x}$) represents the description of material, which follows the RAMP (Rational Approximation of Material Properties) [30] method for both the elasticity modulus (E) and the Poisson ratio ($\nu$) and the SIMP (Simplified Isotropic Material with Penalization) method [31] for the density ($\rho$), as follows:

$$E\left(x_j\right)=E_0+\frac{x_j}{1+p\left(1-x_j\right)}\left(E_m-E_0\right);$$

(5)

$$\nu \left(x_j\right)={\nu }_0+\frac{x_j}{1+p\left(1-x_j\right)}\left({\nu }_m-{\nu }_0\right);$$

(6)

$$\rho \left(x_j\right)={\rho }_0+x_j^p\left({\rho }_m-{\rho }_0\right).$$

(7)

In the above equations: j denotes an arbitrary patch or element j; p is the penalization factor; the subscript m represents the material chosen for the scaled model; $E_0$, ${\nu }_0$ and ${\rho }_0$ are properties corresponding to a material with low stiffness and mass. Each $x_j$ can assume a value within 0 (void) and 1 (solid), although to avoid singularities a minimum value of 0.0001 is established. In all the optimizations performed, the initial value for each $x_j$ is set to 0.5. The penalization factor is defined by default as 1.

Since the mode shape sequence can change during the topology optimization process [32], a mode tracking procedure based on the MAC should be employed. To this end, the first eight eigenvectors (instead of using only five) are used to compute the MAC metric and sort the modes such that on the diagonal terms have the highest correlation possible.

The Interior Point OPTimizer (IPOPT) algorithm is used to solve this topology optimization problem. IPOPT is a gradient-based algorithm developed by [33] that is able to deal with large nonlinear optimization problems, such as the one of this work. Since IPOPT is a gradient-based optimizer it requires computing the sensitivities with respect to the design variables. The natural frequencies are computed by solving the eigenvalue problem:

$$\left(\mathbf{K}-{\omega }_i^2\mathbf{M}\right){\Psi }_i=0\mathrm{with}i=1,\cdots ,m.$$

(8)

In the above equation both the stiffness matrix $\mathbf{K}$ and mass matrix $\mathbf{M}$ depend on the design variables $\mathbf{x}$, thus also the eigenvalues $\omega$ and eigenvectors $\Psi$ depend on the design variables. By differentiating the eigenvalue problem with respect to the design variables, the following is obtained:

$$\left(\frac{\partial \mathbf{K}\left(\mathbf{x}\right)}{\partial \mathbf{x}}-\left({\omega }_i^2\frac{\partial \mathbf{M}\left(\mathbf{x}\right)}{\partial \mathbf{x}}+\frac{\partial {\omega }_i^2\left(\mathbf{x}\right)}{\partial \mathbf{x}}\mathbf{M}\right)\right){\Psi }_i+\left(\mathbf{K}-{\omega }_i^2\mathbf{M}\right)\frac{\partial {\Psi }_i}{\partial \mathbf{x}}=0.$$

(9)

After some algebraic manipulation the following expression emerges:

$$\frac{\partial {\omega }_i^2\left(\mathbf{x}\right)}{\partial \mathbf{x}}={\Psi }_i^T\left(\frac{\partial \mathbf{K}\left(\mathbf{x}\right)}{\partial \mathbf{x}}-{\omega }_i^2\frac{\partial \mathbf{M}\left(\mathbf{x}\right)}{\partial \mathbf{x}}\right){\Psi }_i.$$

(10)

3. Results

In order to evaluate the proposed methodology adequacy for performing dynamic scaling using a single material topology optimization strategy, eight test cases are conducted. These test cases are also important to identify the proposed dynamic scaling methodology limitations. The first two cases, presented in Section 3.1, consists of simply introducing a difference between the material properties of both models (full and scaled), which increases from the first to the second case. Then, to assess the impact of design space size the two different grids mentioned before are compared in Section 3.2. Finally, the strategy of increasing the penalization factor to avoid intermediate densities is discussed in Section 3.3 using 4 different values.

3.1. Material Influence Test Case

Before defining this test case, it is worth to recall that the full scale model is composed of two parts, an external shell (Part-2) and an internal solid part (Part-1) made of two materials (A and B), as illustrated in Figure 2; while the scaled model consists of an internal part made of a single material and an external shell. Furthermore, the scaled model shell remains unmodified during the optimization process.

As mentioned earlier the goal of this sub-section is to evaluate the influence of varying the difference between the material properties of both models. Two case studies, with increasing difficulty, are conducted for this purpose, using the most coarse design grid (18, 24, 1). In the first case study (Case 1), the proposed dynamic scaling methodology is tested considering small differences between the material properties of the full scale model and the scaled one, as presented in

Table 1. Comparison between the material properties of the target and scaled models for Cases 1 and 2.

		Case 1			Case 2
	Material Property	Target	Scaled Model	Difference [%]	Target	Scaled Model	Difference [%]
	Density $\rho [\mathrm{kg}/{\mathrm{m}}^3]$	2400	2700	11	2400	2700	11
Part-1, Solid A	Elasticity modulus $E\left[\mathrm{Pa}\right]$	$60\times {10}^9$	$70\times {10}^9$	14	$60\times {10}^9$	$70\times {10}^9$	14
	Poisson Ratio $\nu [-]$	0.33	0.30	9	0.33	0.30	9
	Density $\rho [\mathrm{kg}/{\mathrm{m}}^3]$	2200	2700	19	2000	2700	26
Part-1, Solid B	Elasticity modulus $E\left[\mathrm{Pa}\right]$	$66\times {10}^9$	$70\times {10}^9$	6	$55\times {10}^9$	$70\times {10}^9$	21
	Poisson Ratio $\nu [-]$	0.36	0.30	17	0.36	0.30	17
	Density $\rho [\mathrm{kg}/{\mathrm{m}}^3]$	2200	1800	18	1210	1800	33
Part-2, Shell	Elasticity modulus $E\left[\mathrm{Pa}\right]$	$66\times {10}^9$	$40\times {10}^9$	39	$66\times {10}^6$	$40\times {10}^9$	60,506
	Poisson Ratio $\nu [-]$	0.36	0.33	8	0.36	0.33	8

. These difference are then increased for the second case (Case 2) as can be observed also from Table 1. In what concerns solid A the same differences between material properties as for Case 1 are maintained, however for solid B these differences are increased and consequently making it more difficult to achieve dynamic similarity. Since unitary scaling factors are defined for both density and pressure the material properties of the target and full size models are the same.

For Case 1, the highest differences in terms of elasticity modulus is observed for the shell part ($39\%$); while for the density and Poisson ratio, the solid part made of material A presents the highest differences ($19\%$ and $17\%$, respectively). It is important to mention that although the shell thickness is 1 mm in the full size model, its reduced size counterpart is four times higher (i.e., 4 mm) to further difficult the dynamic scaling process. This last consideration is also applied in Case 2. Even though a difference of 60,506% in the elasticity modulus is considered for the shell in Case 2, its influence on the results is reduced when compared with the internal solids A and B given its small thickness value.

The optimization for Case 1 reached a reasonable good solution just after performing 16 iterations, as one can see from the convergence history in Figure 6. After these 16 iterations, despite the solution fluctuates slightly the best result is reached at the 26th iteration with an objective function value of $2.80\times {10}^{-4}$. As expected, Case 2 presents a slightly slower convergence than Case 1 given the higher differences in material properties between not only target and scaled models but also between solids A and B. The lowest objective function value ($3.78\times {10}^{-4}$) is reached at 49th iteration. For both cases, the maximum number of iterations is set to 60.

These low objective function values indicate that a very small difference between natural frequencies is reached, which is demonstrated in

Table 2. Comparison between the target frequencies and obtained ones for Cases 1 and 2. F denotes frequency.

	Case 1			Case 2
	Target F [Hz]	Obtained F [Hz]	Relative Error [%]	Target F [Hz]	Obtained F [Hz]	Relative Error [%]
Mode 1	61.00	61.06	0.01	59.24	58.90	0.06
Mode 2	263.77	264.76	0.04	257.61	257.54	0.00
Mode 3	378.81	381.53	0.07	368.66	368.77	0.00
Mode 4	656.28	655.54	0.01	643.87	640.55	0.05
Mode 5	918.68	917.47	0.01	897.06	891.73	0.06

. In fact, the relative error between the target and obtained frequencies is lower than 0.1% for both case studies.

To verify if the mode shapes from both models also match, one can use the previously mentioned MAC metric. Despite not directly expressed in the objective function, the MAC is used to track the right mode shapes to compare the associated natural frequencies. This metric is illustrated for the best result of Case 1 in Figure 7a, where the less correlated mode shape (5th mode shape) presents $94\%$ of correlation concerning its peer (diagonal values), and higher correlation between each other (off-diagonal values) was $10.6\%$. A similar correlation is found for Case 2 as depicted in Figure 7b.

The results obtained indeed indicate that the matching of both the natural frequencies and mode shapes is successfully accomplished even for the most coarse design grid (18, 24, 1). This proves that the developed methodology can be used for dynamic scaling with similar assumptions to those considered in this case study. For this case study, the material properties of the solid parts in the scaled model differ only slightly from the target material properties, which are equal to ones of the full scale model since unitary density and pressure scaling factors are defined. The correlation between mode shapes degrades slightly with the increase of these differences in material properties of the solid parts.

Regarding the topologically optimized internal structures, the solutions for both cases are very similar as illustrated in Figure 8.

3.2. Design Space Test Case

The objective now is to evaluate the impact of increasing the number of design variables on dynamic similarity for a more challenging test case. For this purpose the two design grids presented in Section 2.2 are used. In the previous case studies, the material properties were fictitious and adapted to analyze differences between models. On the other hand, in the current test case existing materials are set for both full and scaled models, which increases the differences in material properties. To increase the difficult of reaching dynamic similarity, the full scale model is designed to be made of a 5000 Series Aluminum Alloy for Part-1 Solid A, Nylon 12 (PA12) for Part-1 Solid B and Nylon 610 (PA610) for the shell (Part-2). Regarding the scaled model, both solids are composed of a general Aluminum alloy and the shell is set to be made of Nylon (PA12). It is worth to note that Nylon (PA12) is considered for both full and scaled models, although for different parts (Solid B and shell for the former and latter, respectively) and considering slightly differences in the properties. Furthermore, it is also important to mention that the material selected for the skin of the scaled model presents higher mechanical properties than those of the equivalent full size counterpart. Even though the shell is expected to have a slightly higher impact on the dynamic scaling process than before, the internal structure still presents higher density and elasticity modulus. The material properties used for designing both models and design grids are shown

Table 3. Comparison between the material properties of the target and scaled models for Cases 3 and 4.

		Target	Scaled Model	Difference [%]
	Density $\rho [\mathrm{kg}/{\mathrm{m}}^3]$	2400	2700	13
Part-1, Solid A	Elasticity modulus $E\left[\mathrm{Pa}\right]$	$60\times {10}^9$	$70\times {10}^9$	17
	Poisson $\nu [-]$	0.33	0.30	9
	Density $\rho [\mathrm{kg}/{\mathrm{m}}^3]$	1800	2700	50
Part-1, Solid B	Elasticity modulus $E\left[\mathrm{Pa}\right]$	$30\times {10}^9$	$70\times {10}^9$	133
	Poisson $\nu [-]$	0.36	0.30	17
	Density $\rho [\mathrm{kg}/{\mathrm{m}}^3]$	1210	1800	49
Part-2, Shell	Elasticity modulus $E\left[\mathrm{Pa}\right]$	$66\times {10}^6$	$40\times {10}^9$	60,506
	Poisson $\nu [-]$	0.36	0.33	8

. Cases 3 and 4 denote coarse and fine design grids, respectively. The highest difference in E is by far for the shell (60,506%). In regard to the $\rho$ and $\nu$, the Part-1 Solid B presents the highest differences. For this internal part, the scaled model is defined with a material $50\%$ denser and $133\%$ stiffer than the full scaled model.

The convergence history of the objective function from both design grids is not very different, as shown in Figure 9. None of the optimization problems converged and both cases stopped by the maximum number of iterations criterion (200). Nevertheless, reasonably good solutions are reached after 90 and 70 iterations when using coarse (18, 24, 1) and fine (54, 72, 1) design variable grids, respectively. After these iterations, small fluctuations are observed. The best solution for the Case 3 (18, 24, 1) is reached at 199th iteration with value of $7.78\times {10}^{-4}$, while for Case 4 the lowest objective function value ($4.99\times {10}^{-4}$) is obtained at the 138th iteration. If more iterations were performed, a better objective function value could have been obtained and even convergence could have been reached. Nevertheless, since the objective function values are already low (in the same magnitude as before) and given the higher computational cost, the maximum number of iterations was maintained.

The low objective function values are an indicative of a good matching between full and scaled models in terms of natural frequencies, which is confirmed in

Table 4. Comparison between the target frequencies and obtained ones for Cases 3 and 4. F denotes frequency.

		Case 3		Case 4
	Target F [Hz]	Obtained F [Hz]	Relative Error [%]	Obtained F [Hz]	Relative Error [%]
Mode 1	52.31	52.72	0.87	51.94	0.71
Mode 2	234.32	235.43	0.47	235.51	0.50
Mode 3	321.93	324.67	0.85	321.07	0.27
Mode 4	592.69	595.44	0.46	593.84	0.19
Mode 5	786.21	785.53	0.09	791.05	0.61

. In relation to the previous test case, now the relative differences are higher, but still lower than $1\%$. This outcome is associated with the higher discrepancies between the material properties, i.e., further way from the target ones. As expected, the higher design freedom of the finer grid (Case 4) allowed for a lower relative errors.

Regarding the mode shapes a good correlation between models is also verified as illustrated in Figure 10. As noticed for the natural frequencies, the mode shape correlation is also reduced for these case studies when compared with the previous ones. For Case 3 the smallest on-diagonal correlation is $87\%$ for the 4th mode shape and the highest off-diagonal value is $15\%$. The correlation slightly decreased for the finest design grid used.

These results demonstrate that the proposed methodology is able to successfully design a model numerically for dynamic scaling. However, this almost perfect matching of the natural frequencies and very good agreement between mode shapes were possible without using a filter in the topology optimization process and the penalization factor was set 1. As consequence, the solution presents intermediate densities as can be seen in Figure 11.

For Case 4, most of the intermediate densities were close to 1 or $0.5$, as one can see from Figure 12a. When analyzing the optimized internal structure, it is worth to recall that the scaled model has a shell made of a material substantially stiffer than the target model and its thickness is 4 times higher. The solution was thus affected by the shell in such a way that the amount of material near the root is smaller compared to other wing locations. Despite the obtained solution does not present checkerboard problem, it is possible to see some “isolated” cells in Figure 12a.

As noted before, it is possible to observe from Figure 12a that most of the cells from Case 4 solution present values either close to 1 or $0.5$. This particular outcome is very interesting since it allows to consider a post-processing, such that these intermediate densities can be assigned to another material. The good matching of the natural frequencies and mode shapes makes this approach even more promising since it would eventually allow it to be manufactured. Moreover, if performing more iterations the results could have been improved. In this sense, in Figure 12b is presented as a post-processing, where black and red cells correspond to the original and second material, respectively. Such intermediate material was found possible for this case after making a search in [34], although it might be infeasible for other cases. The material in question is the Nylon 12 (PA12). From the 3D printing perspective, the challenge consists of defining materials that can be printed together. However, the topic of metal-plastic hybrid 3D printing is evolving [35,36].

Despite the task of designing scaled model with the internal structure made of two materials is feasible, it is not the focus of this paper. Here, the strategy followed to avoid intermediate solutions is the usage of a penalization factor in the next sub-section.

3.3. Penalization Factor Test Case

The penalization factor influence on the topologically optimized structures is now the focus of study in this test case. Four different penalization values are used for the finest considered design grid (54, 72, 1): 1 (Case 5), 2 (Case 6), 3 (Case 7) and 6 (Case 8). Despite the same full size model as for Cases 3 and 4 is used in these case studies, the materials and shell thickness of the scaled model are now different. The material properties of both models are summarized in

Table 5. Comparison between the material properties of the target and scaled models for Cases 5, 6, 7 and 8.

		Target	Scaled Model	Difference [%]
	Density $\rho [\mathrm{kg}/{\mathrm{m}}^3]$	2400	4430	46
Part-1, Solid A	Elasticity modulus $E\left[\mathrm{Pa}\right]$	$60\times {10}^9$	$113.8\times {10}^9$	47
	Poisson $\nu [-]$	0.33	0.342	4
	Density $\rho [\mathrm{kg}/{\mathrm{m}}^3]$	1800	4430	59
Part-1, Solid B	Elasticity modulus $E\left[\mathrm{Pa}\right]$	$30\times {10}^9$	$113.8\times {10}^9$	74
	Poisson $\nu [-]$	0.36	0.342	5
	Density $\rho [\mathrm{kg}/{\mathrm{m}}^3]$	1210	1210	0
Part-2, Shell	Elasticity modulus $E\left[\mathrm{Pa}\right]$	$66\times {10}^6$	$66\times {10}^6$	0
	Poisson $\nu [-]$	0.36	0.33	8

. A Titanium alloy (Ti-6Al-4V (Grade 5)), often used for manufacturing through 3D printing [37], is considered for the internal part (solids A and B) of the scaled model. This pose an even higher challenged than before since the difference between properties increases. Regarding the shell, it is now 1 mm thick, i.e., equal to the full size model. This choice is related to the trend of reducing material near the root observed in Figure 12.

Given the high computational cost of these numerical optimizations, the maximum number of iterations is set to 135. A stable convergence history of the objective function can be seen in Figure 13 for the 4 considered penalization factors. As expected, the proposed dynamic scaling methodology takes more iterations to reduce the objective function when a higher penalization factor is considered since the design freedom is reduced. The lowest obtained objective function values are: $1.7412\times {10}^{-4}$ at iteration 120 for Case 5; $6.1544\times {10}^{-5}$ at iteration 128 for Case 6; $8.6974\times {10}^{-5}$ at iteration 64 for Case 7; and $1.1874\times {10}^{-3}$ at iteration 97 for Case 8.

These relatively low objective function values indicate that the first five natural frequencies are very close to the target ones, which is confirmed in

Table 6. Comparison between the target frequencies and obtained ones for Cases 5, 6, 7 and 8. F denotes frequency. The errors are relative to the target.

		Case 5: Pen = 1		Case 6: Pen = 2		Case 7: Pen = 3		Case 8: Pen = 6
	Target	Obtained	Relative	Obtained	Relative	Obtained	Relative	Obtained	Relative
	NF [Hz]	F [Hz]	Error [%]	F [Hz]	Error [%]	F [Hz]	Error [%]	F [Hz]	Error [%]
Mode 1	52.31	52.58	0.52	52.47	0.31	52.42	0.21	52.67	0.68
Mode 2	234.32	234.61	0.12	234.67	0.15	234.87	0.23	236.63	0.98
Mode 3	321.93	322.79	0.27	322.47	0.17	322.61	0.21	323.22	0.40
Mode 4	592.69	591.37	0.22	592.27	0.07	593.74	0.18	592.33	0.06
Mode 5	786.21	787.46	0.16	786.57	0.04	787.84	0.21	795.34	1.16

. Even though higher objective function values and relative errors of frequencies are in general obtained by increasing the penalization factor, the relative errors are still low. In fact, the highest error (1.16%) is determined for the fifth vibration mode considering the highest penalization factor (Case 8). Thus, showing the proposed methodology effectiveness in matching the natural frequencies.

The correlation between target and scaled mode shapes is again assessed through the MAC metric, which is depicted in Figure 14 for the considered cases. This correlation is clearly observed to degrade as the penalization factor is increased. When using the lowest value (Case 5), a very good agreement between mode shapes is found since the lowest correlation is $94\%$ and the highest off-diagonal term is 16%, as illustrated in Figure 14a. Both these values are within the range usually deemed as acceptable, $80\%$ and $20\%$ for correlation and off-diagonal terms, respectively, [38]. By increasing the penalization factor to the double (Case 6, Figure 14b), despite an off-diagonal term of $24\%$ is calculated, the lowest correlation ($91\%$) is still considerably higher than the above mentioned threshold. For this reason, Case 6 is still considered a good result by the authors. However, for the Case 7 (Figure 14c) these metrics continues degrading and the lowest correlation falls to $82\%$, which is barely above the threshold value, but the maximum off-diagonal term slightly raises to $26\%$. Given the fact that at least the correlation is still above the threshold value, this case is considered acceptable by the authors although not desirable. For Case 8 (Figure 14d) both these thresholds are not satisfied for the first five mode shapes, thus the mode shape matching is deemed not achieved even though the frequencies present low relative errors.

The topologically optimized internal structures for the different cases are shown in Figure 15. It is possible to identify areas with different intermediate density in all the cases. Despite increasing the penalization factor has not allowed for completely removing all intermediate densities, its amount has reduced. The identification of well-defined areas with intermediate densities makes the post-processing with different materials an attractive strategy. However, this approach might not be possible for other cases, e.g., the resulting intermediate material properties can be nonexistent. Furthermore, even if the second material exists it might not be compatible with the original one. Thus, a multi-material topology optimization can be a promising strategy for dynamic scaling.

The studies considering different penalties enabled important remarks. First, the increasing of the penalization negatively affected the quality of results in such a way that the optimizer failed in performing the matching of 5th mode shape for Case 8. The same effect was observed for the matching of the frequencies, although the worst result was $1.16\%$ which is fully accepted. Even the highest penalty factor not allowed achieving a gradient-free solution. However, observing how the penalty affected the results, a suitable solution of such optimization problems considering just gradient-free designs is really a very hard task. For instance, the matching of both natural frequencies and mode shapes with the same quality as reached with $\mathrm{Pen}=1$ maybe impossible to be achieved through a fully gradient-free approach.

4. Concluding Remarks

In this paper, a dynamic scaling methodology using topology optimization was introduced with the aim at designing scaled models that can be used to meaningfully evaluate unconventional designs for different engineering applications (e.g., aerospace, civil and mechanical engineering). Eight case studies of increasing difficulty were conducted to assess its effectiveness in representing the same dynamic behavior, in terms of natural frequencies and associated mode shapes, of a full size model in a reduced scale one. This methodology was found very promising in matching both the natural frequencies and modes shapes when the differences between material properties are small. The MAC metric was found to be close to an identity matrix and the relative difference between natural frequencies is under 1%. As these differences increase, even though the correlation between mode shapes from full and scaled models degrades (some diagonal terms become lower than 90%), a very good agreement was still reached. Nevertheless, the obtained solutions were characterized by presenting intermediate densities, even after increasing the design space, which are not suitable for additive manufacturing. Thus, a penalization factor was added to reduce these intermediate densities. When adding this factor, despite very low discrepancies were found for the frequencies (lower than $1.2\%$), the correlation between mode shapes noticeably degrades, with off-diagonal terms above 20% for penalization factors above 2. In fact, for the last considered case the thresholds usually deemed acceptable for the MAC metric were passed.

A solution that could mitigate this issue is using a multi-material topology optimization strategy for the dynamic scaling process. Alternatively, the MAC metric could be included in the objective function, while accounting for manufacturing feasibility constraints.