**Numerical Analysis of the Influence of Empty Channels Design on Performance of Resin Flow in a Porous Plate**

**Glauciléia M. C. Magalhães 1, Cristiano Fragassa 2,\*, Rafael de L. Lemos 3, Liércio A. Isoldi 1,3, Sandro C. Amico 4, Luiz A. O. Rocha 5, Jeferson A. Souza 1,3 and Elizaldo D. dos Santos 1,3**


Received: 19 March 2020; Accepted: 17 April 2020; Published: 11 June 2020

**Featured Application: The influence of empty channel design on resin flow in a porous plate is studied. Constructal design is used for geometric investigation of I and T-shaped channels. The multiphase resin**/**air flow in the porous plate is solved numerically. T-shaped configurations with long branches lead to the best impregnation performance. Permanent air voids arise in the domain for very thin T-shaped channels.**

**Abstract:** This numerical study aims to investigate the influence of I and T-shaped empty channels' geometry on the filling time of resin in a rectangular porous enclosed mold, mimicking the main operating principle of a liquid resin infusion (LRI) process. Geometrical optimization was conducted with the constructal design (CD) and exhaustive search (ES) methods. The problem was subjected to two constraints (areas of porous mold and empty channels). In addition, the I and T-shaped channels were subjected to one and three degrees of freedom (DOF), respectively. Conservation equations of mass and momentum for modeling of resin/air mixture flow were numerically solved with the finite volume method (FVM). Interaction between the phases was considered with the volume of fluid method (VOF), and the effect of porous medium resistance in the resin flow was calculated with Darcy's law. For the studied conditions, the best T-shaped configuration resulted in a filling time nearly three times lower than that for optimal I-shaped geometry, showing that the complexity of the channels benefited the performance. Moreover, the best T-shaped configurations were achieved for long single and bifurcated branches, except for configurations with skinny channels, which saw the generation of permanent voids.

**Keywords:** constructal design; resin flow; porous media; numerical simulation; filling time

#### **1. Introduction**

Liquid composite molding (LCM) processes were developed to enable the production of large, wide, and complex structural components at low cost. Among the important LCM processes, it is

important to mention liquid resin infusion (LRI), which consists of the injection of a polymer resin through empty channels mounted on a fibrous mold, favoring the global propagation of the resin over the mold domain [1]. One of the most considerable difficulties in applying LRI is related to complete mold filling and void formation, ensuring that the fibrous reinforcement is completely impregnated by the resin within the mold [2–5].

To tackle the reported difficulties, several efforts have been made to improve process comprehension, as well as to avoid traditional trial-and-error approaches. One possibility is the use of computational tools to simulate multiphase resin/air flow into the enclosed mold [3–8]. The study of parametrical investigations in the LRI process (as the influence of design over the behavior of resin/air flow, voids formation, and filling time) also represents an important subject. Therefore, studies have been carried out in the analytical, experimental, and numerical scope to improve the understanding of resin/air flow in this kind of domain [3,4].

The development of computational methods for the representation of LCM processes is an important point of concern. For example, Isoldi et al. [9] performed a numerical study of resin flow in porous domains that mimics the operational principle of the Resin Transfer Molding (RTM) and Light Resin Transfer Molding (LRTM) processes. Different resin flow cases were solved numerically (linear, radial, and complex three-dimensional domain) and compared with analytical, experimental, and numerical results available in the literature. The maximum difference between these solutions was about 8.0%, showing the validity of the developed numerical method. Wang et al. [1] conducted a numerical and experimental study of the resin infusion process under industrial conditions. The authors verified that mold filling time determined by numerical simulation brought about results in agreement with experimental solutions. Afterward, Grössing et al. [10] compared the technical feasibility of using two different numerical methods (PAM-RTM® and OpenFOAM) to obtain numerical solutions of resin flow. Firstly, it was found that the numerical methods led to valid solutions in comparison with experimental results. It was also verified that both numerical methods had very similar results for the simulation of resin flow in porous media. Sirtautas et al. [11] developed a numerical model based on the Finite Element Method (FEM), considering the effects of compressibility and orthotropic permeability dependent on pressure in the domain for different fibrous materials. A comparison between the numerical model and experimental results was performed for resin flow in the two-dimensional domain of an aerospace piece. The results showed good agreement, recommending the use of the proposed computational method. Pierce et al. [12] developed and validated a multi-physics model to improve the simulation of the infusion process in the complex preform. The results are compared with traditional models that do not consider fabric deformation. The combination of deformation-dependent permeability properties and the preform draping model led to more realistic predictions for local infusion flow. Chebil et al. [13] developed a computational model for the simulation of three-dimensional resin flow in laminated preform composed of multiple layers with different permeabilities. The method proposed the use of shell elements for flow simulation instead of solid elements, leading to a reduction in computing time. Rubino and Carlone [14] proposed an analytical/numerical methodology to consider the effects of preform compaction on resin flow for a deformable porous media, which is important for processes such as vacuum assisted resin transfer molding (VARTM) and Seeman composites resin infusion molding process (SCRIMP), where a flexible bag is generally used as part of the mold. The results showed that preform compaction affected resin flow and filling time.

The numerical approach has also been used to define strategies for several LCM processes. Examples of the evaluation of gate placements in RTM processes seeking to prevent void content and minimize filling time had been widely reported in different works [15–20]. The work of Brouwer et al. [3] proposed strategies to improve the positioning of injection channels in one LRI process using computational modeling. More precisely, the authors presented developments in two large structures to be filled with resin, a rotor blade, and a boat hull 20 m and 16 m long, respectively. The strategy adopted by the authors consisted of a main injection channel in the central region of the domain

(the keel for the case of injection in the hull domain) with several bifurcated channels spread along the area of the domain (from the center to the periphery). Afterward, Gomes et al. [21] performed a numerical and experimental study for better characterization of a VARTM (vacuum assisted resin transfer molding) process, indicating the best strategies for the filling of a C-shaped stringer. The results indicated that the consideration of two exit regions (with vacuum) and one entering region of resin placed in the extremes of domain length was the best strategy to minimize filling time and prevent void formation. Wang et al. [1] developed an iterative method based on a centroidal Voronoi diagram for optimization of the distribution of several channels to reduce the number of required tests.

Recently, Struzziero and Skordos [22] performed a numerical multi-objective optimization focusing on the choice of an optimal temperature profile, which minimizes the filling time and the risk of hindering resin flow due to excessive curing. For the optimization, the authors employed the genetic algorithm (GA) method. It achieved reductions of nearly 66% and 15% in filling time and final degree of cure, respectively, in comparison with standard solutions. Geng et al. [23] investigated the behavior of resin impregnation in curved porous plates mimicking a VARTM process. They performed an experimental study investigating the influence of curvature angles and preform layers on the advancement of the resin front line. The results indicated that due to interaction between the preforms and void formation, the curved regions led to an augmentation of 15% to 30% of filling time in comparison with a linear region. Shevtsov et al. [24] performed an experimental and numerical study for the manufacturing of three-dimensional composite parts of complex shape. For the computational model, a numerical technique based on liquid resin flow into porous preform is coupled with a model that describes the fluid motion in unsaturated soils. The developed method allowed the prediction of inner and outer dry spots depending on the resin flow movement. Despite several important works reported, few studies can be found on the evaluation of the influence of the geometry of empty channels in the behavior of resin/air fluid flow into the porous domain. In the present study, a geometrical optimization of I and T-shaped empty channels over the required time to fill a porous rectangular plate is performed with the constructal design (CD) and exhaustive search (ES) methods.

Constructal design is a method for assessment of design in any finite flow system based on the objectives and constraints principle. This method has been widely applied to show that shape and structure in nature (trees, rivers, animal body, and others) are deterministically generated following a physical principle named Constructal Law of Design and Evolution [25–29]. According to Bejan [28], constructal law states that the design of the finite-size flow system to persist in time (to live) must evolve freely to maximize access to its internal currents. This method has also been applied as a powerful tool to improve the design of several engineering problems as cooling systems, heat exchangers, renewable energy, transport systems, and even solid mechanics [30–36]. As per the authors' knowledge, except for the works of Refs. [37,38], the constructal design method has not been employed in the literature for investigation on the influence of geometry over impregnation of resin in a porous medium, which mimics the LRI process, a novelty in the present work.

Here, the application of computational modeling, constructal design, and exhaustive search (ES) in geometrical optimization of I and T-shaped empty channels inserted into the rectangular porous medium is proposed. The main aims here are to minimize the filling time of resin impregnation in the porous medium, prevent the formation of permanent voids in the domain, and investigate the influence of geometry over the performance indicator of the system. The solution domain is presented in a two-dimensional rectangular plate with a porous medium, which simulates a fibrous reinforcement. This work is a sort of continuity of the study by Magalhães et al. [38] with the difference being that here the T-shaped channel is completely optimized, i.e., all three degrees of freedom are investigated. The conservation equations of mass, momentum, and one equation for the transport of volumetric fractions are solved with the finite volume method (FVM) [39–41]. To deal with the resin/air interface flow, the volume of fluid (VOF) method is used [42]. In the region of the porous medium, it is also considered a body force modeled with Darcy's law. More precisely, the simulations were performed using the computational fluid dynamics code FLUENT [41].

#### **2. Mathematical Modeling**

In this study, incompressible, laminar, and transient flow of a resin/air mixture in a two-dimensional domain is considered. Moreover, the two phases are treated as immiscible. For the prediction of this kind of flow, conservation equations of mass and momentum for the mixture of resin/air and one transport equation for prediction of resin volume fraction are numerically solved. The conservation equations of mass and momentum for the resin/air mixture are given by [43]:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho \overrightarrow{V}\right) = 0 \tag{1}$$

$$\frac{\partial \left(\rho \overrightarrow{V}\right)}{\partial t} + \nabla \cdot \left(\rho \overrightarrow{V} \overrightarrow{V}\right) = -\nabla P + \nabla \cdot \overrightarrow{\tau} + \overrightarrow{F} \tag{2}$$

where ρ is the mixture density (kg/m3), → *V* is the velocity vector (m/s), *t* is the time (s), μ is the dynamic viscosity of the fluid (Pa.s), <sup>∇</sup>*<sup>P</sup>* is the pressure gradient (Pa/m), <sup>→</sup> *F* is an external force vector per unit volume (N/m3) and <sup>=</sup> τ is the stress tensor of the fluid (N/m2), given by:

$$
\overline{\tau} = \mu \left( \nabla \overrightarrow{V} + \nabla \overrightarrow{V}^T \right) \tag{3}
$$

The effect of the porous medium is included in the mathematical model by the insertion of a body force in the momentum equation based on Darcy's law, as given by [43–45]:

$$
\overrightarrow{F} = -\frac{\mu}{K}\overrightarrow{V} \tag{4}
$$

To tackle the resin/air mixture, the volume of fluid (VOF) method is employed [42]. In this approach, an additional transport equation for resin volume fraction (*f*) is necessary to define the quantity of resin along with each cell of the domain. This transport equation is given by [42]:

$$
\frac{\partial f}{\partial t} + \nabla \cdot \left( f \overrightarrow{V} \right) = 0 \tag{5}
$$

With the definition of volume fraction, density and dynamic viscosity in each cell of the computational domain can be calculated by [46]:

$$
\rho = f \rho\_{\text{res}} + (1 - f) \rho\_{\text{air}} \tag{6}
$$

$$
\rho = f \mu\_{res} + (1 - f)\mu\_{air} \tag{7}
$$

#### *2.1. Problem Description*

As previously mentioned, an incompressible, laminar, and transient multiphase flow of resin/air was considered, moving through a two-dimensional plate domain composed of a porous medium, as shown in Figure 1a,b. Empty channels are inserted along the porous domain to facilitate impregnation of the resin throughout the mold. In Figure 1a,b, these channels (I and T-shaped, respectively) are represented by the dark gray region, while the porous region is represented by light gray.

**Figure 1.** Computational domain of the problem, which mimics the liquid resin infusion (LRI) process with two empty channels geometries: (**a**) I-shaped, (**b**) T-shaped.

In both cases, flow is imposed by a pressure difference between the resin inlet (lower surface of the empty channel) and the exit region of multiphase flow (upper surface of the porous medium). Inlet pressure of *Pin* <sup>=</sup> <sup>1</sup> <sup>×</sup> 105 Pa is considered, and pressure of *Pout* <sup>=</sup> 0 Pa is imposed in the exit. It is worth mentioning that, in a real LRI process, resin flow is driven by vacuum imposition at the exit of the fibrous reinforcement. Despite this, a pressure difference with the imposition of positive pressure in the resin inlet surface is represented. For computational modeling, it is important to predict the pressure difference. Hence, the imposition made here led to the same results as that achieved when the same pressure difference is imposed with the vacuum in the domain exit. In the remaining domain, bottom, and lateral surfaces, a non-slip and impermeability boundary condition is imposed (*u* = *v* = 0 m /s). For the resin volume fraction (Equation (5)), *f* = 1 at the inlet section and zero gradients (normal to the surface) on all the other surfaces are imposed. Concerning the physical properties of the resin and air, densities of ρ*res* = 916 kg/m<sup>3</sup> and ρ*air* = 1.225 kg/m3 and dynamic viscosities of μ*res* = 0.06 kg/(ms) and <sup>μ</sup>*air* <sup>=</sup> 1.7894 <sup>×</sup> <sup>10</sup>−<sup>5</sup> kg/(ms) are adopted. For the porous media, a permeability of *K* = 3.0 <sup>×</sup> 10−<sup>10</sup> m2 and porosity of ε = 0.88 is imposed.

#### *2.2. Geometrical Evaluation*

Geometrical optimization was done by using an association between constructal design and exhaustive search, as explained in recent literature [47]. Figure 2 illustrates the main steps employed for optimization of the problem. Steps 1 to 3 presents the flow system and identifies a performance indicator (here, it is filling time). Steps 4 and 5 show the constraints of the problem, degrees of freedom, and the main parameters of the physical problem. It is also necessary to understand the physical problem and define a method for its solution (numerical solution in this case).

With the defining of geometry and understanding of the physical solution method, it is possible to predict the performance indicator (minimization of filling time of resin impregnation along with the porous domain) for each configuration studied. After Step 6, it is necessary to define if the design should be optimized or not.

If the problem is not optimized, a geometrical evaluation is performed, seeking to understand the effect of degrees of freedom over performance indicator, simulating the total number of cases defined. However, if the problem is optimized, it is necessary to specify an optimization method. Here, exhaustive search, which consists of evaluating several geometrical possibilities taking into consideration fixed increments of variation for each degree of freedom (as depicted in Figure 3) for the T-shaped case.

**Figure 2.** Flow chart illustrating the main steps of the geometrical optimization of I and T-shaped empty channels inserted into the rectangular porous plate for resin impregnation along the domain.

**Figure 3.** Flowchart of the performed simulations in geometric evaluation of the T-shaped empty channel in resin infusion simulation.

For complete optimization, Steps 6b to 9b, illustrated in Figure 2, are executed. It should be emphasized that heuristic-based optimization methods, e.g., genetic algorithm and simulated annealing, have also been associated with constructal design [48]. These kinds of tools are more used in software packages such as FLUENT [41] and Magma5 [49].

These methods required a significantly lower number of cases to be solved, but need verification with exhaustive search for suitable reproduction of the effects of DOFs over the performance indicator.

Each case under analysis is submitted to two restrictions, the total area of the computational domain:

$$A = HL\tag{8}$$

and the area of the empty channel, which for I-shaped and T-shaped channels, respectively, are represented by:

$$A\_{\rm I} = H\_0 L\_0 \tag{9}$$

$$A\_T = H\_0 L\_0 + H\_1 L\_1 \tag{10}$$

The dimensions of the external domain with porous medium are: *H* = *L* = 0.5 m for all simulations, i.e., the outer area is a geometric constraint of the problem, and *H*/*L* = 1.0 is taken into account. It should be noted that a spacing of 0.05 m was taken into consideration for the sides of the fibrous reinforcement mold, and an additional height of 0.25 m was inserted to prevent void regions in the locality defined by the dashed line shown in Figure 1a,b, which is the main analysis region of this study. The areas *AI* and *AT* can be represented by the ratio between the empty channel and plate areas (fraction area), written by:

$$
\phi = \frac{A\_I}{A} = \frac{A\_T}{A} \tag{11}
$$

For the I-shaped empty channel, one degree of freedom (*L*0/*H*0) is considered, whereas the T-shaped empty channel has three degrees of freedom (*L*0/*H*0, *H*0/*H*1, and *L*1/*H*1). To optimize the geometry of the I-shaped channels, a set of 60 simulations were performed with different *L*0/*H*<sup>0</sup> ratios and different values of φ, where φ = 0.005; 0.01; 0.03; 0.05, and 0.1. In the case of T-shaped channels, only φ = 0.05 was evaluated, due to the high number of simulations required for each magnitude of φ. In this specific case, a total of 189 simulations were performed. For this case (T-shaped channel), the optimization process was divided into three steps, as shown in Figure 3.

Figure 3 illustrates that in the first optimization step, the ratio *L*1/*H*<sup>1</sup> is varied, while the remaining parameters (*L*0/*H*<sup>0</sup> and *H*0/*H*1) are kept fixed. To clarify how the terms are assigned, note that the minimum value found for the filling time is named the once minimized filling time (*t*m), and the corresponding ratio *L*1/*H*<sup>1</sup> is called the once-optimized ratio of *L*1/*H*1, (*L*1/*H*1)o. In a second step, the same process is repeated by varying *L*1/*H*<sup>1</sup> for the other ratios of *L*0/*H*<sup>0</sup> and maintaining *H*0/*H*<sup>1</sup> constant. At this stage of the process, two degrees of freedom are optimized (*L*1/*H*<sup>1</sup> and *L*0/*H*0). In this

case, the minimum filling time obtained is twice minimized (*t*mm), the *L*1/*H*<sup>1</sup> ratio is the twice-optimized ratio (*L*1/*H*1)oo, and the *L*0/*H*<sup>0</sup> ratio is once optimized (*L*0/*H*0)o. Finally, the second step is repeated for the different ratios of *H*0/*H*<sup>1</sup> evaluated. Thus, the minimum filling time is the thrice-minimized (*t*mmm), and the optimal shapes are the thrice-optimized ratio of *L*1/*H*1, (*L*1/*H*1)ooo, the twice optimized ratio of *L*0/*H*0, (*L*0/*H*0)oo, and the once-optimized ratio of *H*0/*H*1, (*H*0/*H*1)o, completing the geometric optimization process.

#### **3. Numerical Modeling and Code Validation**

For simulation of the infusion process, the conservation equations of mass, momentum, and transport equation of volume fraction were solved using a commercial CFD (computational fluid dynamics) code based on FVM [39–41]. Numerical simulations were carried out using computers with 4 Intel Xeon processors and 8 threads of 3.30 GHz clock and 8.0 Gb of RAM. In all the simulations, the second-order upwind discretization scheme was used for the treatment of the advective terms. For pressure interpolation, the PRESTO scheme was used. Pressure-velocity coupling was performed with the PISO method, while the Geo-Reconstruction method was used to reconstruct the interface between the two fluids. In addition, sub-relaxation factors of 0.3 and 0.7 are imposed for the conservation equations of mass and momentum, respectively. The simulations are considered converged when the residuals of conservation equations of mass and momentum are lower than *R* < 10<sup>−</sup>6.

For time discretization, a study of time step independence was performed, since a transient problem was being investigated. Five simulations were performed with a varying time step, and its influence over the volume of domain completely impregnated by the resin was investigated, i.e., the volume of domain with volume fraction of *f* = 1.0. For all simulations, the volume of resin injected up to time of *t* = 7.0 s was evaluated. Table 1 shows, for the analyzed injection time, the volume percentage of the domain filled with resin and the processing time required to carry out the simulations. The analyzed case consists of a I-shaped empty channel with area fraction of φ = 0.05, *H*/*L* = 1.0 and *<sup>L</sup>*0/*H*<sup>0</sup> <sup>=</sup> 4.0. It is possible to observe that for values of <sup>Δ</sup>*<sup>t</sup>* <sup>≤</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>3</sup> s, all employed time steps lead to identical results. For the values of <sup>Δ</sup>*t* = 1.0 <sup>×</sup> 10−<sup>2</sup> s and <sup>Δ</sup>*t* = 1.0 <sup>×</sup> 10−<sup>1</sup> s, there is no convergence in the simulations. Consequently, the time step of <sup>Δ</sup>*t* = 1.0 <sup>×</sup> 10−<sup>3</sup> s is adopted in subsequent geometric evaluation simulations, as it leads to a lower computational effort compared to smaller time steps.


**Table 1.** Study of sensibility of time step over resin impregnation for an I-shaped empty channel with φ = 0.05, *H*/*L* = 1, and *L*0/*H*<sup>0</sup> = 4.

\* Non-converging solutions.

For spatial discretization, the domain was divided into finite rectangular volumes, and a mesh independence test was performed to define the number of volumes used for the simulations. The time step used is defined in Table 1 (Δ*t* = 1.0 <sup>×</sup> 10−<sup>3</sup> s). Table 2 shows the number of volumes, the filling time for complete impregnation of resin in the mold, and the processing time required to perform the simulations for each evaluated mesh and for the same I-shaped channel studied in the time independence study, i.e., φ = 0.05, *H*/*L* = 1.0 and *L*0/*H*<sup>0</sup> = 4.0. The following equation presents the criterion for achievement of an independent mesh:

$$Difference = \frac{100 \cdot \left| t^j - t^{j+1} \right|}{t^j} < 1.0\% \tag{12}$$

where *t <sup>j</sup>* represents the minimum value of the filling time calculated with the coarser mesh, and *t j*+1 corresponds to the magnitude calculated with the refined successive mesh. Successive refinements determine appropriate mesh size until the magnitude of the difference variable for two successive meshes is lower than 1.0%. In this sense, a mesh with 18,271 finite rectangular volumes is used in this study.


**Table 2.** Mesh independence test considering an I-shaped empty channel with φ = 0.05, *H*/*L* = 1.0 and *L*0/*H*<sup>0</sup> = 4.0.

To show the reliability of the present computational model, a verification of the numerical method used here was done by comparing the present numerical solution and a classical analytical solution for the rectilinear case. More precisely, resin advancement as a function of infusion time obtained with the present method and the analytical solution described in Jinlian et al. [50] and Rudd [45] are compared. This comparison can be made for empty channels that lead to rectilinear behavior of resin in the porous medium. For the sake of comparison, a simulation in an I-shaped empty channel with φ = 0.05, *H*/*L* = 1.0, and *L*0/*H*<sup>0</sup> = 19.0, where *L*<sup>0</sup> is equal to *L*, is chosen. The case with these configurations represents a rectilinear case.

An analytical solution for front-line resin advancement as a function of time is given by [50]:

$$X\_f = \sqrt{\frac{2KP\_{\text{int}}t}{\mu\varepsilon}}\tag{13}$$

where *Xf* is the position of the resin front-line (m), *t* is the time (s), μ is the resin viscosity (Pa/s), ε is the porosity, *K* is the permeability (m2), and *Pin* is the injection pressure (Pa). It should be noted that the present formulation is only applicable to a constant *P*<sup>0</sup> situation.

A monitoring line was created in the center of the domain to obtain the numerical results. More precisely, a line that overlaps the *y*-axis was defined by the following points: P1 (*x*<sup>1</sup> = 0.0 m, *y*<sup>1</sup> = 0.0256 m) and P2 (*x*<sup>2</sup> = 0.0 m, *y*<sup>2</sup> = 0.7500 m). The coordinates *x1* = *x*<sup>2</sup> represent the center of the domain in the *x*-direction, while the coordinate *y1* represents the interface placement between the open channel and the porous medium, and *y*<sup>2</sup> is the exit region of the resin. Figure 4 shows the advancement of the resin front line as a function of the injection time for both solutions (analytical and numerical) and an illustration of the simulated domain with monitoring points for resin advancement in the porous medium. A maximum difference of 0.66% was observed between the analytical and numerical results, which is an indication of very good agreement. Therefore, the verification of the numerical model presented was satisfactorily achieved, and the code can be used for theoretical recommendations on the geometric configurations of I and T-shaped empty channels inserted in a rectangular porous plate.

It is also worth mentioning that the present computational model was previously verified and validated in the work by the authors of [9]. The resin propagation obtained with this computational model was verified with the analytical data for linear and radial cases and with numerical results obtained with the PAM-RTM® software for three-dimensional resin flow in a rectangular box and spherical shell. Moreover, the code was validated with the experimental results obtained by Schmidt et al. [51] for a case of resin flow in a cavity mold with an inlet nozzle. More precisely, it was compared to the results of resin flow advancement as a function of time: differences of nearly 8.0% were found. The cases in Reference [9] mimic resin transfer molding (RTM) or light-resin transfer molding (LRTM) processes. In spite of the simulation of different processes, the fluid dynamic behavior of the resin flow of Reference [9] is the same as performed here. These results were reproduced for

code validation in this study, but are not provided here since they have already been presented in Reference [9].

**Figure 4.** Comparison between the placement of resin flow as a function of time by analytical solution.

#### **4. Results and Discussion**

To eliminate residual air from a porous medium, it is necessary to perform successful molding. The formation of voids due to partial impregnation of resin inside the fibrous reinforcement can lead to failures in the use of the final product [45]. In this sense, geometric optimization is essential because it not only increases productivity by reducing the time of injection, but also enhances the quality of the pieces produced. Thus, the prevention of the formation of voids in the porous medium is a restriction in the analysis and definition of the liquid resin infusion (LRI) process.

A reference line was created in the domain to measure filling time (see the horizontal dashed lines in Figure 1a,b, defined by the following points: P1 (*x*<sup>1</sup> =−0.25 m, *y*<sup>1</sup> =0.50 m) andP2 (*x*<sup>2</sup> = 0.25 m, *y*<sup>2</sup> = 0.50 m)). The position of the front line of the resin flow is monitored along this line. When the resin flow completely crosses the reference line, i.e., when the volumetric fraction along the entire line is equal to *f* = 1.0, the infusion process is considered complete.

To evaluate the geometric influence of the I-shaped empty channel over the filling infusion time of the resin along the rectangular porous plate, the fraction area inside the domain is analyzed. Simulations performed with different *L*0/*H*<sup>0</sup> ratios are used to determine the case with lower infusion time for five different fraction areas of the empty channel, φ = 0.005; 0.01; 0.03; 0.05 and 0.1. It has been found that the effect of the ratio *L*0/*H*<sup>0</sup> over filling time (*t*) is similar for all investigated fractions φ. The only exception happened for φ = 0.1 due to the restriction of the geometry in the upper limit of *L*0/*H*0, as can be seen in Figure 5. In general, it is noticed that the lowest *L*0/*H*<sup>0</sup> ratios lead to the best results, i.e., there is a reduction in filling time (*t*) when the channel has higher penetration in the *y*-direction of the domain. However, the progressive increase of *L*0/*H*<sup>0</sup> to the upper limit results in a decrease in injection time too. Thus, there is a globally optimized ratio placed in the lowest region of *L*0/*H*<sup>0</sup> and a local optimized ratio region, placed in the superior limit of the ratio *L*0/*H*<sup>0</sup> investigated. The worst performance, for almost all cases of φ, is reached for intermediate ratios of *L*0/*H*0. Different behavior is noticed only for φ = 0.100, where there is an almost asymptotic growth of filling time as a function of ratio *L*0/*H*0, showing that restriction influences the best shape and effect of geometry over the performance indicator.

The results also showed that the best configurations for lower magnitudes of φ can have a superior performance than those achieved for intermediate ratios of *L*0/*H*0, even for higher magnitudes of φ. For example, the best configuration obtained for φ = 0.01, (*L*0/*H*0)o = 0.015 led to a filling time of *t*<sup>m</sup> = 105.3 s, which is lower than the filling time reached for several configurations (not optimized) with φ = 0.03; 0.05 and 0.1. This behavior indicates that geometric rationalization can lead to smaller

dimension channels for a performance superior to that reached for larger channels, where geometric configuration is not optimized.

**Figure 5.** Effect of the *L*0/*H*<sup>0</sup> ratio over resin infusion filling time for different area fractions of the I-shaped empty channels.

The optimal results obtained in Figure 5 for each fraction area, φ, are compiled in Figure 6a. Figure 6a shows the influence of the ratio between the channel area and the area of the plate (φ) on the once minimized filling time (*t*m) for the I-shaped empty channel. Results in Figure 6a indicate that for the largest area fractions (φ) of the I-shaped channel, the filling time once minimized (*t*m) increases. Figure 6b shows the influence of φ on the optimal ratio of *L*0/*H*<sup>0</sup> named (*L*0/*H*0)o. It is also noted that the behavior of the geometric ratio (*L*0/*H*0)o is modified with the variation of the area occupied by the channel, i.e., there is no universal optimal geometry that leads to the best performance for all area fractions of the open channels.

**Figure 6.** Effect of empty I-shaped channel area fraction (φ) over: (**a**) once minimized filling time (*t*m), (**b**) once optimized *L*0/*H*<sup>0</sup> ratio, (*L*0/*H*0)o.

To illustrate the influence of φ and the *L*0/*H*<sup>0</sup> ratio over resin flow considering the I-shaped channel, Figures 7 and 8 show the behavior of the resin advancement for three different instants of time at φ = 0.01 and 0.1, respectively. The red region represents the resin, i.e., when the volume fraction is *f* = 1.0, while the blue region represents the air (*f* = 0). Regions with different colors represent resin/air mixture with intermediary volume fractions between resin and air (0.0 < *f* < 1.0). It is worth mentioning that this description is applied to all figures where the resin volume fraction is illustrated. Channel configuration is illustrated with the black line.

**Figure 7.** Resin volume fraction for φ = 0.01 and three different ratios of *L*0/*H*0: (*L*0/*H*0)o = 0.015, *L*0/*H*<sup>0</sup> = 2.0, *L*0/*H*<sup>0</sup> = 90.0 as a function of time: (**a**) *t* = 20.0 s, (**b**) *t* = 60.0 s, (**c**) *t* = 100 s.

**Figure 8.** Resin volume fraction for φ = 0.1 and two different ratios of *L*0/*H*0: (*L*0/*H*0)o = 0.11 (optimal), *L*0/*H*<sup>0</sup> = 7.0 (worst scenario) as a function of time: (**a**) *t* = 20.0 s, (**b**) *t* = 40.0 s, (**c**) *t* = 60.0 s.

Figure 7 shows resin advancement for φ = 0.010 and three different magnitudes of the ratio *L*0/*H*0: (*L*0/*H*0)o = 0.015, *L*0/*H*<sup>0</sup> = 1.0 and *L*0/*H*<sup>0</sup> = 90.0, which represents the two extreme magnitudes and an intermediate value of the ratio *L*0/*H*0, where the worst performance was reached. Three different instants of time are illustrated to show resin impregnation along time *t* = 20.0 s (Figure 7a); 60.0 s (Figure 7b) and 100.0 s (Figure 7c). Figure 8 shows the resin flows for φ = 0.100 and two different ratios of *L*0/*H*0: (*L*0/*H*0)o = 0.11 and *L*0/*H*<sup>0</sup> = 7.0, representing the optimal geometry and worst case for three different instants of time: *t* = 20.0 s (Figure 8a), 40.0 s (Figure 8b), and 60.0 s (Figure 8c).

Results of Figures 7 and 8 indicated that for the smallest ratios of *L*0/*H*0, there is a higher advancement of resin in the *y*-direction. After the spread of resin along the entire empty channel, resin flows from the central region of the domain towards the lateral surfaces of the porous plate. Furthermore, in the superior region of the empty channel, the behavior of the resin flow is similar to that found in radial configurations. This behavior makes the distribution of the resin in the fibrous medium more efficient, minimizing resin filling time in the porous medium. It is also noticed that there are no void formations that could delay the infusion process or cause problems in the manufacturing process. For the highest magnitudes of *L*0/*H*0, the channel is placed in the lower region of the porous plate, leading to an almost constant front-line of resin.

Moreover, resin advances towards the exit of the porous domain with a behavior similar to the rectilinear case. In a general sense, geometrical investigation of I-shaped empty channels has attested that the elongated channels in resin propagation direction led to the best performance due to the increase of pressure gradient imposed inside the mold. This behavior is in agreement with the constructal principle of the optimal distribution of imperfections.

Subsequent investigation consisted of the evaluation of the influence of the T-shaped channel aspect ratios *L*1/*H*1, *L*0/*H*0, and *H*0/*H*<sup>1</sup> on the filling time for resin impregnation in the rectangular porous plate. For this configuration, the optimization process was performed only for φ = 0.05 due to the high number of simulations needed. Firstly, the effect of the ratio *L*1/*H*<sup>1</sup> on resin infusion time along the mold was evaluated for different values of the ratio *L*0/*H*0, keeping *H*0/*H*<sup>1</sup> fixed. Figure 9a–d present the effect of the ratio *L*1/*H*<sup>1</sup> over the filling time for various constant ratios of *L*0/*H*<sup>0</sup> and four different magnitudes of *H*0/*H*<sup>1</sup> = 10.0, 20.0, 30.0 and 40.0, respectively. It should be noted that for *H*0/*H*<sup>1</sup> = 30.0 and *H*0/*H*<sup>1</sup> = 40.0, the *L*0/*H*<sup>0</sup> search space becomes more restricted due to the generation of very thin channels. For example, when *H*0/*H*<sup>1</sup> = 40.0, the variation is viable only in the range 0.05 ≤ *L*0/*H*<sup>0</sup> ≤ 0.1 because for higher magnitudes of *L*1/*H*1, the thickness of the bifurcated channel would become insignificant, which disfigures the T-shaped geometry of the channel and leads to the generation of permanent voids.

**Figure 9.** Effect of ratio *L*1/*H*<sup>1</sup> over resin infusion filling time for several values of ratio *L*0/*H*<sup>0</sup> and fixed values of *H*0/*H*1: (**a**) *H*0/*H*<sup>1</sup> = 10.0, (**b**) *H*0/*H*<sup>1</sup> = 20.0, (**c**) *H*0/*H*<sup>1</sup> = 30.0, (**d**) *H*0/*H*<sup>1</sup> = 40.0.

In Figure 9a–d, the continuous curves represent valid geometries, i.e., those geometries with no permanent void inside the porous plate or empty channel. The dashed curves represent cases where permanent voids are generated, even if the monitoring line indicates complete impregnation of the mold (with *f* = 1.0 along the whole measurement line). In this sense, these results are disregarded for optimization evaluation, since void formations are a critical problem in manufacturing. It is to be noted that some cases in the dashed region of the curves achieved lower filling times, indicating that the resin has a higher impregnation in the porous plate. However, the mold was not completely filled.

Results in Figure 9a–d also show that the higher magnitudes of *L*1/*H*<sup>1</sup> for all investigated ratios of *L*0/*H*<sup>0</sup> lead to the best performance, i.e., there is a reduction in resin filling time (*t*) when the bifurcated channel has higher penetration into the mold toward its lateral surfaces. However, some T-shaped channels with higher ratios of *L*1/*H*<sup>1</sup> and lower magnitudes of *L*0/*H*<sup>0</sup> lead to the formation of permanent voids, mainly when the geometry is composed of skinny channels. Therefore, the results recommend the building of elongated channels, since this does not lead to skinny channels. The generation of permanent voids occur, in general, in the corner that connects the single stem to the bifurcated branches. One possible reason for this behavior is the detachment of resin flow in the channel corner with not enough resin flow coming from the inlet to degenerate the void. Therefore, the permanent void remains during all the simulations in the area mentioned above.

To illustrate the influence of ratio *L*1/*H*<sup>1</sup> on mold filling time, Figure 10 shows the advancement of resin for two different ratios of *L*1/*H*1, *L*1/*H*<sup>1</sup> = 6.0 and (*L*1/*H*1)o = 35.0, which represents the worst and best shapes when the constant ratios of *H*0/*H*<sup>1</sup> = 40.0 and *L*0/*H*<sup>0</sup> = 0.1 are assumed. For both cases, the volume fraction fields for three different instants of time—*t* = 10.0 s, 30.0 s, and 60.0 s—are presented in Figure 10a–c, respectively. For *L*1/*H*<sup>1</sup> = 6.0, filling time of *t* = 104.7 s is obtained, while the ratio (*L*1/*H*1)o = 35.0 leads to a resin injection time of the once minimized of *t*<sup>m</sup> = 63.5 s, which is about 65.0% faster than the worst-performing geometry. For the lower magnitudes of *H*1/*L*1, resin distribution occurrs in a radial form in the upper region of the T-shaped channel and linearly from the center to the side surfaces of the mold at the bottom region of the T-shaped channel. This behavior is quite similar to that seen for I-shaped empty channels with high intrusion in the porous plate. As the ratio of *H*1/*L*<sup>1</sup> is changed to (*H*1/*L*1)o = 35.0, the radial advancement of the resin at the upper region of the T-shaped channel is found to be intensified, when compared with the previous case (*H*1/*L*<sup>1</sup> = 6.0). This effect also intensifies the linear advancement of resin in the inferior region of the porous plate.

**Figure 10.** Resin volume fraction for φ = 0.05, *H*/*L* = 1.0, *H*0/*H*<sup>1</sup> = 40.0, *L*0/*H*<sup>0</sup> = 0.1 and two different ratios of *L*1/*H*<sup>1</sup> = 6.0 and (*L*1/*H*1)o = 35.0 as a function of time: (**a**) *t* = 10.0 s, (**b**) *t* = 30.0 s, (**c**) *t* = 60.0 s.

Consequently, the filling of the resin in the horizontal direction (from the open channel towards the side walls) is completed in a shorter amount of time, also leading to an improvement in resin advancement towards the porous plate exit. However, significant increase in bifurcated length (augmentation of the ratio *H*1/*L*1) leads to the formation of permanent voids inside the plate. It is worth mentioning that the channel area is constant; as a consequence, augmentation of the channel bifurcated branch leads to a decrease in channel thickness, which is the main factor responsible for the formation of voids, especially in the corners between the single and bifurcated regions of the channel.

The variation of the *L*1/*H*<sup>1</sup> ratio provides the once minimized filling time (*t*m) of the resin in the porous plate for each investigated ratio of *L*0/*H*<sup>0</sup> and different ratios of *H*0/*H*1. In the second optimization level, the results of Figure 9 are compiled and presented in Figure 11a,b.

**Figure 11.** Effect of the ratio *L*0/*H*<sup>0</sup> for different ratios of *H*0/*H*<sup>1</sup> on (**a**) infusion time (once minimized), (*t*)m; (**b**) ratio (once optimized), (*L*1/*H*1)o.

Figure 11a illustrates the influence of the ratio *L*0/*H*<sup>0</sup> over the once minimized filling time, *t*m, for different *H*0/*H*<sup>1</sup> ratios studied here. It can be noted that for *H*0/*H*<sup>1</sup> = 20.0; 30.0 and 40.0, the effect of the *L*0/*H*<sup>0</sup> over *t*<sup>m</sup> is similar, except for cases with *H*0/*H*<sup>1</sup> = 10.0. In general, it is noticed that lower ratios of *L*0/*H*<sup>0</sup> present a predominance in achieving the twice-minimized filling time (*t*mm), i.e., the best shapes are obtained when the main channel has a higher penetration in the direction of the exit region of the domain. For the ratio *H*0/*H*<sup>1</sup> = 10.0, it can be observed that intermediate ratios of *L*0/*H*<sup>0</sup> lead to the best performance of the problem. For the lowest magnitudes of *L*0/*H*0, the optimal ratios of (*L*1/*H*1)o suffer a strong decrease (see Figure 11b), degenerating the T-shaped channel into an I-shaped one. For higher magnitudes of *L*0/*H*0, the T-shaped channel is strongly restricted in the lower region of the porous plate. Figure 11b highlights the influence of the ratio *L*0/*H*<sup>0</sup> over the once-optimized ratio of *L*1/*H*1, (*L*1/*H*1)o, for the four different magnitudes of *H*0/*H*<sup>1</sup> investigated. The effect of *L*0/*H*<sup>0</sup> over (*L*1/*H*1)o varied for different magnitudes of *H*0/*H*1. For the ratios *H*0/*H*<sup>1</sup> = 20.0; 30.0 and 40.0, the best shapes are obtained for (*L*0/*H*0)o = 0.05 and for high magnitudes of (*L*1/*H*1)oo in the range 20.0 ≤ (*L*1/*H*1)oo ≤ 25.0, corroborating previous findings that the best performance is achieved for the most elongated possible single and bifurcated branches of T-shaped channel. For *H*0/*H*<sup>1</sup> = 10.0, the best configurations changed for (*L*0/*H*0)o = 0.075 and the optimal ratio (*L*1/*H*1)oo dropped dramatically to (*L*1/*H*1)oo = 2.0. In general, the results show that changes in one degree of freedom have a strong influence on the effect of other geometric ratios on the performance indicator. In this case, it was also noticed that there is no optimal universal configuration that leads to the best performance in this problem.

The last optimization step consists of investigation on the influence of the ratio *H*0/*H*<sup>1</sup> over the twice-minimized filling time (*t*mm) and the respective optimal configurations: (*L*0/*H*0)o and (*L*1/*H*1)oo. To obtain this evaluation, the best results achieved in Figure 11a,b are compiled in Figure 12.

These results show that the behavior of the (*t*mm) is largely affected by ratio *H*0/*H*1, and it has a strong dependence on ratios (*L*1/*H*1)oo and (*L*0/*H*0)o. The thrice-minimized filling time (*t*mmm) is obtained for the rate (*H*0/*H*1)o = 30.0, which is only 5.0% lower than that reached for *H*0/*H*<sup>1</sup> = 40.0. In general, it is possible to state that the higher magnitudes of the ratio *H*0/*H*<sup>1</sup> led to the best performance for this problem. The lowest magnitude of *H*0/*H*<sup>1</sup> = 10.0 obtained a performance 183.0% inferior to that reached for the optimal ratio (*H*0/*H*1)o = 30.0. Moreover, differences found for *t*mm between the ratios *H*0/*H*<sup>1</sup> = 20.0 and (*H*0/*H*1)o = 30.0 was nearly 23.0%.

**Figure 12.** Effect of the ratio *H*0/*H*<sup>1</sup> for twice-minimized infusion time, *t*mm, and the respective optimal shapes, (*L*0/*H*0)o, (*L*1/*H*1)oo.

Concerning optimal shapes, it was noticed that the ratio (*L*0/*H*0)o suffered a strong decrease in its magnitude in the region 10.0 ≤ *H*0/*H*<sup>1</sup> ≤ 20.0, followed by a smooth decrease in the range 20.0 ≤ *H*0/*H*<sup>1</sup> ≤ 30.0 and a stabilization for *H*0/*H*<sup>1</sup> ≥ 30.0. For (*L*1/*H*1)oo, an increase in magnitude was noticed, until stabilization for *H*0/*H*<sup>1</sup> ≥ 30.0. In general, the results show the importance of geometric investigation in the reduction of filling time of resin impregnation along a porous domain.

To illustrate the transient behavior of resin flow in the mold, Figures 13 and 14 show the distribution of resin volume fraction as a function of time for the worst performing ratio, *H*0/*H*<sup>1</sup> = 10.0, and for the optimal ratio (*H*0/*H*1)o = 30.0, respectively. Topologies are presented for the following times: *t* = 20.0 s; 40.0 s and 60.0 s, respectively. Figure 7, Figure 8, and Figure 10 illustrate the topology of the configuration of the empty channel.

**Figure 13.** Resin volume fraction for φ = 0.05, *H*/*L* = 1.0, *H*0/*H*<sup>1</sup> = 10.0, for three different ratios of *L*0/*H*<sup>0</sup> = 0.075, (*L*0/*H*0)o = 0.100 and *L*0/*H*<sup>0</sup> = 0.300 as a function of time: (**a**) *t* = 20.0 s, (**b**) *t* = 40.0 s, (**c**) *t* = 60.0 s.

**Figure 14.** Resin volumetric fraction for φ = 0.05, *H*/*L* = 1.0, *H*0/*H*<sup>1</sup> = 30.0, for two different ratios of (*L*0/*H*0)oo = 0.050 and *L*0/*H*<sup>0</sup> = 0.200 as a function of time: (**a**) *t* = 20.0 s, (**b**) *t* = 40.0 s, (**c**) *t* = 60.0 s.

Figure 13 shows resin progression at the analyzed instants of time for *H*0/*H*<sup>1</sup> = 10.0 and *L*0/*H*<sup>0</sup> = 0.075, (*L*0/*H*0)o = 0.1, and *L*0/*H*<sup>0</sup> = 0.3 corresponds to the lowest ratio of *L*0/*H*0, the optimal configuration, and the superior magnitude, respectively. For the lowest ratio of *L*0/*H*0, it is observed that the T-shaped channel almost degenerates into an I-shaped channel. This configuration leads to a poor distribution of resin in the lower region of the porous plate. The resin propagation is shown for three different instants of time: *t* = 20.0 s (Figure 13a), *t* = 40.0 s (Figure 13b) and *t* = 60.0 s (Figure 13c).

For optimal configuration, despite the T-shaped channel being placed in a lower region of the porous plate in comparison with the case with *L*0/*H*<sup>0</sup> = 0.075, the elongated bifurcated branches led to better distribution of the resin in the lower region of the porous plate and allow easy resin propagation towards the exit of the domain. For the superior ratio of *L*0/*H*<sup>0</sup> = 0.3, the restriction of the T-shaped channel in the lower region of the porous plate affects the advancement of resin towards the exit of the domain, leading to a poorer performance than the optimal ratio of *L*0/*H*0.

Figure 14 presents the results of resin spread for the optimum ratio of *H*0/*H*1, (*H*0/*H*1)o = 30.0, for two different ratios of *L*0/*H*0, (*L*0/*H*0)oo =0.05, *L*0/*H*<sup>0</sup> =0.2 and three different instants of time: *t* = 20.0 s (Figure 14a), *t* = 40.0 s (Figure 14b), *t* = 60.0 s (Figure 14c). These ratios represent the corresponding twice-optimized and worst-case scenarios for the ratio of (*H*0/*H*1)o = 30.0, respectively. For the ratio (*L*0/*H*0)oo = 0.05, an injection time of *t*mmm = 27.9 s is obtained, and the ratio *L*0/*H*<sup>0</sup> = 0.2 leads to a once minimized injection time of *t*m = 87.9 s, which is, about of 215.0% higher than the previously obtained optimal geometry. It is then possible to significantly improve the performance of the studied system even in the last level of optimization.

In a general sense, the results demonstrate that for the T-shaped channel, all geometrical ratios (degrees of freedom) affect the filling time of resin impregnation in a porous domain. They also show that one degree of freedom has a strong influence over the optimal ratios of other geometric ratios, e.g., the ratio *H*0/*H*<sup>1</sup> affected the optimal ratios of (*L*1/*H*1)oo and (*L*0/*H*0)o.

All the geometrical parameters studied here (*L*1/*H*1, *L*0/*H*0, and *H*0/*H*1) has some sensibility in the infusion process and cannot be neglected. Moreover, a comparison between the performance of the T-shaped and I-shaped empty channels, for the same conditions, demonstrates that the thrice-optimized T-shaped channel leads to performance that is almost three times superior to that achieved by the optimal I-shaped configuration.

In the present fluid dynamic conditions, where resin flow has little advective dominance, the rectangular porous plate can be viewed as an elementary construction of multiple structures generated by the replication of a porous plate in *x-* or *y*-axes. In this case, multiple optimal configurations of I and T-shaped channels can be associated in a series and/or parallel using the optimal configurations reached in this work.

#### **5. Conclusions**

This paper developed a numerical study of the geometric evaluation of empty channels inserted in a porous rectangular plate, mimicking the liquid resin infusion process (LRI), and applying constructal design and exhaustive search (ES) for geometric optimization. The primary purposes here were to analyze the influence of the design of I and T-shaped empty channels on the filling time of resin impregnation in a porous plate, preventing the generation of permanent voids.

The following important recommendations were made:


Other geometrical configurations, such as Y-shaped and fishbone empty channels, or geometries being constructed with no pre-defined form are recommended for future studies. Another interesting possibility of study is the insertion of two or more inlet ports investigating the influence of interactions between different channels on their optimal shapes.

**Author Contributions:** Conceptualization, J.A.S., and E.D.d.S.; methodology, G.M.C.M., L.A.I., and R.d.L.L.; software, G.M.C.M., J.A.S., and E.D.d.S.; validation, L.A.I., C.F., G.M.C.M., and J.A.S.; formal analysis, E.D.d.S., L.A.O.R., J.A.S., and C.F.; investigation, G.M.C.M., E.D.d.S., and J.A.S.; resources, S.C.A., J.A.S., and E.D.d.S.; data curation, E.D.d.S., and J.A.S.; writing—original draft preparation, G.M.C.M., R.d.L.L., and E.D.d.S.; writing—review and editing, E.D.d.S., L.A.I., and C.F.; visualization, J.S. and L.A.O.R.; supervision, G.M.C.M., L.A.O.R., and C.F.; project administration, E.D.d.S. and J.A.S.; funding acquisition, G.M.C.M., J.A.S., and E.D.d.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Brazilian Coordination for the Improvement of Higher Education Personnel (CAPES) (Finance Code 001), the Brazilian National Council for Scientific and Technological Development (CNPq), and the Italian Minister of Foreign Affairs and International Cooperation (MAECI), as part of the 'Two Seats for a Solar Car' international project.

**Acknowledgments:** The authors thank Giangiacomo Minak for supporting this publication.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Nomenclature**


#### **Greek symbols**


#### **Subscripts**


#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Material Characterization for Reliable Resin Transfer Molding Process Simulation**

**Maria Pia Falaschetti 1,\*, Francesco Rondina 1, Nicola Zavatta 1, Lisa Gragnani 2, Martina Gironi 1, Enrico Troiani 1,\* and Lorenzo Donati <sup>1</sup>**


Received: 11 February 2020; Accepted: 3 March 2020; Published: 6 March 2020

**Abstract:** Resin transfer molding (RTM) technologies are widely used in automotive, marine, and aerospace applications. The need to evaluate the impact of design and production critical choices, also in terms of final costs, leads to the wider use of numerical simulation in the preliminary phase of component development. The main issue for accurate RTM analysis is the reliable characterization of the involved materials. The aim of this paper is to present a validated methodology for material characterization to be implemented and introduce data elaboration in the ESI PAM-RTM software. Experimental campaigns for reinforcement permeabilities and resin viscosity measurement are presented and discussed. Finally, the obtained data are implemented in the software and then compared to experimental results in order to validate the described methodology.

**Keywords:** RTM; composites; FEM simulation; permeability characterization

#### **1. Introduction**

Fiber-reinforced composite materials are widely used in several lightweight applications in the nautical, aerospace, and automotive sectors due to high strength-to-density ratios. More extensive applications are limited by materials and technology costs. Even if the cost of carbon and glass fibers follows a decreasing trend, the manufacturing costs will still be almost constant if conventional high-performance technologies are used (i.e., autoclave manufacturing). Moreover, several composite technologies cannot provide a component cost reduction by increasing the production volume. To achieve a reliable component cost reduction, new out-of-autoclave technologies (like infusion technologies, prepreg compression technologies, or short-fiber compound processing) have been recently optimized for industrial use.

Among the others, resin transfer molding (also known as light RTM), and its variants vacuum-assisted resin infusion (VARI) and high-pressure resin transfer molding (HP-RTM), are promising technologies for producing complex components with high performance at a higher volume rate. Moreover, these processes are economically advantageous with respect to the others, such as autoclave curing of pre-impregnated layers, due to cheaper materials and shorter processing times.

RTM technologies consist of placing a dry three-dimensional preform into a mold cavity and injecting resin-hardener mixture into the closed mold. The injection pressure (assisted by a vacuum for VARI or high pressure in HP-RTM processing), together with fiber permeability, is able to produce components with 50% and higher fiber volume fraction, thus providing excellent mechanical properties [1–3]. Nevertheless, in the light RTM process, the low flow rate of resin and low injection pressure result in a long resin injection time, thus hampering the use of fast-cure resins and, consequently, the overall productivity. As a consequence, the light RTM process is mainly limited to the low-volume manufacturing capability and, only to a minor extent, to slightly lower mechanical properties related to the lower fiber volume fraction. High-volume manufacturing with RTM is possible only if the process cycle time is significantly reduced. VARI and HP-RTM variants are improvements in this direction. Due to the differential pressure between inlet and outlet flows (VARI) or to higher flow rates (HP-RTM), together with the choice of adequate resin and hardener, production time can be reduced. A critical point of RTM technologies is the assurance of a full impregnation of the component, which is influenced by reinforcement types, the number of layers, resin viscosity, position of resin injection (ports), and vacuum extraction (vents). As molds can reach very high manufacturing costs when wide components are involved or when metal molds are required, tools supporting the design and critical process selection are needed in order to reduce uncertainties already in the preliminary phase of component development.

In this direction, Finite Element Method (FEM) software has been developed in the last few years to assist in component design and processing, predicting issues such as resin-rich areas, air bubbles, dry spots, zones of high porosity, etc., and for optimization of mold geometry and process parameters, such as port and vent locations [4–7]. In the case of critical molds, therefore, it is advantageous to perform a preliminary study by means of numerical simulation in order to avoid expensive experimental setups and prototypes.

In this context, the ESI group developed a PAM-RTM module, which can solve liquid resin infusion processes, like light and high-pressure RTM, compression resin transfer molding (C-RTM), and VARI. ESI's approach, resulting from several years of collaboration with academics [8–10], consists of a finite element solver based on the coupling of resin flow (governed by Darcy's law) and the preform behavior (considered as a porous medium undergoing deformations). Thus, the solver can provide the filling time and properties (thicknesses, fiber volume contents, geometry) of the final product [11].

However, accurate information on the material properties is needed for reliable simulations of the RTM process: in particular, reinforcement permeabilities and resin viscosity. Several works in the literature are available for the experimental determination of the permeability of RTM reinforcements [12]. Most of them evaluate the in-plane permeability [13,14], concluding that radial flow (liquid injected from transverse direction at well-defined locations to evaluate 2D permeability) and linear flow (liquid injected from one end of the composite preform to evaluate 1D permeability) experiments give consistent results. The through-the-thickness permeability is more complex to determine [15–21] and is usually neglected because of the small thicknesses of the preforms.

Therefore, a well-defined procedure for the experimental determination of the permeability tensor of fiber textiles and resin viscosity has been formulated in this paper in order to provide a reliable method to simulate RTM production of industrial components using PAM-RTM software. The key steps, from material characterization to numerical simulations, are illustrated by referring to the production of a real component in order to compare to a real-life test case. The permeabilities of four different types of textile reinforcements (three carbon fibers and one glass fiber) were experimentally measured. Moreover, the viscosity of a synthetic oil (15W40, used in glass textile permeability tests) and orthophthalic resin (Lavesan LERPOL 666/S RAL 9010) was tested. Permeability tests were modelled in the case of carbon fiber reinforcement tests, while in the case of glass fibers (GFRP), a naval component was reproduced. The analyzed GFRP part represents coverage for an aft-peak, placed aft-ward the cockpit, and it belongs to a series of pieces that are mounted on a 15.85 m (52 ft) cruising sailing yacht. This choice was sustained by the wider experience and higher amount of available data that RTM technology shows in marine industry compared to automotive, where production information represents more sensible data.

#### **2. Materials and Methods**

#### *2.1. Test Material*

The naval component was made of a polyester resin reinforced by glass fibers. In particular, the following materials were used:


The nominal properties of the resin are reported in Table 1. Accurate characterization of the flow properties of the resin is vital for manufacturing control, as well as for numerical modelling of the process. To this end, experimental tests were conducted to assess the resin viscosity and the fiber volume fraction and permeability of the textile.



Moreover, in order to provide references respect to other fibres and styles, carbon fiber reinforcement textiles were tested in permeability tests. These materials are a unidirectional fiber 12K T700 300 g/m2, a triaxial fiber (30/90/−30) 24K T700 300 g/m2, and a Twill 2×2 fabric 12K T700 630 g/m2.

#### *2.2. Methods Used for Material Properties Characterization*

Resin viscosity was measured by using a Thermo Scientific Brookfield Haake 7 plus viscosimeter. Four values of the rotational speed were used for testing, i.e., 20, 30, 50, and 100 rpm, considering both room and production temperature. The viscosity of the oil employed for the permeability tests was also measured by using the same apparatus.

The textile permeability was assessed by means of unsaturated linear flow tests. Assuming that the textile reinforcement in an RTM component is a porous medium through which the fluid flows, the resin flow can be estimated by Darcy's law:

$$V = -\frac{K}{\mu} \nabla P\_{\prime} \tag{1}$$

where *V* is the velocity of the flowing resin, ∇*P* is the applied pressure gradient, μ is the dynamic viscosity of the resin, and *K* is the permeability tensor. Using a coordinate system aligned to the principal directions of the fiber reinforcement, the permeability tensor can be expressed in the following diagonal form:

$$K = \begin{bmatrix} K\_1 & 0 & 0 \\ 0 & K\_2 & 0 \\ 0 & 0 & K\_3 \end{bmatrix}. \tag{2}$$

Thus, measurements of only three distinct scalar coefficients is required, which can be done according to existing procedures reported in the literature [13–21]. A further simplification is possible: as the thickness of the reinforcement is small compared to its length and width, the problem can be assumed to be 2D. This leaves only the two planar components of permeability, i.e., *K*<sup>1</sup> and *K*2, to be determined. The principal directions of the reinforcement, with respect to which the components of permeability are found, are also found from the experimental tests.

By applying a given pressure gradient, an unsaturated fluid flow is obtained. In the case of a mono-dimensional, i.e., linear, unsaturated flow, Darcy's equation can be expressed as

$$\frac{d\mathbf{x}}{dt} = \frac{\mathbf{K}\_{\text{xx}}\Delta P}{\mu \varepsilon \mathbf{x}(t)},\tag{3}$$

where *x*(*t*) is the position of the flow front at time *t*, Δ*P* is the pressure difference between the inlet and outlet, ε = 1 − *Vf* is the void ratio (portion of volume not occupied by the fibres), and *Kxx* is the component of the permeability matrix along the direction of the flow. With the initial conditions *x*(*t* = 0) = 0, the solution of Equation (3) reads:

$$\frac{\mu\_i^2}{2} = \frac{K\_{xx}\Delta Pt\_i}{\mu\varepsilon},\tag{4}$$

which allows one to measure the permeability *Kxx*, provided the position of the flow front is tracked throughout the test.

In order to identify the principal directions and the corresponding components of the permeability matrix using linear injection experiments, tests in three planar directions were performed, i.e., 0, 45, and 90◦ orientations. The principal directions were then found by [22,23]

$$K\_1 = K\_{\exp}^0 \frac{\alpha\_1 - \alpha\_2}{\alpha\_1 - \frac{\alpha\_2}{\cos(2\beta)}},\tag{5}$$

$$K\_2 = K\_{\exp}^{\otimes 0} \frac{\alpha\_1 + \alpha\_2}{\alpha\_1 + \frac{\alpha\_2}{\cos(2\beta)}},\tag{6}$$

$$\beta = 0.5 \tan^{-1} \left( \frac{\alpha\_1}{\alpha\_2} - \frac{\alpha\_1^2 - \alpha\_2^2}{\alpha\_2 K\_{\text{exp}}^{45}} \right) \tag{7}$$

where the parameters α<sup>1</sup> and α<sup>2</sup> are defined as

$$\alpha\_1 = \frac{K\_{\rm exp}^0 + K\_{\rm exp}^{\rm 0}}{2},\tag{8}$$

$$a\_2 = \frac{K\_{\text{exp}}^0 - K\_{\text{exp}}^{90}}{2}. \tag{9}$$

#### *2.3. Experimental Setup*

The permeability tests were performed on a 1D linear flow test bench. Two different molds were used to test glass and carbon reinforcements (test bench used for glass fibers is shown in Figure 1).

**Figure 1.** Linear flow test bench used for measurement of glass textile permeability. The oil flow front is clearly visible in the detailed view.

In order to obtain accurate and reproducible results, permeability tests should be carried out with non-reactive fluids: synthetic oil 15W40 and water were used instead of the resin to measure the textile permeability of glass fiber and carbon fiber, respectively. The oil was selected for its controlled viscosity profile over a broad range of temperature, while water is an option in the case of shortage of textiles to be tested (giving the possibility of reuse). The mold cavity had a size of 440 × 100 × 6 mm [16,23–26]. A free length of 20 mm was left at the inlet and outlet on the two sides of the specimen, in order to initiate a stable flow front and prevent potentially dangerous migrations of the fluid to the vacuum apparatus. Thus, a buffer volume was created, which works additionally as a device to reduce pressure fluctuations at the vent. A Plexiglass cover was positioned on top of the mold and fixed using a set of screws; a sealing rubber assured contact between the two mold halves and avoided leakages. Inlet and outlet ports were connected to the top plate for ease of construction.

The mold cavity thickness was chosen in accordance with literature investigations [23,24] and should be considered as an optimum that minimizes two adverse effects. On the one hand, capillarity at the boundaries becomes significant in shallower cavities, resulting in an incorrect reading of the fluid flow; on the other hand, the thickness effects cannot be neglected if thicker cavities are used, thus invalidating the assumption of planar flow.

An adequate stack of textiles should be used to fill the cavity and realize the same fiber volume content as that on the part to be manufactured. Care must be taken when cutting the textiles to avoid wrinkles or voids at the edges of the mold cavity, as these would affect the planar fluid flow. Three repetitions were taken for each direction, i.e., 0◦, 45◦, and 90◦ orientations, for a total of 9 trials, for each material. A camera was used to monitor the injection process, and time was measured with a chronometer.

#### *2.4. Numerical Model*

A numerical model was developed to simulate the resin infusion process. This consists of modelling the flow of the resin through a porous medium, i.e., the textile. Neglecting the capillarity forces of attraction or repulsion acting at the flow front (as they are deemed sufficiently small in front of the pressure field in RTM), the flow is governed by Darcy's law, e.g., Equation (1). Thus, considering the resin as an incompressible fluid and combining it with Darcy's law result in Richard's equation:

$$\nabla \cdot \left(\frac{1}{\mu} K \nabla P\right) = 0.\tag{10}$$

Equation (10) was solved in the software PAM-RTM using the non-conforming finite element method. Two numerical models were developed: The first one simulates the permeability tests of the carbon textile reinforcements, and the second one simulates the infusion process of the entire fiberglass component. The model of the textile permeability test was made up of 52652 3D tetrahedral elements.

The pressure difference applied between the inlet and outlet was equal to 40 kPa.

The fiberglass component was modelled with a 2D model in order to diminish the computational cost. A schematic of the model is shown in Figure 2. Four distinct regions can be identified: areas A and B are made of glass fiber textile layers, with one and two textile layers, respectively; region C consists of an empty rectangular channel that is filled by the resin at the beginning of the infusion process; region D is the foam core area. The injection valve is positioned in region C, while the vacuum valve is placed in the middle of D.

**Figure 2.** Numerical model description.

The material properties input into the model are given in Table 2. A total of 223,927 triangular elements were used to mesh the body; the core included 862 holes and required a fine mesh with an element size of approximately 2 mm. The empty channel permeability value (1.25 <sup>×</sup> <sup>10</sup>–7 <sup>m</sup>2) was used for the permeability of the holes, neglecting the contribution of the glass fibers filling. A sensitivity study was run to substantiate this assumption: by varying the hole permeability between 10−<sup>3</sup> and 10−<sup>7</sup> m2, the filling time changed by less than 2%. This shows that the permeability of the holes has little influence on the filling time of the entire component.


**Table 2.** Material data of the fiberglass component model.

A pressure of 101.3 kPa was applied at the injection valve, while the pressure at the vent, i.e., in the middle of region D, was set equal to 5 kPa.

#### **3. Results and Discussion**

#### *3.1. Material Characterization Tests*

The values of the resin and test oil viscosity measured in the experimental tests are given in Table 3. Clearly, as the room and production temperatures are almost equal, the resin viscosity is almost unaffected by these temperature changes.


**Table 3.** Resin and test oil viscosity as measured in experimental tests.

Table 4 shows the results of the permeability tests, grouped according to their orientation. As visible from the table, the scatter in the data is limited: for every orientation, the results are within a 5% error range.

**Table 4.** Glass textile permeability measured in the tests.


By using Equations (5)–(7), the principal textile permeabilities can be computed. This gives the values collected in Table 5.


Knowing the values of the principal permeabilities and their corresponding orientation, it is possible to compute the effective permeability *Keff* along a general direction θ, which is defined as

$$\sqrt{K\_{eff}} = \frac{\sqrt{K\_1 K\_2}}{K\_1 \sin(\theta) K\_2 \cos(\theta)}.\tag{11}$$

This allows one to visualize the flow distribution in the given layup as a function of the orientation angle θ. For the tested glass textile layup, the resulting effective permeability is plotted in Figure 3. Effective permeability of the fiberglass textile plotted as a function of the orientation angle θ.

**Figure 3.** Effective permeability of the fiberglass textile plotted as a function of the orientation angle θ.

As for the carbon fiber reinforcements, the permeability of the three different carbon fiber reinforcements was tested, i.e., (ply A) unidirectional fiber 12K T700 300 g/m2, (ply B) triaxial fiber (30/90/−30) 24K T700 300 g/m2, and (ply C) Twill 2 <sup>×</sup> 2 fabric 12K T700 630 g/m2.

Three tests were performed for each textile. Data obtained are listed in Tables 6 and 7:


**Table 6.** Local permeability of the carbon textile.



Data show that the permeabilities of the selected carbon fiber textile styles are at least one order of magnitude lower than those of glass fiber. The achievement was expected as the glass fiber textile comprises a multi-layer of 600 g/m<sup>2</sup> E-glass biaxial fibers [+/−45◦ orientation] and a 450 g/m<sup>2</sup> chopped strand mat stitched by means of a 18 g/m<sup>2</sup> PP for enhancing permeability characteristics. In addition, it is shown that anisotropic permeability occurs even in balanced textiles where no apparent preferential direction is expected. Under deeper analysis, the ratio between major and minor in-plane permeability spans between 1.2 for the glass textile and 2.0 for carbon triaxial fiber (ply B) and 1.3 for the carbon 2 × 2 twill. These results indicate that proper characterization of the textile permeability is crucial in the design phase of an infusion process, as the filling path and time can be significantly affected by the fabric configuration and layup.

#### *3.2. Numerical Models and Experimental Comparison*

The average filling time measured in the carbon textiles permeability tests is reported in Table 8, together with values computed by the numerical model. Comparison between the experimental and numerical results shows a limited difference within experimental data scatter. For better visualization of the numerical data, they are also plotted in Figure 4**.** as a contour map. Clearly, the points located further from the injection valve, i.e., the right-most point in the figure, have the highest filling time.


**Table 8.** Carbon textile reinforcement permeability experimental and numerical filling time results.

**Figure 4.** Contour map of the filling time computed numerically: example of case Ply C 0◦ (color scale in s).

For the fiberglass component, the filling time calculated in the simulation is equal to 304.3 s and the resin flow is uniform in the entire component. Figure 5 shows the position of the flow front at four time instants, while a contour plot of the filling time is given in Figure 6.

As visible in Figure 5, the resin flows through all the empty borders in less than 1 s. After that, the flow front uniformly follows the vertical faces and the flow eventually forms a circular shape around the vacuum valve.

Therefore, the analysis of the flow pattern in the numerical simulation shows that the injection and vent valves are in good positions to assure an adequate filling path of this simple component.

The filling time of the real component was measured to be 303 s. A comparison of this value to that predicted by the numerical model shows an error of 0.44%.

Moreover, by analyzing the component pulled out of the mold after resin polymerization (Figure 7), there was an absence of holes, bubbles, or other defects. This shows that, also in the real condition, valve locations and pressure settings are adequate in order to obtain a good structure.

Therefore, comparing the numerical and experimental results shows a good correspondence in all analyzed cases. This highlights the reliability of the methodology discussed in this paper. In fact, it has been proved that the RTM process simulation, by means of the PAM-RTM software, could give reliable results for optimization of the injection strategy in complex components, provided that realistic values are input to the model. To this end, accurate material characterization is needed, which must take into account the dependence of the material properties on the specific experimental conditions (e.g., temperature). Moreover, the results presented are in good accordance with other studies where good correspondence between RTM (as well as VARI and C-RTM) modelling and experimental results was obtained [1–3,27,28]. In summary, a workflow explaining the phases and procedures for the characterization of the materials and then to simulate the RTM process, as followed in this paper, is reported in Figure 8.

**Figure 5.** Visualization of the flow front calculated by the numerical model.

**Figure 6.** Contour plot of the filling time computed by the numerical model.

**Figure 7.** Fiberglass component manufacturing: (**a**) Dry material before resin infusion, (**b**) finished part.

**Figure 8.** Workflow to model the resin transfer molding (RTM) process of a real component.

#### **4. Conclusions**

In this paper, the procedure for implementing reliable inputs in the simulation of an RTM-based process is presented and discussed. The experimental setup for textiles permeability characterization is described, as well as the procedure for retrieving local and global permeability. Four different types of textiles (three carbon fiber styles and one glass fiber multilayer textile) have been characterized in order to provide to readers reference data for different types and styles of textiles. The permeabilities of the carbon fiber textiles were found to be one order of magnitude lower than that of glass fiber, whereas fiber alignment had a less relevant influence. Moreover, tests show that anisotropic permeability occurs even in balanced textiles where no apparent preferential direction is expected. The characterization of the viscosity of the orthophthalic resin Lavesan LERPOL 666/S and of the synthetic oil 15W40 showed comparable values at room temperature.

ESI PAM-RTM software was used to perform RTM simulations. As a first step, simulation of permeability characterization in the different directions was performed by means of 3D models. Filling time results showed an average good correspondence with experimental data (3.6% average error). As a second step, the resin infiltration of a 15.85 m (52 ft) sailing yacht component, comprised of stitched glass fibers around a perforated core, was implemented in a 2D model. Output data including filling flow and time proved consistent (0.44% error in filling time) with data acquired in industrial environment, validating the procedure proposed in this paper.

**Author Contributions:** Conceptualization, L.D. and F.R.; methodology, L.D.; software, M.P.F., L.G. and M.G.; validation, M.P.F., E.T. and N.Z.; formal analysis, M.P.F. and F.R.; investigation, F.R.; data curation, L.D.; writing—original draft preparation, M.P.F.; writing—review and editing, N.Z and E.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors thank Corset&Co Srl for the naval component production.

**Conflicts of Interest:** The authors declare no conflict of interest

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Investigation of Electrodeposition External Conditions on Morphology and Texture of Ni**/**SiC***w* **Composite Coatings**

#### **Liyan Lai, Hongfang Li, Yunna Sun, Guifu Ding \*, Hong Wang and Zhuoqing Yang**

National Key Laboratory of Science and Technology on Micro/Nano Fabrication, School of Electronic Information and Electrical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China **\*** Correspondence: gfding@sjtu.edu.cn; Tel.: +86-21-34206686

Received: 19 August 2019; Accepted: 9 September 2019; Published: 12 September 2019

**Abstract:** Recently, a strategy of synthesizing SiC whisker-reinforced nickel (Ni/SiC*w*) composites with excellent mechanical properties by electrodeposition has been proposed for exploring its potential applications in micromechanical devices. In this paper, a series of external conditions that affected the content of SiC whiskers in composite films were studied, such as cathode current density, stirring rate and electrolyte temperature. The experimental results indicated that the optimum morphology was obtained at a stirring speed of 300 rpm, a temperature of 50 ◦C, and a current density of 18 mA/cm2. Additionally, the content of SiC whiskers and textural preference were also investigated by varying its external conditions, and the results demonstrated that the composites with high mass percentage whiskers are more advantageous for electrocrystallization of Ni in the (200) orientation. Finally, the relationship between external conditions and intrinsic morphology, composition and texture of Ni/SiC*w* composites was revealed, and it provides a constructive approach to fabricate the high-content SiC whiskers of these composites.

**Keywords:** physical conditions; electrodeposition; SiC whisker; texture; morphology

#### **1. Introduction**

Electrodeposition is one of the most commercially successful and inexpensive superior techniques for fabrication of metallic coatings. However, the conventional electrodeposition methods are very simple, since it has only one phase throughout the process. The pure electroplated metals are difficult to satisfy the current requirements, due to their plastic deformation. In order to resolve this problem, the method of co-deposition insoluble solid particles in metal matrix was established and recognized for fabrication the composite coatings. The reinforcing particles, suspending in the electrolyte, are entrapped and incorporated into the metal matrix by adopting suitable methods during the co-deposition process, such as electrophoresis, adsorption or mechanical entrapment [1]. In 1928, Fink et al. [2] introduced the co-deposition method to produce Cu/graphite composite coatings that applied in car engines for the first time. Besides, the first patent of co-deposition was issued by Grazen, in 1962 [3]. Since then, considerable researches have been focused on different types of composites fabricated by the co-deposition method. These composites include ceramics, metals, polymer or microcapsule/liquid incorporated into metal matrices, which could improve the coating properties as mechanical strength, wear/friction-resistance, high-temperature oxidation-resistant and corrosion-resistance [4–15]. Based on the SiCw with features in high modulus and tensile strength, good wear resistance, as well as dimensional stability, recently, we have successfully synthesized high mechanical strength Ni/SiCw composites by electrodeposition.

Since the performance of composite coatings is highly sensitive to external conditions during the co-deposition process, such as current density, temperature, stirring speed, and current density, etc. Therefore, numerous efforts have been dedicated to improving its microstructure and micromorphology of the deposited coatings for enhancing their properties. Lately, relevant studies have revealed that the content of particles in the composite coating can be controlled by the external conditions during the deposition process, and resulting in a difference in the performance of the composite coatings. H. Gül et al. [16] reported the nickel metal matrix composites reinforced with SiC submicron particles by electrodeposition. Besides, the influence of stirring speed and particles concentration on co-deposited composite coatings was investigated. The results indicated that the maximum of particle concentration was 20 g/L and the high stirring speed was chosen to obtain the wear resistant nickel matrix coatings. P. Gyftou et al. [17] applied both direct and pulse current conditions to produce Ni/ nano-SiC composites, and it proved that pulse electrodeposition significantly improved the hardness of the Ni/SiC composite deposits. C.K. Chung et al. [18] investigated the electrodeposition process of Ni at low electrolytic temperatures (5–20 ◦C). It demonstrated that low temperature could enhance the hardness up to 6.18 GPa and produce the smoothness of films. Tushar Borkar et al. [19] presented the effects of the deposition conditions and the nanoparticle concentration in the electrolyte on the surface microstructure, crystallographic micro-texture, microhardness, and tribological properties of coatings. The reinforcement of Ni-Al2O3 nanoparticles significantly improved microhardness and wear resistance of the composite coatings. Based on the previous studies, it is necessary to clearly understand the influence of various conditions on the co-deposition of SiC whiskers in the nickel matrix, so as to consciously control the whiskers content and effectively manipulate the performance of the composite coatings. The purpose of this study is to achieve high SiC whiskers content in Ni/SiC*w* composite coatings. In addition, the influence of the parameters, such as stirring speed, electrolyte temperature and cathode current density on the SiC*w* content in composites were investigated, which clearly indicated the relationship between the external conditions and the content of SiC whiskers in the composites.

#### **2. Materials and Preparation**

The SiC whiskers (β- SiC*w*, XFNANO Material Co., Ltd., Nanjing, China), with an average length of 50 μ*m* and diameter of 0.2 μ*m*, were selected. Ni/SiC*w* composite coatings were prepared by a constant-current electrodeposition method from a nickel sulfamate electrolyte containing SiC whiskers as a reinforcement phase. Prior to the SiC*w* addition into the electroplating solution, it was first stirred in 3 vol. % hydrofluoric acid and refluxed for 8 h at 95 ◦C [20], then ultrasonically washed in distilled water until the pH becomes neutral, subsequently modified by γ-aminopropyltriethoxysilane (KH550). These treatments can remove the SiC*w* surface impurities (SiO2) and improve its wettability, as well as dispersibility. After the above treatments, SiC whiskers were added in the electrolyte by continuous magnetically stirring with a rate of 300 rpm for at least 24 h, which can prevent whiskers agglomeration and maintain them in suspension state. According to a number of optimization experiments, the optimal whisker concentration in the electrolyte was 0.8 g/L. The electrolyte composition and electrodeposition parameters are listed in Table 1. Analytical reagents and deionized water were used to prepare the plating solution. After electrodeposition, Ni/SiC*w* composite coatings were ultrasonically cleaned in distilled water for 10 min to remove loosely adsorbed SiC*w* from the surface.

The surface morphology and microstructure of the composite coatings were characterized by scanning electron microscopy (SEM; ULTRA55, Zeiss, Germany). The amount of embedded SiC*w* was evaluated by the energy dispersive X-ray spectroscopy (EDS) with the same system of SEM. The phase structure analysis of the coatings was conducted by X-Ray Diffraction (XRD, D8 Advance, Bruker-axs) operating with Cu Kα (λ = 1.54178 Å) radiation at room temperature. The preferred XRD orientation index TC(*hkl*) was calculated, and the texture co-efficient (TC) for each (hkl) reflection is given by [21,22]:

$$TC(hkl) = \frac{I\_{(hkl)} / I\_{0(hkl)}}{\frac{1}{n} \sum I\_{(hkl)} / I\_{0(hkl)}} \tag{1}$$

where *TC* (*hkl*) is the texture coefficient of the specific (*hkl*) plane; *I* (*hkl*) and *I*<sup>0</sup> (*hkl*) are the diffraction intensity of the electrodeposited nickel and the powder pattern in the JCPDS cards, respectively; n is the number of reflection peaks used in the calculation. Three peaks of (111), (200) and (221) are used to calculate the texture coefficient.



#### **3. Results and Discussions**

*3.1. E*ff*ects of Electrodeposition Parameters on Morphology and Composition of the Ni*/ *SiCw Composite*

Figure 1 shows the surface morphology of various current density for the Ni/SiC*w* composite coatings, which was prepared at a stirring speed of 300 rpm, a temperature of 50 ◦C, and a current density from 2 mA/cm<sup>2</sup> to 40 mA/cm2. As shown in Figure 1, obviously, the microstructure and content of SiC*w* in coatings were appreciably affected by the current density. When the current density was lower than 18 mA/cm2, the glossy surface of the nickel matrix with uniformly distributed SiC whiskers could be seen, and the content of SiC*w* increased with the increase of current density (Figure 1a–d). However, when the current density was higher than 25 mA/cm2, the whiskers content reduced, and the surface morphology of coatings became nodular and inhomogeneous (Figure 1e–f). These phenomena made clear that the electrodeposits prepared at low current density tend to form smooth and luminous surface morphology with fine grains.

**Figure 1.** SEM images of Ni/SiC*w* composite coating fabricated from the electrolyte at stirring speed of 200 rpm, and temperatures of 35 ◦C, and current density of (**a**) 2 mA/cm2, (**b**) 6 mA/cm2, (**c**) 10 mA/cm2, (**d**) 18 mA/cm2, (**e**) 25 mA/cm2, (**f**) 40 mA/cm2.

The influences of current density on composite coatings could be explained by the following aspects. Firstly, as the cathode current density increases, the Coulomb force between the whisker-adsorbed metal ions and the cathode increases, which can increase the deposition rate of the matrix metal and save time for SiCw to be embedded in Ni matrix. That is, the content of SiCw in the nickel matrix increases at per unit time. Secondly, as the cathode current density increases, the cathode overpotential and the electric field force increases correspondingly, which facilitates the electrostatic attraction of the cathode to the solid particles adsorbing the positive ions. At higher overvoltage, the transportation of metal ions to cathode becomes an important factor in electrodeposition, and it is well known that

the addition of inert particles could enhance this transportation [23]. Compared with lower current density, the rate of whisker that transfer to cathode surface and embed into the composite coatings increases remarkably when the cathode current density is too high. Thirdly, part of the surface of the cathode is covered by whiskers, due to its embedding behavior, which results in the increase of the cathode overpotential, the enhancement of hydrogen adsorption capacity and the formation of alkaline salts near the cathode [4]. These by-products that adsorbed on the cathode surface or embedded into the coating also prevent the whisker from co-depositing with nickel metal and lead to a decrease of the SiCw content in electrodeposits. Besides, at a higher overpotential, the dendrites or nodules were formed on the coating's surface, due to the discharge of metal ions near the cathode, which cause the discharge at the swelling of deposits and it is consistent with previous literature reported by other researchers [2].

During the electrodeposition process, the deposited crystal and grew behavior are closely related to the electrolyte temperature, which plays an essential role in controlling its texture and grain size. Moreover, the mass transfer process of metal ions and whiskers to the cathode surface could be enhanced at an appropriate electrodeposition temperature [24]. Therefore, it is indispensable to evaluate the effect of electrolyte temperature on the composite for the purpose of obtaining the high content of SiCw in the coating with smooth and uniform morphology. In this experiment, the temperature of electrolyte ranging from 20 ◦C to 60 ◦C were investigated at a stirring speed of 200 rpm and a current density 18 mA/cm2. As illustrated in Figure 2, the dendrites and nodules could be observed clearly at a lower temperature (Figure 2a), while the smooth and luminous surface accompanied by the increasing whisker content was obtained at a temperature rising to 50 ◦C (Figure 2d).In addition, the nickel grains became finer, and the composite coatings were also strengthened after adding the whiskers. However, the whiskers tended to aggregate and decreased its content in composite coatings when the temperature was exceeding 50 ◦C, which resulted in the deposited grains was coarse and oversized.

**Figure 2.** SEM images of Ni/SiC*w* composite coating fabricated from the electrolyte at stirring speed of 200 rpm, and current density of 18 mA/cm2, and temperatures of (**a**) 20 ◦C, (**b**) 30 ◦C, (**c**) 40 ◦C, (**d**) 50 ◦C, (**e**) 60 ◦C.

In general, both the average kinetic energy of ions and the viscosity of electrolytic could be consolidated by increasing the electrolyte temperature, which all consequently facilities the transportation of whiskers. Besides, the diffusion mobility of metal ions adsorbed on the cathode surface would be increased on the electrode surface, which resulted in an increase of the whiskers contents in composites coating. Whereas, the anions diffusion rate, the overpotential and the electric field force of the cathode reduced when the temperature reached to a specific value, which was detrimental to the whiskers embedded into the nickel matrix. In addition, the electrolyte would evaporate seriously at a higher temperature, thus, affecting the deposition of active ingredients onto the cathode and coursing the surface of composites. As a result, the excessively high temperature led to a reduction of the SiC*w* adsorption ability and the cathode electric field force, which ultimately reduced the SiCw content in the composite coatings.

In order to prevent the precipitation and agglomeration of SiCw in the electrolyte, mechanical stirring was adopted as a key technology in electroplating process, due to the high surface energy of the SiCw. Particularly, the content of particles in composite coatings is depended largely on the strength of stirring, which is advantageous to the uniform distribution of particles in electrolyte and the successful transportation to cathode surface [23]. Figure 3 shows the SEM images relating to the morphology of Ni/SiCw coatings fabricated at a different stirring speed. The conditions of current density and temperature were 18 mA/cm<sup>2</sup> and 50 ◦C, and the stirring speed was controlled at 0, 100, 300, and 500 rpm, respectively. Due to the gravity, a large number of whiskers accumulated and covered the major cathode surface without stirring, which resulted in an increase of the cathode overpotential and a decrease of the discharged metal ions near the cathode. Therefore, metal ions discharged at the corner and appeared dendrites or nodules in the composite coatings under the condition without stirring. With the increasing strength of stirring, the mass transfer rate of electrolyte and the effective concentration of whisker were also increased accordingly, which could promote the transportation of SiC*w* and increase the whisker content in the composite coatings. However, the electrolyte flew at a high speed when the stirring strength is higher than 300 rpm, which caused intense shock to the cathode surface, led to high-speed transfer of the electrolyte accompanied by whiskers and suppressed the attachment of whiskers to the cathode surface. More seriously, whiskers that were not fully embedded in the nickel matrix might be released from the cathode surface to the electrolyte, which resulted in a dramatic decrease of the whisker content in the composite coatings.

**Figure 3.** SEM images of Ni/SiC*w* composite coating fabricated from the electrolyte at a current density of 18 mA/cm2, and temperatures of 50 ◦C, and stirring speed of (**a**) 0 rpm, (**b**) 100 rpm, (**c**) 300 rpm, (**d**) 500 rpm.

The cross-sectional SEM images of Ni/SiCw composite coatings prepared with different electrodeposited conditions were shown in Figure 4. From Figure 4a,b, in which the current density ranges from 2 to 8 mA cm−2, an increase in the number of SiCw (pointed by an arrow), which are fairly well dispersed in the nickel matrix, can be clearly observed on the fracture surface. According to the micrographs of the cross section, the fracture surface does not represent the initial interface between the the nickel matrix and SiCw after breaking with external force, so the interface was further illustrated after treatment by an ion beam thinner (Figure 4c,d). From the cross-sectional images of the thinned surface, the whiskers were cut in an identical plane, but in different orientations, due to their random distribution in the the nickel matrix; meanwhile, multiple shapes of whiskers after thinning by ion beam thinner were also noticed.

**Figure 4.** SEM images of the cross section for the fracture (**a**,**b**) and thinned surface (**c**,**d**) of the electrodeposited Ni/SiCw composite coatings (**a**) SiCw 0.8 g L<sup>−</sup>1, 2 mA/cm2, 35 ◦C, 150 rpm; (**b**) SiCw 0.8 g L<sup>−</sup>1, 18 mA/cm2, 35 ◦C, 150 rpm; (**c**) and (**d**) represent the thinned surface of (**a**) and (**b**), respectively.

The EDS spectra of Ni/SiCw composite coatings, prepared at different conditions, are illustrated in Figure 5; and is related to the EDS spectra, which has been listed, as shown in Table 2. As shown in Figure 5a, only Ni peaks were observed in the sample without SiCw coating, whereas the peaks of Ni, Si and C were all collected in Figure 5b–e. The EDS results indicated that the SiC whiskers were incorporated into the nickel matrix successfully. Besides, the obviously increased content of Si and C could be seen from Figure 5b,c, which was resulted from the increased current density from 2 to 18 mA/cm2. Since the co-deposition rate of SiCw was strongly related to that of Ni matrix in the process of electrodeposition, then increasing current density could promote the deposition rate of matrix metals and save the time of SiCw embedding into the Ni matrix. In Figure 5c,d, the EDS spectra of the Ni/SiCw coatings deposited at 35 ◦C and 50 ◦C, which showed that the higher temperature facilitated the increase of whiskers content in electrodeposits. It shows that the higher temperature facilitates the increase of whiskers content in electrodeposits. As to the influence of stirring speed on the content of SiC*w*, the stronger stirring speed can contribute to the whiskers transmitting to the cathode surface from the electrolyte which results in its higher content in coatings (Figure 5e).

**Figure 5.** EDS spectra of electrodeposited Ni/SiC*w* composite coatings (**a**) without SiC*w*, 2 mA/cm2, 35 ◦C, 150 rpm; (**b**) SiC*w* 0.8 g L<sup>−</sup>1, 2 mA/cm2, 35 ◦C, 150 rpm; (**c**) SiC*w* 0.8 g L<sup>−</sup>1, 18 mA/cm<sup>−</sup>2, 35 ◦C, 150 rpm; (**d**) SiC*w* 0.8 g L<sup>−</sup>1, 18 mA/cm2, 50 ◦C, 150 rpm; (**e**) SiC*w* 0.8 g L<sup>−</sup>1, 18 mA/cm2, 50 ◦C, 300 rpm.


**Table 2.** Weight fraction of the composite coatings as deposited at different conditions.

#### *3.2. E*ff*ects of Electrodeposition Parameters on Crystallographic Texture of the Ni*/ *SiCw Composite Coatings*

For all the electrodeposited coatings, the nickel crystals are face-centered cubic structures (FCC). The three peaks at 2θ = 44.5, 51.85 and 76.37◦ that corresponds to (111), (200) and (220) crystallographic planes of nickel respectively, are in agreement with the standard XRD pattern JCPDS 04-0850. Note that the peaks corresponding to the SiC*w* were not detected in the XRD pattern of the composites. This might be due to the high content of Ni and its high scattering factor leading to higher peak intensities compared to those of SiC*w* [25]. Figure 6 shows the XRD pattern for the applied SiC whiskers. Comparing to the standard card, SiC particles was fit to No. 29-1129 exactly.

**Figure 6.** XRD pattern of the SiC whiskers.

Figure 7 shows the XRD patterns and its crystallographic orientation indexes of the Ni/SiC*w* composite coatings fabricated at a different current density. As shown in Figure 7a, at low current density (2 mA·cm<sup>−</sup>2), the composite coating exhibits a strong (111) plane preferred orientation and a relatively low intensity of peaks (200) and (220). However, the intensity of (111) and (220) peak gradually reduce with the increase of current density, while the (200) peak increased obviously. Particularly, when the current density increased to 40 mA·cm−2, the value of (200) peak intensity reaches the strongest, and the other peaks are almost disappeared. These results mean that the current density alters the texture of the composites and enhances the (200) reflections, which is in good agreement with the literature reported observation (the current efficiency of Ni deposition reached high values of over 97%) [26]. XRD analyses show that the texture coefficient of the deposits varies with current density. From Figure 7b and Table 3, the calculated results exhibit that the texture coefficient of (111) and (220) orientation decreases from 0.91 to 0.05 and 0.43 to 0.18, respectively; while the (200) orientation index increases from 0.66 to 2.9.

The variation of preferred orientation implied two possible reasons. One is the growth axis of the nickel grain mainly depends on the applied current density. There is a trend for dominating surface textures with increasing current density, as follows—(110) < (211) < (100). This means that (110), (211) and (100) texture are preferred orientation at very low current, medium and high current, respectively [27]. This also can be explained by the theory of nucleation energy (W*hkl*) [28,29] as Equation (2).

$$\mathcal{W}\_{hkl} = \frac{B\_{hkl}}{\frac{zF}{\text{N}\text{a}}\eta - A\_{hkl}} \, ^{\prime} \tag{2}$$

where A*hkl* and B*hkl* are material parameters of the crystallographic direction (*hkl*); Na and F are the constant of Avogadro and Faraday, respectively; z is the number of electrons, and η is the over-voltage of cathodic polarization. For face-centered cubic metals, at low over-voltage, W111 < W100 < W110, while W100 < W111 at high over-voltage [26]. This phenomenon indicates that when particular nucleation energy is lower than other types of nucleation energy, the crystals of the deposited metal will exhibit particular preferred orientation. The other is the incorporated whiskers provides the occurrence of new nuclei and limits the growth of the original crystal grains [30,31], which facilitates more nucleation sites and the grain size refined. As mentioned above, the crystalline structure of Ni-W/SiC deposits is also affected by the factor of the SiC*w* content in the composites.

**Figure 7.** (**a**) XRD patterns of composites fabricated at a different current density, (**b**) The crystallographic orientation indexes of (111), (200) and (220).

**Table 3.** Texture coefficients of various (*hkl*) for the composite coatings as deposited at different current density.


Figure 8 shows the XRD patterns and its crystallographic orientation indexes of the Ni/SiC*w* composite coatings electroplated at 20, 30, 40, 50, 60 ◦C, respectively. With the temperature elevated, the intensity of the (200) increases, while (111), (220) peaks increases. The calculated orientation index

of (200) plan increases from 1.08 to 2.32, while the (111) and (220) orientation index decreases from 1.15 to 0.32 and 0.73 to 0.36 (Figure 8b and Table 4), respectively. These results indicate that relative high electrolyte temperature is benefit for the preferred orientation of (200) plan, which is in agreement with previous studies [32].

In the view of the observations made by E. Budevski et al. [28], the energy barrier of nucleation showed a positive linear dependence with 1/η and 1/η2 for the 2D and 3D model, respectively; the nucleation kinetics of nickel crystallites depends on the applied overpotential η. At low overvoltage, the particles near the cathode suppress metal ion reduction; while at high overvoltage, the transmission of the metal ions to the cathode becomes an important factor, moreover, the incorporated of inert whiskers enhance this transmission [23]. Therefore, the grains size of composites at high temperature becomes larger than that at low temperature (Figure 2).

**Figure 8.** (**a**) XRD patterns of composites fabricated at a different temperature, (**b**) The crystallographic orientation indexes of (111), (200) and (220).

**Table 4.** Texture coefficients of various (*hkl*) for the composite coatings as deposited at different temperature.


Figure 9 shows the XRD patterns and its crystallographic orientation indexes of the Ni/SiC*w* composite coatings electroplated under different stirring conditions. The XRD patterns (Figure 9a) shows that, under static conditions, the intensity of (200) peak is slightly higher than that of (220) peak, and (111) peak possess the highest intensity. When given a certain stirring speed at 100 rpm, there is about the same intensity of (111) and (200) crystal planes, while the intensity of (220) peak is very weak and almost disappeared. With increasing stirring speed to 300 rpm, the intensity of (200) peaks gradually increase, along with a relative decrease in intensity of (111) and (220) peaks. Notably, at high stirring speed (500 rpm), the intensity of (200) peak presents a weakening trend. This is because the excessive stirring speed takes a large impact force for the cathode surface, it can, therefore, caused the whiskers transferred with the electrolyte at a high speed. Ultimately, the content of whiskers incorporated into the nickel matrix is reduced. The XRD patterns further analysis is shown in Figure 9b and Table 5. The preferred orientation of the composite coating at different stirring speeds is mainly related to the content of SiC*w* in coatings, and thereby the embedded SiC whiskers into the nickel matrix change the microstructure of the composites.

**Figure 9.** (**a**) XRD patterns of composites fabricated at a different stirring speed; (**b**) the crystallographic orientation indexes of (111), (200) and (220).

**Table 5.** Texture coefficients of various (*hkl*) for the composite coatings as deposited at different stirring speed.


#### **4. Conclusions**

In summary, the influences of the external conditions of electrodeposition on the crystallographic texture, morphology and composition were investigated. For the deposited Ni/SiCw composite coatings, the content of SiC*w* in coatings increased with the increasing current density, stirring rate and electrolyte temperature in a certain range. Particularly, glossy Ni matrix surface with uniformly distributed whiskers was obtained at a current density of 18 mA/cm2, whereas it became nodular and inhomogeneous with the widely differed spherical particles. Moreover, the detailed XRD analysis indicated that the current density altered the texture of the composites and enhanced the reflections of (200) plane. Also, the optimized temperature results revealed that the flat and bright morphology, as well as high content of whiskers in the composite coatings, was obtained at a temperature of 50 ◦C. Furthermore, the effects of stirring speed on the composition and structure on the coatings were studied. It was found that when the stirring speed was too low, the composite coatings appeared dendrites or nodules, while the whisker content was drastically reduced at a high speed. Consequently, these external conditions not only affect the surface morphology and the SiCw content of the composites, but also modified their preferred orientation.

**Author Contributions:** Conceptualization, G.D. and L.L.; methodology, L.L. and H.W.; software, L.L. and Y.S.; validation, G.D.; formal analysis, L.L. and H.W.; investigation, Z.Y. and Y.S.; resources, G.D.; data curation, L.L. and G.D; writing—original draft preparation, L.L.; writing—review and editing, L.L. and G.D.; visualization, H.W.; supervision, Y.Z.; project administration, G.D.; funding acquisition, G.D.

**Funding:** This research was funded and supported by the National Natural Science Foundation of China (No. 61571287).

**Acknowledgments:** The authors would like to thank supports from the Shanghai Professional Technical Service Platform for Non-Silicon Micro-Nano Integrated Manufacturing.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Optimization of Processing Parameters for Water-Jet-Assisted Laser Etching of Polycrystalline Silicon**

#### **Xuehui Chen, Xiang Li, Chao Wu, Yuping Ma, Yao Zhang, Lei Huang \* and Wei Liu**

College of Mechanical and Electrical Engineering, Anhui Jianzhu University, Hefei 230601, China; xhenxh@163.com (X.C.); lix1117@163.com (X.L.); 15212428852@163.com (C.W.); wxlmyp@163.com (Y.M.); nianmin134@163.com (Y.Z.); wliu@hfcas.ac.cn (W.L.)

**\*** Correspondence: huangl75@ahjzu.edu.cn; Tel.: +86-137-2110-2968

Received: 11 April 2019; Accepted: 1 May 2019; Published: 8 May 2019

**Abstract:** Liquid-assisted laser technology is used to etch defect-free materials for high-precision electronics and machinery. This study investigates water-jet-assisted laser etching of polysilicon material. The depths and widths of the etched grooves were investigated for different water-jet incident angles and velocities. To select optimal parameters for a composite etching processing, the results of many tests must be compared, and at least one set of good processing parameter combinations must be identified. Herein, the influence of different parameters on the processing results is studied using an orthogonal test method. The results demonstrate that the depths and widths of the processing grooves were nearly identical at water-jet angles of 30◦and 60◦; however, the 60◦ incidence conferred a slight advantage over 30◦ incidence. The section taper, section depth, and surface topography were optimized at a water-jet velocity of 24-m/s, 1.1-ms laser pulse width, 40-Hz frequency, and 180-A current. Under these conditions, the section taper and groove depth were 1.2◦ and 1.88 mm, respectively. The groove surfaces exhibited no splitting, slagging, or other defects, and no recast layers were visible.

**Keywords:** laser etching; water jet; polycrystalline silicon; orthogonal test

#### **1. Introduction**

Polysilicon material is a brittle semiconductor material with high wear resistance, hardness, chemical stability, and low thermal conductivity. The superior performance of silicon materials has secured their use in photovoltaics, aerospace, electronics, and other industries. However, owing to their physical properties, polysilicon and similar semiconductor materials are difficult to process effectively via conventional machining, i.e., the quality of processing may be poor. This problem can be resolved by employing laser processing technology. Such technology is widely used in electronic and other industries for precision machining of difficult-to-machine metal and non-metal material with high hardness, strength, toughness, and brittleness [1–3]. Unfortunately, under ambient air or inert gas conditions, the transient action of laser heat in the laser processing of materials inevitably produces microcracks and slag on the surface of the material, thereby reducing the qualification rate and processing quality of precision device products. To ameliorate the shortcomings of traditional lasers in processing such materials, liquid-assisted laser technology has been proposed herein. In this technology, with a hard and brittle material, the liquid reduces the temperature gradient at the processing site, thereby avoiding cracking and heat-affected areas of the material and improving the processing quality. Over time, liquid-assisted laser processing technology has branched into water-guided, underwater, and water-jet-assisted laser processing technology. Li et al. [4] etched trenches into silicon using a water-guided laser micromachining method. They simulated the melting and removal of silicon

under the cooling effect of the water jet and verified their model via comparative simulations and experiments. Their model considers the effects of cutting speed and other factors on the surface quality and heat-affected zone of the workpiece. Porter et al. [5] applied water-guided laser processing technology to cut metal sheets and summarized the influence of some processing parameters on the results. Hock et al. [6] cut stainless-steel plates and brass sheets (thickness ≤ 100 μm) using conventional laser and water-guided laser cutting technologies and compared the cutting widths, heat-affected zones of the two products. Compared to the conventional technique, the water-guided laser cutting system achieved a lower heat-affected zone, slag accumulation height, and slit width. Adelmann et al. [7] cut aluminum, titanium, steel, and other metals using a water-guided laser and investigated the effects of laser power, pulse repetition frequency and etching times on the processing quality. The maximum depths of the cut aluminum, titanium, and steel were 8, 4.7, and 1.5 mm, respectively, and the maximum aspect ratios were 66.7, 39.2 and 12.5, respectively. Using Nd:YAG pulse fiber lasers, Choubey et al. [8] conducted stainless-steel plate experiments in both air and underwater environments. Their results demonstrate that compared to air cutting, underwater cutting improves the sample and slightly reduces its slit width, heat-affected zone, and adhesion slag. Mullick et al. [9] established a model for investigating the energy-loss mechanism in underwater laser cutting. The model accounts for the interaction between the laser, material, and water. Most of the energy is lost by the scattering of the laser light in water vapor. The losses account for 40–50% of the total laser energy. The correctness of the model results was verified via experiments. Tsai [10] processed LCD glass materials using underwater and conventional laser techniques and investigated the effects of the laser parameters, material size, and liquid type on the removal rate, slit width, incision depth, taper, and surface roughness of the processed materials. Their results demonstrated the advantages of underwater laser processing over conventional laser processing. Charee et al. [11] processed silicon wafers via a flowing underwater laser composite etching process. To study the effects of underwater laser etching, they controlled the rate and direction of the water flow in a closed water chamber. With a high water flow rate, the grooves were deepened and few re-agglomeration layers were formed. The etching effect was superior to that obtained in still water. Kalyanasundaram [12] combined the parameters of laser water jets with those of hard and brittle processing materials and significantly improved the quality of the cut. Tangwarodomnukun et al. [13,14] compared the surface morphologies of samples prepared via hybrid laser–water jet micromachining, traditional silicon processing, and laser composite processing technologies. They varied the processing parameters and analyzed their effect on the heat-affected zone and processing quality. Bao et al. [15] established a fluid dynamics model based on smoothed particle hydrodynamics experimentally investigated the water-jet-assisted laser cutting of silicon materials. Experimental analysis revealed that the laser ablation of silicon mainly occurs via explosive melting. Zhu et al. [16] established a numerical model of heat transfer and material ablation in the water-jet-assisted laser etching of single-crystal germanium. They reported that the irradiated material can be discharged through the water jet and that during the non-pulsing periods, the water-jet cooling effectively removes heat buildup in the workpiece, thereby minimizing the thermal damage caused by laser heating. They also analyzed the influence of the processing parameters on the ablation process. Increasing the laser pulse energy deepened the grooves, whereas increasing the applied water pressure reduced the threshold workpiece temperature for material removal. Feng et al. [17] established a three-dimensional analytical model of the temperature field during the hybrid laser water-jet micromachining process. Their model considers the interaction between the laser, water jet, and workpiece. The absorption of laser light by the water, the formation of laser-induced plasma in water, the formation of bubbles, and the laser refraction at the air–water interface were discussed. To evaluate two machining methods, Hao Zhu et al. [18] compared the direct and chemical-assisted picosecond laser trepanning of single crystalline silicon. In their study, an orthogonal test design scheme was adopted to consider the relevant parameters affecting the trepanning process. They found that the direct laser trepanning results are associated with significant thermal defects, whereas the chemical-assisted method can process microholes with negligible thermal damage.

Conventional gas-assisted laser processed materials have defects such as a large heat-affected zone and many re-casting layers. Compared with traditional gas-assisted laser processing, the laser processing area is cooled and impacted by the water jet. Specifically, during laser processing, the water jet cools the material and its impact can wash away the slag and recast a layerafter laser processing. The cooling action can reduce excess heat acting on the surface of the material, reduce the heat-affected zone, and reduce the generation of microcracks, thereby improving the surface quality of the material. In this paper, water jets with different angle and velocities were used to assist laser etching of polysilicon. Previously, we studied the effect of different water jet velocities on the processing grooves when the water jet angle was 30◦ [19]. However, to the best of our knowledge, the effects of changing both the water-jet angle and velocity on the depth and width of the etched grooves have not been investigated. This paper builds on previous studies by varying the jet angles and velocities in water-jet-assisted laser etching of polysilicon materials and examining its effect on the depth and width of the etched grooves. Subsequently, after fixing the water-jet angle at 60◦, the optimal processing parameters in the experimental grooves are investigated in an orthogonal optimization experiment.

#### **2. Theoretical Analysis of the Influence of Water-Jet-Assisted Laser Etching**

Water-jet-assisted laser processing introduces a water jet to traditional laser processing. The auxiliary gas removes slag and debris generated in the laser processing, thereby improving the processing quality of the materials and reducing slag accumulation, the generation of microcracks, formation of recast layers during processing, and the heat-affected zone. Figure 1 schematizes the water-jet-assisted laser technique. Note that the water jet consumes part of the laser energy. In contrast, when jetted at appropriate velocity, the water jet washes away slag and the recast layer in an appropriate time, which reduces the heat-affected zone and improves processing quality. Water-jet-assisted laser processing also reduces photothermal efficiency and suppresses secondary adhesion of the etched material via laser-induced liquid cavitation, water cooling, and water flow.

**Figure 1.** Schematic of water-jet-assisted laser processing.

Conventional gas-assisted and water-jet-assisted laser processing techniques remove material via laser irradiation, which smelts or vaporizes the material surface. In the conventional technique, most laser energy is absorbed by the material. The material's thermal properties will transfer excess heat around the laser action zone, which enlarges the heat-affected zone and increases slag accumulation. In water-jet-assisted laser etching technology, conventional laser etching is supplemented by water jets applied at a certain angle and velocity. This technology combines the high efficiency of traditional gas-assisted laser processing with the advantages of water jets. The laser primarily acts as a softening

material that melts or vaporizes the material as it is heated. The water jet strikes the laser processing area and cools it by heat transfer. Generally, the specific heat capacity of water is greater than that of the material matrix, which implies that water absorbs excess heat from the material surface and provides a cooling effect. In addition, increasing the water-jet velocity increases the impact force on the material. Within a certain range, the impact force is sufficient to flush the slag and other defects generated by the laser action. Then, the flowing water removes the slag or recast layer from the tank in sufficient time. Note that reducing slag and the recast layer from the bottom and walls of the trench can improve processing quality [20].

In this study, the effects of water jet pressure at different angles and velocities on the material were compared in fluent simulations. Herein, the water-jet angle was set to 30◦ or 60◦, and the velocity was varied between 16, 20, 24, and 28 m/s (Figure 2). Note that many factors were simplified during simulation; thus, the results are theoretical rather than practical. As shown in Figure 2, increasing jet velocity at a fixed incident angle gradually increased the impact force of the water jet on the material. In addition, for a fixed water-jet velocity, the impact force was ~1.5 times higher at an impact angle of 60◦ than at 30◦. Thus, it is theoretically demonstrated that 60◦ water-jet-assisted laser etching of the polysilicon material demonstrates groove depth and width greater than that of the 30◦ water-jet-assisted laser etching polysilicon material. In actual machining, greater jet velocity is advantageous for material removal; however, if the water flow velocity is excessive, the water jet will shake the experimental device, which will also deflect the impact center of the jet on the material. A significantly diverted water beam will cause splashing, which will generate "water mist" that affects the focused laser spot and incurs large laser energy loss.

**Figure 2.** Impact force of a water jet striking a material at different velocities and impact angles.

#### **3. Experimental Methods and Analysis**

An HGL-LMY500 solid laser cutting and welding processing system (Wuhan Huagong Laser Engineering Co., Ltd., Wuhan, China) was used in the test. The HGL-LMY500 comprises an Nd:YAG solid-state pulsed laser (average output: 500 W), a power supply system, an optical system, a control system, a cooling system, a numerical control table, and a computer (Figure 3). During processing, this system produces a laser with a wavelength of 1064 nm, a multimode circular spot, and a spot diameter of 0.2 mm. By adjusting the current, pulse width, and repetition frequency, the laser output energy, laser light intensity, beam quality, etc. satisfy processing requirements. The main parameters of this laser system are shown in Table 1. Here, the water jet device is a water jet processing system designed by our team. Note that the water-jet velocity and angle be adjusted. In addition, the water jet nozzle (water jet nozzle diameter is 0.7 mm) can be fixed to ensure that the water jet can impact and cool the

material in the laser processing area. According to the actual work requirement, the motor, liquid tank, filter device, inlet pipe, plunger pump, overflow return control valve, pressure gauge, accumulator, one-way control valve, and jet nozzle and its fixing device.

**Figure 3.** Schematic diagram of the experimental Nd3<sup>+</sup>: YAG laser etching apparatus.

**Table 1.** Laser system processing parameters.


The experimental sample was a polycrystalline silicon plate (length × width × height: 20 <sup>×</sup> 20 <sup>×</sup> 2 mm3). The effect of different water jet angles and velocities on the depth and widths of the processing grooves in water-jet-assisted laser etching was investigated. Prior to experiment, it is necessary to ensure that the following laser parameters remain unchanged: 40-Hz laser repetition frequency, 0.6-ms pulse width, 1-mm/s scanning speed, and 140-A current. Under this processing parameter, the laser irradiation intensity reaches 5.5 <sup>×</sup> 105 W/cm2. When the irradiation intensity is 5.5 <sup>×</sup> 105 W/cm2, it is sufficient to achieve the energy density required for irreversible damage to the surface material, and the material is removed by the action of the laser. In water-jet-assisted laser processing, the jet is impacted slightly behind the laser beam (1 mm behind processing area) to avoid direct contact between the jet and laser beam, avoid excessive loss of laser energy, and boost cooling and impact performance. The cross-sections and widths of the resulting grooves are shown in Figures 4 and 5, respectively. Figures 6 and 7 show the groove depths and widths under different water-jet velocity (horizontal axis) and impact angle (colors), respectively.

**Figure 4.** Cross-sections of grooves obtained by etching polycrystalline silicon at different impact angles (top row: 30◦; bottom row: 60◦) and jet velocities (left to right: 16, 20, 24, and 28 m/s) (25×).

**Figure 5.** Widths of grooves obtained by laser etching at different impact angles (top row: 30◦; bottom row: 60◦) and jet velocities (left to right: 16, 20, 24, and 28 m/s) (30×).

**Figure 6.** Groove depths obtained at different water-jet velocities (horizontal axis) and impact angle (colors).

**Figure 7.** Groove widths obtained at different water-jet velocities (horizontal axis) and impact angles (colors).

As shown in Figures 4–7, the depth and width of the laser-etched grooves in the polysilicon material are dependent on water-jet velocity but not impact angle (the water-jet velocity trends were very similar at both 30◦ and 60◦). The depth of the processing groove initially increased with increased jet velocity from 16 to 24 m/s but decreased between 24 and 28 m/s. In addition, the width of the processing groove gradually increased across the entire range of tested jet velocities. For the same water-jet velocity, the groove depths and widths were slightly larger at an impact angle of 60◦ than at 30◦. The slight increase at 60◦ is attributable to the higher impact force at 60◦ than at 30◦. As the water-jet velocity increased from 16 to 24 m/s, both the impact and cooling effect of the water-jet increased; however, the increasing trend of the cooling effect slowed down. As a result, the total increasing trend was dominated by the impact of the water jet, whereas the cooling effect played a secondary role. The molten silicon and slag deposited on the bottom of the tank were removed to reduce laser energy absorption; thus, the depth and width of the processing grooves increased gradually with water-jet velocity. However, at a water-jet velocity of 28 m/s, the large impact force was compromised because the experimental device was disturbed. At this velocity, the impact of the water jet deviated from the center of the jet impact, and the water beam diverged when striking the material surface, which caused splashes and "water mist" that affected the focused laser spot. In water-assisted laser processing, the material surface is covered with a thin, flowing layer of water, which absorbs some of the laser energy, thereby weakening the laser's action on the material. Simultaneously, the enhanced convective heat transfer between the laser thermal energy and surrounding medium (i.e., the water) incurs large laser energy loss, which reduces the processing depth.

#### **4. Orthogonal Test Design**

#### *4.1. Design Table of Orthogonal Test Plan*

The orthogonal test method was employed to assess the quality of the water-jet-assisted laser etching of the polysilicon samples. In this experiment, the processing quality was primarily evaluated by examining the section taper, section morphology, and processing depth. At the largest processing depth, the processing tank must have the smallest taper and fine surface morphology. The cross-section of the sample trough was observed using a VM-3020E 2D imager, and the taper and depth of the section were measured. The results are shown in Figure 8.

**Figure 8.** Cross-sectional profile of the groove in the polysilicon sample after processing.

Section 3 discussed the effect of water-jet angles and velocities on the depth and width of the grooves fabricated by laser-assisted etching of polysilicon material. Note that many interacting processing parameters are involved in the actual processing. The parameters that affect water-jet-assisted laser etching are the laser and water-jet parameters. To improve the etching effect, one or several sets of good combinations of processing parameters must be determined, which requires numerous experiments. The orthogonal test method optimizes the selection of processing parameters, which reduces the number of tests and the time required to select suitable processing parameters. As shown in Figure 2, the impact effect of the water jet at any incident velocity was slightly higher at 60◦ compared to that at 30◦. Therefore, in the orthogonal test, the incident angle of the water jet was set to 60◦ and the other parameters were set according to those of the actual situation. Here, the main research objects were the laser processing parameters (laser pulse width, repetition frequency, and pulse energy) and pulse energy (represented by input current).

In the orthogonal test, the water-jet velocity, laser pulse width, laser repetition frequency, and laser input current were labeled *A*, *B*, *C*, and *D*, respectively, and each factor was assigned four levels (1, 2, 3, and 4), as shown in Table 2.


**Table 2.** Factor/level table.

The water-jet-assisted laser etching process was tested according to the L16 (45) orthogonal table (presented as Table 3). The design included 16 orthogonal experimental schedules. The subscripts in each plan indicate the level number (first column in Table 1). For example, A1 represents jet velocity of 16 m/s and B4 indicates a pulse width of 1.1 ms.


**Table 3.** Orthogonal test table (L16 (45)).

#### *4.2. Experimental Results and Optimization Options*

Following the test scheme shown in Table 3, an orthogonal test was conducted on the water-jet-assisted laser etching processing platform. Here, the measured data were the depth and taper of the groove of the given evaluation index. The 16 test results and their corresponding data are shown in Figures 9 and 10, respectively.

**Figure 9.** Orthogonal test results (15×). Label "a" represents the depth and taper of the section, and "b" represents surface topography.

**Figure 10.** Experimental values of orthogonal tests 1–16. (**a**) Taper value; (**b**) Depth value.

The test results for each evaluation index in Table 4 were calculated and analyzed using the range analysis method. After computing test indices (i.e., section taper and section depth, corresponding to the *m*-th level of the *j*-th factor) and their average values (*Kjm*), the range differences *Rj* of the two evaluation indexes were calculated for the *j*-th factor. The results are shown in Table 4.


**Table 4.** Range calculation values of orthogonal test results.

From the average values of the indicators corresponding to each level of each factor, the appropriate level of each factor can be summarized as follows:


To analyze the test results visually and intuitively and to obtain comprehensive conclusions, the relationships between the evaluation indicators and various factors (effect curves) were derived from the range analysis data. In this experiment, the average corresponding experimental indicator was plotted for each level in a given factor. The results of factors *A*, *B*, *C*, and *D* are shown in Figures 11–14, respectively.

**Figure 11.** Effect of water-jet velocity on (**a**) cross-section taper and (**b**) cross-section depth.

**Figure 12.** Effect of laser pulse width on (**a**) cross-section taper and (**b**) cross-section depth.

**Figure 13.** Effect of laser repetition frequency on (**a**) cross-section taper and (**b**) cross-section depth.

**Figure 14.** Effect of current on (**a**) cross-section taper and (**b**) cross-section depth.

The range analysis results given in Table 5 demonstrate that the section taper was most influenced by factor *B*, followed by factors *D*, *C*, and *A*. In addition, the actual processing requires that a smaller section taper will result in a larger section depth.

**Table 5.** Range analysis results of orthogonal test.


The primary and secondary factors, optimal level, and optimal combination (i.e., the optimal solution) were obtained (Table 6) by combining the range calculation results given in Table 4 with the relationships between evaluation indicators and factors shown in Figures 11–14.

**Table 6.** Processing parameters in optimized schemes of each evaluation index.


As shown in Table 6, the optimal process that minimizes the cross-section taper combines a 1.1-ms laser pulse width, 40-Hz frequency, and 180-A current with the 28-m/s water-jet velocity. Under this condition, the cross-section taper after composite processing was 0.3◦ (Figure 15a). The optimum process that maximizes the tank depth combines a 1.1-ms laser pulse width, 40-Hz frequency, and 150-A input current with the 20-m/s water-jet velocity. With this configuration, the etched tank had a depth of 2.05 mm but was shaped like a teardrop with a non-ideal taper (Figure 15b).

**Figure 15.** Cross-sectional morphology of samples after composite processing at processing parameters that optimize each evaluation index (×30): (**a**) section taper and (**b**) section depth.

Based on the above analysis, the final combination of optimal processing parameters was determined to be *B*4*A*3*C*3*D*4. The combination of processing parameters that optimized the cross-section taper and cross-section depth of the etched tank were as follows: jet velocity of 24 m/s, laser pulse width of 1.1 m/s, laser frequency of 40 Hz, and input current of 180 A. The cross-section and surface topography of the groove after composite processing under these conditions are shown in Figure 16. As can be seen, the surface of the groove body is free of splitting and slagging defects and has no recast layer. Here, the taper and depth of the groove were 1.2◦ and 1.88 mm, respectively. Although the etching result of each index was slightly worse than the result of optimizing a single indicator, it generally satisfies the etching effect requirements.

**Figure 16.** Surface and section morphologies of groove formed by composite machining under final optimized combination of processing parameters (×30): (**a**) sectional topography and (**b**) surface topography.

#### **5. Conclusions**

A water jet is frequently employed to assist laser etching of polysilicon materials. This study assessed the influence of water-jet incident angle and velocity on the depth and width of the etching groove. At a water jet incident at 60◦, the effects of processing parameters (i.e., water-jet velocity, laser pulse width, repetition frequency, and input current) on the surface quality of the composite etching were investigated using the orthogonal test method. The processing parameters were optimized and verified. By comparing the verification test results with those of the orthogonal test, the processing parameters that optimize each evaluation index were obtained. Our primary findings are summarized as follows.

(1) The effects of 30◦ and 60◦ water-jet-assisted laser etching of polysilicon materials on the depth and width of the etched trenches were investigated. It was found that when the water jet angle was maintained at one value, the depth of the machining tank initially increased and then decreased as water-jet velocity increased, and the width of the machining tank increased gradually. As the velocity increased from 16 to 24 m/s, the depth of the processed groove increased, but a further velocity increase (28 m/s) reduced the groove depth. The width of the processing groove was a gradually increasing function of the water-jet velocity. When the water-jet velocity was maintained at one value, the 60◦ water jet assisted the laser etching of the polysilicon material to obtain a tank with depth and width greater than that obtained by the 30◦ etching.

(2) Taking the 60◦ water-jet-assisted laser etching of polysilicon as an example, the orthogonal experimental method was used to optimize the processing parameters. In this study, the processing quality was evaluated relative to the depth, taper, and surface topography of the processing tank (the deeper the depth of the tank, the smaller the taper and the better the surface topography).

The experimental results demonstrated that the section taper was minimized under a 1.1-ms laser pulse width, 40-Hz frequency, and 180-A current when assisted by 28-m/s water-jet velocity. Here, the minimum section taper was 0.3◦. The tank depth was maximized with a 1.1-ms laser pulse width, 40-Hz repetition frequency, and 150-A input current when assisted by a 20-m/s water-jet velocity. Here, the maximum tank depth was 2.05 mm; however, the section taper was not ideal. The optimal processing parameters were a laser pulse width of 1.1ms, repetition frequency of 40 Hz, water-jet velocity of 24 m/s, and input current of 180 A. Under these conditions, the cross-section taper and groove depth were 1.2◦ and 1.88 mm, respectively, and the groove surface demonstrated no defects.

**Author Contributions:** X.C. designed the study. Y.M., Y.Z., C.W., W.L. and L.H. developed the methodology. X.L. wrote manuscript.

**Funding:** This research was funded by Provincial Key Project of Natural Science Research Anhui Colleges (grant number: KJ2015A013 and KJ2015A050), Anhui Province Outstanding Young Talent Support Program Key Project (grant number: xyqZD2016153), National Natural Science Foundation of China (grant number: 51175229).

**Acknowledgments:** We thank all the authors for their joint efforts to complete the experiment. We also thank Jiangnan University for providing experimental equipment.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Multi-Objective Design Optimization of the Reinforced Composite Roof in a Solar Vehicle**

#### **Ana Pavlovic, Davide Sintoni, Cristiano Fragassa \* and Giangiacomo Minak**

Department Industrial Engineering, University of Bologna, via Fontanelle 40, 47121 Forlì, Italy; ana.pavlovic@unibo.it (A.P.); davide.sintoni@studio.unibo.it (D.S.); giangiacomo.minak@unibo.it (G.M.) **\*** Correspondence: cristiano.fragassa@unibo.it; Tel.: +39-347-6974046

Received: 16 March 2020; Accepted: 10 April 2020; Published: 12 April 2020

#### **Featured Application: A photovoltaic roof for vehicles, in sandwich-structured composite, is designed for optimizing static sti**ff**ness and dynamic response, but also energy e**ffi**ciency.**

**Abstract:** A multi-step and -objective design approach was used to optimize the photovoltaic roof in a multi-occupant racing vehicle. It permitted to select the best combination of design features (as shapes, widths, angles) in composite structures simultaneously balancing opposite requirements as static strength and dynamic stiffness. An attention to functional requirements, as weight, solar cells cooling and solar energy conversion, was also essential. Alternative carbon fiber-reinforced plastic structures were investigated by finite elements using static and modal analyses in the way to compare several design configurations in terms of natural frequencies, deformations, flexural stiffness, torsional stiffness, and heat exchange surfaces. A representative roof section was manufactured and tested for model validation. A significant improvement respect to the pre-existing solar roof was detected. The final configuration was manufactured and installed on the vehicle.

**Keywords:** design optimization; solar vehicles; photovoltaic roof; lightweight structures; carbon fiber-reinforced plastic (CFRP); natural frequencies; stiffness; heat exchange; Ansys ACP

#### **1. Introduction**

The design process in engineering is often based on a cyclical path aiming to improve existing solutions. However, this step-by-step approach becomes rather complex to be applied when the need exists to consider different optimization criteria at the same time. Such a situation can occur in many real applications [1], but it comes to be essential in the case of competition vehicles.

Sports car design surely represents one of the most exciting test benches for designers to set up efficient design methodologies and solutions. The desire of competitiveness commonly drives the designer to direct the design action towards a multi-objective optimization (MOO).

In the case of a Formula One car [2], for instance, the opportunity to lighten the vehicle comes up against the need to offer higher lift [3] and, then, requires a proper balance between weights and loads [4]. In these terms, the field of solar races is maybe even more interesting, where the power involved is extremely low and the optimization criteria must be even more shrewd [5].

Solar cars are electric vehicles where the energy for powertrain is provided by the Sun thanks to the installation of photovoltaic panels (usually on the roof) and a battery pack for energy storage. Thus, energy efficiency definitely represents the key concept driving each design choice in the case of solar prototypes [6]: aerodynamics [7], dynamics [8], kinematics [9] and even manufacturing [10].

However, weight minimization probably stands for one of the most deemed design policies with the scope to minimize the inertial masses and energy losses. In accordance with [11] based on electric cars, in fact, an amount of 13% can be saved in term of energy if a 10% decrease in weight is achieved. Weight reduction strategy in vehicle design is often associated to the massive use of composite materials

like carbon fiber-reinforced plastics (CFRP) that allow us to make lightweight structures [12]. Thanks to CFRP, the designers have at their disposal two specific tactics for structural optimization [13,14], directly linked to the mechanic properties of composites and able to permit lighter structures: (1) high stiffness to weight ratio; and (2) strong anisotropy in the mechanical properties.

The possibility to design the material properties in different directions leads to a greater designers' opportunity in optimization. However, as mentioned, this optimization may involve multiple aspects, not always related to the same category. Precisely for this reason, the general concept of multi-objective optimization is spreading in vehicle engineering [15].

The scientific term of 'optimization' refers to that branch of applied mathematics which studies theory and methods for the research of the maximum and minimum points of a mathematical function [16]. The scope is to translate into mathematical terms a given problem which can be related to different disciplines, as physics, finance, social sciences and, of course, engineering [17–28].

In [18] and [19], for instance, the structural improvements in vehicle components using a design procedure not more based on a trial-and-error process, but linked to optimization methods are discussed, showing their potential, in the case of weight reduction in automotive chassis design.

Far from intending to represent all the general aspects of the optimization problem or its methods and processes in engineering (in part available in [24]), it is possible to here summarize some fundamental elements of specific interest for the present study, as:


Although widely used in engineering (as for structural design, product design, shape optimization, topology optimization, processing planning, and so on), the complexity of applying optimization methods as generally developed by mathematicians has often convinced vehicle designers to additionally consider hybrid and/or simplified approaches.

In [25], for instance, optimized solutions in structural design were achieved merging the Taguchi's method for robust design and the particle swarm algorithm for optimization.

The general state of the art in the field of composite laminates and sandwich optimization is available in two recent comprehensive reviews [26,27]. Results from analyses carried out since 2000 are discussed based on the type of structures, objective functions, design variables, constraints and applied algorithms. Addition factors as boundary conditions, orientation of fibers, design variables are also considered respect to improvements in mechanical behavior such as buckling resistance, stiffness and strength along with reducing weight, cost and stress under various types of loadings.

Closer to the current case, in [28] an aerodynamic shape optimization with the scope to reduce the drag coefficient is performed based on a hybrid process that coupled a genetic algorithm with a quite common iterative method for solving unconstrained nonlinear optimization problems. As often happens, this investigation merged potentialities of a commercial code (a Computational Fluid Dynamics code in the case) with optimization methods and validation test cases.

Furthermore, in [29] a simpler methodological approach is proposed and used to develop an automotive structure able to balance technological, economic and ecological aspects at the same time. It is focused on the application of foams as a core material in sandwiches for floor panels in a concept car. The problem of multi-objective optimization is solved by a procedure of weight minimization that

also ponders the structure behavior respect to static loads in terms of stiffness, strength and buckling constraints. The comparison is formulated for each material application, including considerations on mass and materials' price, and environmental impact.

Strictly in line with that last two cases, the current work was based on a multi-objective optimization where a compound objective function was adopted in the form of a sum of performance indexes balanced by proper weighting coefficients. The specific scope of the paper and its novelty are to show the development of a practical real-life application from the stage of the conceptual design till the actual embodiment and the manufacture of the component. In this framework, the multi-objective optimization is constrained from the beginning by geometrical and functional requirements, moreover also the raw material and the manufacturing technology cannot be changed.

#### **2. Materials and Methods**

#### *2.1. Composite Materials and Structures*

This study aimed to redesign a new photovoltaic roof for the solar vehicle, designed and built by the University of Bologna to take part in solar competitions worldwide [30].

The vehicle is a four-passenger quadricycle [31] where every technical solution was designed with the scope to improve the overall car performance (Figure 1). The replacement of the out-of-date photovoltaic panels, as well as changes in the safety structure eliminating the metal roll bar [32], driven by the need to respect rules of a different race, gave the chance to intervene on the roof through this optimization.

**Figure 1.** The multi-occupant solar vehicle: (**a**) general layout; (**b**) roof support structure.

As a racing car, each part of the vehicle should be as light weight as possible, but also with excellent mechanical properties, permitted by the large employment of high stiffness and strength Carbon fiber reinforced plastics (CFRP). The same occurred for the roof, which is the structure of the photovoltaic panel, where unidirectional (UD) and bidirectional (twill) CFRP were layered in the form of sandwich structured laminates. Specifically, Toray T1000 UD and T800 twill pre-impregnated fabrics with, respectively, 0.15 mm and 0.30 mm thicknesses, were chosen considering their remarkable mechanical properties. Polyvinyl chloride (PVC) foam is used as sandwich core. The outermost layer is a T800 fabric, with fibers along 0◦/90◦ direction. Internally there are two unidirectional layers in T1000, at 0◦. The PVC core is not present everywhere, but only in the central section, having a roof with variable thickness, between 1.0 and 10.0 mm (and 4.0 mm in this case). Table 1 reports the main mechanical properties of materials and Table 2 their layout. It is evident the marked anisotropy (unidirectional or bidirectional) in properties of composites, but also a symmetrical layout respect to the core.

With reference to the direction of the reinforcing fibers, it should be noted that the angle in the table is to be considered with respect to the construction method (better specified below). In grid structures the main directions (0◦) are along the direction of the different beams. In the perforated panels there is a single main direction that coincides with the longitudinal direction of the roof.




#### **Table 2.** Composite layout.

#### *2.2. Optimization Parameters*

In general terms, the design optimization of the roof in the case of this solar vehicle was based on a combination of mechanical and functional features, as reported in Table 3. The first ones, i.e., the mechanical features, are related to the structural stiffness respect to static and dynamic loads and are required to assure the best conditions in terms of integrity, safety and modal response. They can be measured by flexural stiffness, torsional stiffness and (first) resonance frequency in the way that the higher their values, the better for the design solution is. The second ones, i.e., the functional features, are also desired to assure the best performances in terms of energy efficiency, including the reduction of inertial masses and possibility of maximizing the area exposed to heat exchange.


**Table 3.** Optimization objectives and related weighting coefficients.

At the same time, since these features are clearly interconnected, it was preferred to reduce the size of the system by fixing one parameter, the weight (= 0.250 ± 0.010 kg) as common design target, and performing the optimization based on the other factors. Since each factor has to be maximized, it was possible to express the performance index, the index for comparison (IC) as:

$$\text{ICC} = \frac{1}{\sum\_{i=1}^{n} p\_i} \sum\_{i=1}^{n} p\_i \left( \frac{y\_i - y\_{i\text{ max}}}{y\_{i\text{ min}} - y\_{i\text{ max}}} \right) \tag{1}$$

where *pi* represents the weighting coefficients (with *i* = 1 ... 4 and Σ *pi* = 1) and *yi* is the *i*-th optimization parameter which can vary between *yimin* and *yimax*, representing the minimum and maximum value assumed by each parameter *yi*.

#### *2.3. Overall Methodology*

The present investigation was also based on the following concepts and phases (Appendix A):

	- A composite laminate made in a single piece where rectangular holes, with rounded edges (a), or elliptical ones (b) where shaped through;
	- A grid-based structure, made up of a series of intersecting straight (vertical, horizontal and angular) lines (grid lines), forming a rectangular (c) or triangular (d) texture of beams.

**Figure 2.** Shape definition: four geometries were considered: (**a**) rounded rectangles and (**b**) ellipses holes, or (**c**) rectangular and (**d**) triangles in quadridirectional grids.

It is noteworthy that, actually, Table 4 reports only a subset of parameters suitable to describe in their generality the geometric properties of the represented shapes (i.e., rounded rectangles, ellipses, rectangular and triangles) especially as regards the grid configurations since the possibility to overlap beams in a large variety of combinations (angles, distance between elements' recurrences, etc.).

**Table 4.** Geometrical parameters for optimization.

#### *2.4. Optimization Space Reducing*

However, these options already create a quite large *n*-dimensional optimization space (with n ≥ 10 based on how certain geometric similarities are taken into account) which must be deemed by four optimization objectives (output), or even five, if the weight is reintroduced in the comparison. Each point of this *n*-dimensional space should represent a potential solution of the optimization problem, that deserves to be investigated (in terms of stiffness, resonances and so on) by three separate numerical analyses (one modal and two static). At the same time, it is evident it is not possible to reiterate these calculations for an amplified number of situations. Then, as commonly happens in each optimization analysis, a way to decrease the size of the solution space was achieved by:


**Figure 3.** Different scales of analysis: basic unit (500 × 500 m), experimental mock-up (2500 × 800 mm) and half (front) section of the roof (~5200 × 1.600 mm).

#### *2.5. Design and Simulation Tools*

3D CAD modelling was performed in SolidWorks (by Dassault Systèmes, Vélizy-Villacoublay, France). Ansys Workbench ver. 18.1 platform (by ANSYS, Inc., Canonsburg, Pennsylvania, USA) was used for structural analysis of composites (by an ACP toolkit) and for design optimization (Design Exploration) by the direct optimization features and the screening and response surface method.

Specifically, in 3D modelling solid parts were suppressed with the scope to manage 3D surfaces. Discretization was done by shell elements (FE), quadrangular of quadratic order with 8 knots, preferred for precision. After mesh convergence tests, the maximum size for FE was set at 3 mm (in the case of basic units), employing 12,000–18,000 FE for meshing the different configurations.

In accordance with reality, the thickness of sandwiches was made variable along the section by creating cut-off selection rules in ACP and assigning them to the PVC core. The same technique was adopted to manage edges and intersections. A visual example of effects of expedients used for the correct discretisation of composite structures is reported in Figure 4 in the case of rectangular grid. The 500 × 500 mm base section, consisting of 4 × 4 squares (representing the solar cell frames) are shown highlighting differences in thickness and the final mesh with implemented these changes.

**Figure 4.** Details in sandwich discretization (in the case of rectangular grid): (**a**) the 500 × 500 mm base section; (**b**) vertical beams over respect to the horizontal ones; (**c**) sides where the thickness of the sandwich must decrease; (**d**) final mesh (also showing the thickness changes).

The so-called single-layered shell method [33], instead of many others (e.g., stacked shell [34]), was chosen and implemented for discretizing the layout. This simplification quite common in analyzing composites in the case of large and complex structures, permits to quickly investigate the main phenomena at the level of macroscale. With a proper conversion in properties, it reduces a multilayered laminate to an equivalent single layer laminate and use shell elements all along the surface with integration points (IP) throughout the thickness [33]. Specifically, in the case an IP was set per each of the seven layers: no additional IP was used to investigate the six interfaces between layers. In such a way it was possible to drastically limit the number of FE speeding up the simulation focusing the attention on in-plane phenomena [33–35].

However, this methodological limit has no practical effect on the present study since the roof has no global structural functions [35]. The vehicle was designed around a lower monocoque and an upper structure (Figure 1b). A rigid frame on the composite structure guarantees the fixing and solidity of the roof which, practically, in addition to its own weight, must only support the solar cells (~ 1.2–1.5 kg/m2).

From the Ansys internal library, material models labelled as "Epoxy Carbon UD (230 GPa) Prepreg", "Epoxy Carbon Woven (230 GPa) Prepreg" and "PVC Foam (80 kg m<sup>−</sup>3)" were chosen for, T1000, T800 and PVC respectively, but default properties were changed in accordance to Table 1. ACP toolkit permitted to build exact layouts as shown with Table 2. In the grid configuration, in order to correctly orient the fibers according to the direction of the single beam, it was necessary to divide the grid in section subgrouping beams characterized by a same orientation in the fibers.

#### *2.6. Simulation Procedures*

Three different simulations were carried out on each configuration under investigation.

#### 2.6.1. Resonance Frequency

The first simulation was a modal analysis and concerned the resonance frequency. It was carried out in free-free conditions: no constraints or forcing were set. The first six modes are rigid modes, corresponding to the structure's six degrees of freedom. Therefore, the seventh mode and frequency correspond to the first modal shape of interest (Figure 5a).

**Figure 5.** Numerical evaluation of (**a**) resonance frequency, (**b**) flexural and (**c**) torsional stiffness.

#### 2.6.2. Flexural Stiffness

The second simulation was a static analysis and concerned the flexural stiffness. Load conditions were set equivalent to a simply supported beam with distributed load. Constraints were assigned to two opposite external sides of the beam and a force equal to 1 N (to simplify calculations) to upper face. The model was simplified as a beam considering the complete roof structure, consisting of various repetitions of it, has a prevalent dimension over the others. The required output is the flexural displacement (total deformation). Flexural stiffness can be calculated by dividing the force by the maximum displacement (Figure 5b).

#### 2.6.3. Torsional Stiffness

The third simulation was a static analysis, as the previous one, but concerned torsional stiffness. As a system, it was considered a cantilever beam with a pure twisting moment applied at the end. A fixed support was therefore assigned to one side and a torque of 1 N·mm to the opposite side. Torsional stiffness can be calculated as:

$$\begin{array}{c} K = \frac{T}{D} \\ \theta = \tan^{-1} \left( \frac{z\_1 + z\_2}{2b} \right) \end{array}$$

where *T* is the twisting moment, θ is the rotation, *z*<sup>1</sup> and *z*<sup>2</sup> are the vertical displacements at the two ends and *b* is the length of the beam. The rotation θ was determined, in the case, directly using the Flexible Rotation Probe function instead of the analytical formula after verified the results by the two methods were similar (Figure 5c).

#### **3. Results and Discussion**

#### *3.1. Shape Topological Optimisation*

In Table 5 a first comparison between geometries at the level of basic unit is reported in terms of performance index (IC) respect to their main dimensional characteristics (sizes).


**Table 5.** Effect of changes in geometry (at the level of basic element).

Specifically, in the case of laminates dimensions refers to the hole sizes as width (W), height (H) and fillet radius (r) for a rounded rectangle or 1st and 2nd (A, B) axes for ellipsis.

In the case of the grids, the dimensions refer to the width of beams that make up the structure: an upper orthogonal (0–90◦) grid, reinforced by a lower diagonal (±45◦) one. A value of 30mm in width for the upper grid is fixed considering the need to sustain the singular solar cell. Both for laminates and grids 5.2 mm of thickness is fixed since the predetermined composite layout (Table 2).

In terms of results, it should be noted that, together with the torsional and flexural stiffness, the first (lowest) natural frequency is also reported in the table. This value represents an indication of how sensitive the structure is to the effects of vibrations. The maximum displacement respect to the first natural frequency should be also considered representing the effect of vibrations in terms of intensity (not sensibility).

By comparing the values of the performance index (IC), it is immediately evident that the best configurations are, in general, grid configurations (as Figure 2c,d) reinforced by several crosses (IC ≥ 0.48) and, specifically, a quadridirectional configuration (IC = 0.79). Despite this, the option of manufacturing a laminate in a single body and then lightening it through holes could be taken into consideration provided the choice of rectangular (IC = 0.61), almost squared (~95 mm vs. 93 mm), and rather large holes (~8800 mm2) with rounded borders (~24mm). In making these considerations, it should be also remembered that all these configurations are obtained keeping the same weight (w = 0.264 ± 0.01 kg) as target and, consequently, the same amount of material.

Outcomes from Table 5 also take into count, as a preliminary optimization (improved below), the way an upper orthogonal grid (orthogrid) can be reinforced by lower transversal crosses leading to double grid configurations. The most interesting options, between many others under consideration, are displayed in Table 6 where at the investigation level of the basic unit, the impact of a crescent number of crosses is analyzed. The same table also exhibits the effect of changes in cross angles.


**Table 6.** Effect of changes in the crosses (at the level of basic unit).

Results, as said, assumed that:


For each configuration, the table reports stiffness in the cases of application of (flexural or torsional) loads along X or Y axes. Differences are evident related the geometrical anisotropies of grids respect to these changes in the axes.

#### *3.2. Grid Topological Optimisation*

In Table 7, a second comparison between configurations is reported. It also deals with the best solutions detected in terms of geometry and shapes at the level of basic units (Table 5) but extends results on the mock-up case. In particular, the laminate with rectangular holes where compared with quadridirectional grids in different configurations. Having demonstrated a certain appropriateness (measured through IC) for a specific configuration at the level of basic unit does not necessarily mean, in fact, that this rank remains the same on the mock-up. Such a fact immediately emerges from the table making clear that the laminate with rectangular holes (IC = 0.25) is no longer convenient.

Having also demonstrated the convenience of using, as general design concept, a grid structure where an upper orthogonal grid (on which the solar panels are laid) are reinforced by a series of diagonal crosses, a different topological optimization study was carried out with respect to the definition of the (lower) grid in terms of number of elements and their arrangement.

Without entering in unnecessary details, Table 8 exhibits the impact of diagonal crosses at the level of mock-up, both as number, distribution and crosses angles.


**Table 7.** Effect of changes in grid shape (at the level of mock-up).

**Table 8.** Effect of changes in grid topology (at the level of mock-up).


The beam width in upper grids remains 30 mm (to conveniently host the solar cells), but changes affect the lower ones (from 14 to 50 mm) according to the condition of weight conservation. Spacing parameter (in Tables 8 and A2) provides information on the distance between beams on the lower grid: reducing the distance increases the number of beams (and crosses as reinforce) but, on the other side, as said, reduce their width saving the total weight. Finally, unlike the first analyses on the geometric parameters, where an automatic change for (almost) all the parameters was possible, in these further investigation configurations often involve the development of specific CAD models. Thus, several actions were done manually.

#### *3.3. Multi-Objective Optimization Results at a Glace*

In line with results from similar investigations (as [36]), it is clear how the grid structures are preferable to perforated panels. It is also evident that a grid with a denser mesh is better than those of sparse one with an equivalent weight. However, although the comparisons give unambiguous indications when only mechanical properties are considered, additional criteria should be included.

In fact, the optimized option as here detected, i.e., quadridirectional grid with crosses every 250 mm, could present (compared to less dense grids) additional challenges in its applicability such as:


These aspects were not including in the optimization since the beginning and are not considered.

A preliminary estimation of the free surface carried out by Ansys with respect to two borderline grids showed a quite low variability (<10%): using a denser grid does not significantly affect the panels cooling. It also depends to the fact that, denser grids are necessarily made by tighter beams (since the precondition of equivalent weight).

Regarding the composite roof construction, with an area of approx. 8.0 m2, its manufacturing can represent a laborious task. A larger number of beams and crosses, to be made and glued, would be preferable only when mechanical properties were significantly better: it is not the present case.

Table 9 reports an update in design solutions evaluation, taking into count of additional parameters and objectives of optimization. In particular, the maximum displacement when the lowest natural frequency occurs is included: those parameters have to consider in combination for better analyze the dynamic behavior of the structure. Furthermore, the table also introduces parameters, not strictly related to mechanical properties, as the area available for heat transfer and an indicator of producibility. While the exposed area may be directly detected by Ansys functionalities, an empirical index had to be defined for estimating the producibility, it was done in accordance with information from manufacturer, making its value proportional to the number of operations necessary to build composite structures.


**Table 9.** Multi-objective optimization of composite structure.

As a synthesis, a quadridirectional grid with crosses every 500 mm was adopted. In accordance with previous results, in fact, this structure exhibited good mechanical properties, better than both perforated sandwich panels (especially regarding the resonance frequency [37]), but also compared to an orthogonal grid, thanks to the presence of reinforcing beams in the diagonal directions [38]. However, its medium sparse grid does not entail the construction problems that can occur during production in the case of denser grids.

Finally, since it was noticed that a slight increase in the mass of the angular grid provided positive benefits to the quadridirectional grids, it was also preliminary checked the effect of minor changes in the constraint of the minimum width for the orthogonal beams (=30 mm).

The basic unit of 500 mm sides was examined, performing a further optimization in Ansys by varying the width of the orthogonal beams between 25 mm and 30 mm and diagonal beams between 30 mm and 50 mm. Since the condition of equal total weight, an increase in the mass of one grid, orthogonal or diagonal, decreases the other. In general terms, the best situation is present when widths

of both grids are quite similar. A specific optimization can be obtained when the choice of the solar cells, defined the lower value of width for elements in the upper grid.

#### *3.4. Validation*

A modal analysis was used for validating the FE model by experiments.

The mesh consisted of 56163 shell elements and allowed to identify the first 20 modes of the structure (of which, six are rigid modes). Deformations along main directions and Cartesian axes were also evaluated. On the other side, a modal test was also carried out both with accelerometers (as in [39–41]) and Bragg grating fibers (FGB) sensor (as done in [42]). Figure 6 and Table 10 summarize this comparison for a range lower than 100 Hz, equivalent to the five lowest natural frequencies. A good accuracy is clear with a 7.3% average deviation between predictions and experiments. Besides, a constant underestimation is also evident which suggests a refining in the FE model discretization.

**Figure 6.** Modal analysis with results from (**a**) experiment and (**b**) simulation.


**Table 10.** Comparing the results from modal analysis in terms of natural frequencies.

#### *3.5. Results Implementation*

The same design approach was used to finalize the roof design. The 3D geometry of the roof was filled by structural modules with design parameters based on the previous optimization. Specifically, a quadridirectional grid with a medium density of beam element was used as valid compromise between the different targets. It was characterized by an orthogonal grid on the top and a lower diagonal grid on the bottom with, respectively, 300 mm and 200 mm widths. This design pattern was replicated along the roof dimension and shape, partitioned in front and back sections (Table 11).

**Table 11.** Overall geometrical dimensions for roof sections.


The dynamic behavior of these parts was examined by a further modal analysis considering unloaded bodies (as previously) but constrains as in real case. First frequencies were 21 and 12 Hz with total displacements of 42 and 69 mm for, respectively, front and back sections (Figure 7)

**Figure 7.** Numerical modal analysis of (**a**) front, (**b**) back section of the roof.

Then, the composite structure was manufactured using manual layout and autoclave molding and hot gluing with two-component for part assembly. The solar cells, E60 bin Me1 by Sunpower, were directly laminated on the panel with ethylene vinyl acetate (EVA) films. Flexible layers on both surfaces (front and back) for a total of five layers (including 2 EVAs), for a 1.5 mm of overall thickness, were overlapped and cured in autoclave. In particular, solar cells were positioned with the scope to maximize the energy yield of vehicle-integrated photovoltaics (VIPV) [43].

Lastly, the solar roof was installed and finalized with other vehicle components. Figure 8 shows several imagines from the solar vehicle at the end of the present investigation. Functional tests were performed running over 800km on roads open to traffic. No complications emerged related during this first trial in terms of structural design of the roof which appeared stable and functional.

**Figure 8.** Solar roof manufacturing and installation on the solar vehicle.

#### *3.6. Further Considerations and Novelty*

Following the use of a multi-objective optimization in this design action, several general considerations can be introduced. As first, it is essential to note how the outcomes strongly depends on the weight assigned to the objectives: even minor changes in their values can lead to very different design indications. Thus, it would be relevant to find criteria for an empirical definition of these parameters. By an analysis simply based on mechanical aspects (e.g., Eigen value optimization), it is quite hard, in fact, to detect very uncommon design solutions or issues respect what already available as technical know-how such as, for instance, the superiority of reticular structures. Moreover, when further objectives are introduced in optimization, but always of the same type (such as, e.g., the minimization of an area keeping weight or stiffness unchanged), results do not change. Reticular structures seem the best solutions, especially when characterized by a certain geometrical complexity. However, as soon as a goal not strictly related to the 'structural engineering' is introduced, optimal solutions start to evolve along unpredicted directions.

For instance, a sandwich panel, made in composite by a single stratification that also considers the lighting holes, can be much simpler in terms of producibility respect to a composite grid. Then, when an additional parameter related to the producibility is introduced in the evaluation, moving the analysis from a 'topologic' to 'multi-object' approach, this perforated panel significantly increases its ranking respect to grid solutions. Similarly, lower density grids start to be more attractive.

However, without the possibility of objectively validating the weight of each objective, a deeper level of analysis risks to be inconsistent respect to the real applicability of results, which represents the essence and novelty of the present investigation. In fact, this study was intended to be a first attempt in the contest of solar vehicles to move from redesign action based, as tradition, on a trial and error approach toward a multi-step and -objective one method. To the knowledge of the authors, in fact, no other design studies are available dealing with the multi-objective optimization in an automotive context that involve fiber-reinforced sandwich structures of such large size and geometrical complexity.

#### **4. Conclusions**

The application of a multi-objective approach to the complete design process of an actual composite structure was shown. Specifically, a photovoltaic roof for solar vehicles was designed following a multi-objective approach in the way to balance divergence structural criteria, as static stiffness and dynamic response, with additional functional targets. Heat transmission and energy efficiency were also considered.

As first, several alternative shapes (circles, squares, triangles) were compared as pattern. Then, for selected shapes, the optimization was carried out respect to their main parameters (e.g., lengths, angles) searching the optimal points inside Ansys FEA software. A performance index was properly defined to represent the best compromise and a large number of configurations were compared. This index combined aspects as flexural stiffness, torsional stiffness, resonance frequency and heat transfer surface by the definition of weight parameters for each target. To simplify the study, a multistage approach was preferred. As first, a 500 × 500 mm section, equivalent to 4 × 4 solar cells, was adopted as base for an initial comparison between fundamental shapes (e.g., triangle, rectangle, ellipses). A total of 49 designs, each one characterized by a specific combination of shape and geometrical parameters (as widths, angles) were considered and compared, limiting the input for the next stage to few (4) optimal options. Hence, this base unit was used as design modulus to build a larger geometry (2500 × 800 mm) able to better predict and compare the structural behavior in a case study closer to reality. This second step of optimization, performed by 51 designs, permitted to recognize the quadridirectional grid as the best solution and to define proper combinations of geometrical parameters.

The FE model was validated comparing results from simulations and experiment respect to a modal analysis: a good accuracy with a 7% deviation in predicting the lowest natural frequencies was detected. The structure design, as here optimized thanks to a quadridirectional grid, was applied to the case of the real roof, characterized by larger dimensions (~5200 × 1.600 mm) and a double curvature geometry. A FE modal analysis of the roof, done in accordance with real loads and constrains, was carried out to determine the lower frequencies (higher than 12 Hz) and modes.

The composite structure was produced using autoclave technology; solar cells were also direct laminated on it. The solar panels, with an overall thickness lower than 5.2 mm, were installed on the vehicle and functionally tested on the road with valid results.

**Author Contributions:** Conceptualization, C.F. and G.M.; methodology, C.F. and G.M.; software, D.S. and A.P.; validation, A.P. and G.M.; formal analysis, A.P.; investigation, D.S. and A.P.; resources, G.M.; data curation, D.S. and C.F.; writing—original draft preparation, C.F. and D.S.; writing—review and editing, G.M. and C.F.; visualization, C.F.; supervision, G.M.; project administration, G.M.; funding acquisition, C.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research has been carried out within the international collaboration project 'Two Seats for a Solar Car', an intervention funded by the Italian Ministry of Foreign Affairs and International Cooperation (MAECI), aimed at transforming a racing solar car into a solar road vehicle. The experimental activity was funded by the Italian Ministry of University and Research (MIUR) within the project PRIN2015 "Smart Composite Laminates".

**Acknowledgments:** The construction and installation of the composite roof were carried out by members of Onda Solare Sports Association at the MetalTig Srl workshops (in Castel San Pietro Terme, Italy). A special recognition for their commitment is dedicated to Mauro Sassatelli, Morena Falcone, Marco Berdoldi and Luigi Russi. The authors also thank to Marco Troncossi, Francesco Falcetelli and Alberto Martini for their support during the experimental sessions.

**Conflicts of Interest:** The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A**

This section contains details on the optimization procedure and parameters. The optimization was based on a three-step process where, as first, a selection of shapes was compared and optimized by 49 designs on a basic unit (500 × 500) (Tables 5 and 6 and Table A1), then, the comparison between geometries was improved by 51 designs at a level of mock-up (2500 × 800) (Tables 7 and 8 and Table A2), finally, the most promising solutions were evaluated in terms of multi-objective targets (Table 9). During the optimization, parameters were considered as fixed or modified inside the procedure.


**Table A1.** Preliminary (hole) shape optimization on basic unit.

#### **Table A2.** Grid optimization on the mock-up.


#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18