Sage Journals: Discover world-class research

Abstract

In this paper, a novel multiscale and multi-stage structural design optimization procedure is developed for the weight minimization of hopper cars. The procedure is tested under various loading conditions according to guidelines established by regulatory bodies, as well as a novel load case that considers fluid-structure interaction by means of explicit finite elements employing Smoothed Particle Hydrodynamics. The first stage in the design procedure involves topology optimization whereby optimal beam locations are determined within the design space of the hopper car wall structure. This is followed by cross-sectional sizing of the frame to concentrate mass in critical regions of the hopper car. In the second stage, hexagonal honeycomb sandwich panels are considered in lower load regions, and are optimized by means of Multiscale Design Optimization (MSDO). The MSDO drew upon the Kreisselmeier–Steinhausser equations to calculate a penalized cost function for the mass and compliance of a hopper car Finite Element Model (FEM) at the mesoscale. For each iteration in the MSDO, the FEM was updated with homogenized sandwich composite properties according to four design variables of interest at the microscale. A cost penalty is summed with the base cost by comparing results of the FEM with the imposed constraints. Efficacy of the novel design methodology is compared according to a baseline design employing conventional materials. By invoking the proposed methodology in a case study, it is demonstrated that a mass savings as high as 16.36% can be yielded for a single hopper car, which translates into a reduction in greenhouse gas emissions of 13.09% per car based on available literature.

Keywords

Freight rail lightweight design multiscale design optimization sandwich panels topology Genetic Algorithm (GA)stochastic optimization Fluid-Structure Interaction (FSI)

Introduction

The increasing need to transport and deliver greater quantities of key economic resources at lower cost coupled with proliferation of environmental concerns has induced greater pressure on freight companies to reduce costs and emissions. It has been demonstrated that for an 10% reduction in vehicle mass, 5% energy savings can be attained, which can address these concerns simultaneously.^1,2 As an added benefit, it is further found that structural weight reduction in rail vehicles can reduce yearly energy consumption while simultaneously improving life-cycle metrics.³ While many efforts have been made to reduce car-body weight in passenger rail, freight rail remains scarcely investigated. In commercial freight rail, reduced weight of car-bodies can potentially result in a greater yield of transported cargo, which can significantly reduce costs or bolster profits.

Implementation of honeycomb sandwich panels has found much success in passenger rail. Magnesium alloy corrugated sandwich panel extrusions were proposed by Lee et al., where extrusion dimensions of each panel was optimized a priori.⁴ Seo et al. developed a hybrid composite bodyshell for the Korean Tilting Train Express, which met all safety requirements and reduced car-body mass by 30%.⁵ Cho et al. employed CFRP-Aluminum honeycomb sandwich composites to the underframe and roof of a K-EMU car-body, resulting in a weight reduction as high as 29%.⁶ Using a multiscale approach where material data for composite sandwich structures was characterized at the coupon level and then gradually extrapolated to the component level, Zinno et al. developed a novel lightweight roof for a rail car, which abides by European Code.⁷

To design for minimal mass, it is often advantageous to employ numerical optimization methodology. In this regard, MSDO has garnered much attention, and is used to bolster capabilities of numerical optimization. Catapano and Montemurro optimized sandwich panels for strength against buckling by integrating a FEM directly into the MSDO loop.⁸ Wennberg et al. employed a topology optimization method to design a novel car-body shell at the mesoscale, and then optimized a sandwich panel to replace a characteristic section of the car.^9,10 Molavitabrizi and Laliberte investigated the multiscale design optimization of hopper cars using octet core lattices optimized at the microstructural level and then extrapolated to the mesoscale.¹¹ Lee et al. employed a multiscale approach whereby topology optimization was used to determine where corrugation is to be concentrated in a rail car and then subsequently those panels were sized.¹² Harte et al. employed a multiscale approach to determine ideal CFRP shape and layer sizing of a passenger LRV.¹³

This paper aims to develop a novel methodology whereby optimization of the structure at the microscale and its implementation toward the mesoscale are coupled in a MSDO procedure to ensure a successful global design. The proposed methodology employs a multi-stage study where, upon determination of the loads, a baseline structure is produced using conventional materials. Upon establishing the baseline, the novel MSDO is used to design a freight car-body of minimal mass employing honeycomb sandwich panels. While much of the available literature discusses static loading, a broader range of load cases is considered in this paper, which implements a dynamic load case that is modeled using Fluid-Structure Interaction (FSI). The Discrete Element Method (DEM) has gained prominence in FSI modeling of structures carrying a granular medium.^14,15 However, until recently DEM remained computationally costly to implement.^16,17 On the other hand, SPH is already used in prior work of Ushakov et al. to model granular cargo, and is known to be successful in implementation.¹⁸ Accordingly, SPH is implemented in this paper to determine cargo dynamic load on hopper car structure, as it is available in a commercial form compatible with techniques for reduced solve times such as mass scaling.¹⁹

This paper is organized into five body sections. The first section is an analysis of the available literature on sandwich panels, whereby the parametric studies from established literature are employed to determine the structural candidate of superior strength-to-weight characteristics. It is demonstrated in the review that hexagonal honeycomb sandwich panels are the candidate with highest potential for weight minimization when implemented in the roof and sidewalls of a hopper car under expected loads. In the second section, the load on a hopper car is characterized in further detail using explicit Finite Element Analysis (FEA) employing SPH to simulate sloshing of cargo under a dynamic load. Further load cases are introduced per regulatory standards in the implicit model as well. In the third section, a multi-stage optimization is introduced, where the hopper car is first optimized globally using topology optimization under all load cases of concern, which provides a car frame design which has been globally optimized for a high stiffness-to-weight ratio. Then, the second stage is introduced, whereby a MSDO is performed on a characteristic panel of the hopper car. In the fourth section, a case study employs the multi-stage design optimization for a generic hopper car. Results of the case study are discussed in the fifth section, and broader implications are summarized in the conclusions.

Analysis of literature and selection of structural candidate

This section aims to perform a literature analysis to determine the sandwich core topology of highest strength-to-weight capacity under expected hopper car loads. Rigidity per unit width of a sandwich panel with very thin faces is given by the equation²⁰

D = \frac{E_{f} t_{f} {(H_{c} + t_{f})}^{2}}{2}

(1)

where D is flexural rigidity of the panel, E_f is elastic modulus of the face, t_f is thickness of the face sheets and H_c is core thickness. Though design of a panel for minimum mass and maximum rigidity might produce potentially optimal results, the strength of the panels must be considered. Failure of hexagonal honeycomb sandwich panels in particular were studied by Petras and Sutcliffe, who expressed failure modes of the hexagonal honeycomb core sandwich panel as²¹

\begin{array}{l} σ_{f Y} & = σ_{Y} & (FaceYielding) \end{array}

(2a)

\begin{array}{l} σ_{i c b} & = \frac{E_{f}}{1 - ν_{f}^{2}} {(\frac{t_{f}}{s})}^{2} & (Intracellular Buckling) \end{array}

(2b)

\begin{array}{l} σ_{w r} & = 0.33 {(\frac{E_{c} t_{f}}{E_{f} H_{c}})}^{\frac{1}{2}} E_{f} (FaceWrinkling) \end{array}

(2c)

\begin{array}{l} σ_{c r} & = \frac{{(H_{c} + t_{f})}^{2}}{2 t_{f} H_{c}} \sqrt{G_{c, x z} G_{c, y z}} & (Crimping - Compressive) \end{array}

(2d)

\begin{array}{l} τ_{c y} & = 1.7 E_{f} {(\frac{ρ_{c}}{ρ_{f}})}^{3} & (Crimping - Shear) \end{array}

(2e)

where s is diameter of the largest possible circle in the hexagonal cell, t_f is face thickness, H_c is core thickness, ρ is density, E is elastic modulus, and G is shear modulus. Subscripts c and f denote core and face properties, respectively, while subscripts x, y, and z denote the two in-plane directions and the out-of-plane direction, respectively. In this paper, steel is assumed for all structural components. Therefore, the face elastic modulus of E_f is also used to denote the properties of solid steel.

Equations (2d) and (2e) reference two forms of crimping, namely, compressive and shear. Though the failure modes are similar, shear crimping occurs due to out-of-plane shearing of the core, while compressive shear occurs when buckling of the panel results in earlier failure of the core than the face, resulting in out-of-plane shearing. Both forms of failure are investigated in this paper to ensure a sandwich panel which is structurally sound in all load cases.

Rathbun et al. investigated strength-to-weight ratios of various metallic sandwich panels of varying topologies in a parametric study using the failure characteristics laid out above.²² To characterize the load experienced by sandwich panels under bending in a parametric manner, Rathbun et al. introduce the load index, Π, of a sandwich panel as^22,23

Π = \frac{V^{2}}{E M}

(3)

where M is maximum bending moment, V is maximum shear, and E is material elastic modulus of the sandwich material. Equations for weight indices are also derived for various core topologies and vary depending on core geometry and load index regime. Derivations of all but the hexagonal honeycomb core are found in the work of Rathbun et al., while Al-Sukhon derives that of the hexagonal honeycomb panel in their thesis.^22,24 Unitless load indices derived according to the Association of American Railroads’ Manual of Standards and Recommended Practices (AAR MSRP) are now introduced for the roof, wall, and floor of a hopper car as^22,23

\begin{array}{l} {(Π)}_{roof} & = 2.5 \frac{2 {(w_{o})}_{roof}}{E_{f}} \end{array}

(4a)

\begin{array}{l} {(Π)}_{wall} & = 1.25 \frac{W h_{w}}{\sqrt{3} E_{f}} (\frac{1 - \sin (ϕ)}{1 + \sin (ϕ)}) \end{array}

(4b)

\begin{array}{l} {(Π)}_{floor} & = 1.25 \frac{2 W h_{w}}{E_{f}} \end{array}

(4c)

where

{(w_{o})}_{roof}

is pressure applied to the roof of the car, W is specific weight of the cargo, h_w is depth of the cargo and ϕ is static angle of repose of the cargo. Note the equations above are derived considering a simply supported sandwich panel under even distributed loading in the case of the roof and floor. A linearly distributed load decreasing to zero at the top is considered loaded for the walls. The wall load is derived in the Appendix as an example of how the load indices are determined.

Next, Rathbun et al. establish the parametric relationship between the load and weight indices, Ψ, of a metallic sandwich panel.²² Here, the weight index is a unitless property related to the load index based on the geometry of the core and the dominant sandwich failure mode. The development of the weight index is discussed extensively in the work of Rathbun et al.²² and in the Thesis of Al-Sukhon.²⁴ This relationship is plotted in Figure 1. For an example calculation to show the anticipated load index for each component, a wall height of 2.5 m and a cargo specific weight of 850 kg/m³ are used to give numerical value to the load index, while considering an elastic modulus of 210 GPa. From the plot, it is evident that the relationship between weight and load index is discontinuous. This discontinuity is because for each panel, there exists a weight index equation for each panel failure mode. The equation providing the highest weight index is taken to be the panel of conservative design. As the load index varies, so too does the dominant failure mode for each panel, and by extension the equation dictating the weight index. The result is a parametric comparison between the weights of core topologies of different forms of sandwich panels given the same loading and material.

Figure 1.

Weight-to-Load Index of sandwich panels with varying core topologies on a logarithmic scale based on the work of Rathbun et al.²² Load indices of hopper car structural components are highlighted by the dashed vertical lines.

As shown in Figure 1, hexagonal honeycomb sandwich panels demonstrate potential to offer a design of minimum mass in the walls and roofs of hopper cars. Meanwhile, the floor of the hopper car is found to be in the high load regime where the weight index for a sandwich panel approaches that of a solid metallic plate, implying its implementation is not feasible. Therefore, the present study focuses on implementation of hexagonal honeycomb core sandwich panels in the wall and roof of a hopper car.

Numerical methodology

Scope of analysis

Design loads calculations are based on a combination of both explicit and implicit FEM, where explicit FEM is employed for determination of dynamic loading by FSI, and linear-static implicit FEM is leveraged for all other load cases. Four load cases from Section C, Part II – Design, Fabrication, and Construction of Freight Cars, M-1001 of the AAR MSRP are identified and used for analysis in the present work²³:

• Live Vertical Loading (LVL): The LVL is meant to test the structural integrity of the hopper car when fully laden. This load case comprises:

– A floor loading with a pressure P = Wh_w to simulate the cargo

– A side wall lateral hydrostatic load according Equation:

{(w_{o})}_{wall} = 1.25 W h_{w} (\frac{1 - \sin (ϕ)}{1 + \sin (ϕ)})

(5)

– A roof pressure of 720 Pa, and a concentrated load of 1335 N on the middle of the roof.

To simulate inertial loads resulting from a sharp turn, a lateral acceleration of 0.3 g is considered.²³ All loads are multiplied by a safety factor of 1.25 per requirements of the MSRP.

• Draft Tensile Load (DTL): The draft tensile load simulates the hopper car under conditions where excessive pulling might occur. The same loads as the LVL are applied without the factor of safety, along with a tensile load of 2,805,660 N on the coupler.²³

• End Compressive Load (ECL): The compressive load test simulates a situation in which one hopper car crashes into another, resulting in compressive loading on the underframe. Here, it is assumed that the hopper car is unladen and a 4,448,000 N compression load distributed along the underframe is applied.

• Modal Analysis (MA): The MSRP requires that damping be applied for natural flexible body modes with a frequency under 12 Hz.²³ A study of normal modes is considered in this paper to evaluate natural frequencies of the structure throughout the optimization.

Though the above load cases can produce adequate stress results for sizing and design optimization of the car structure, an additional FSI by means of explicit FEM is proposed. In this FSI load case, the entire structure and its contents are excited by a 0.3 g acceleration over a time of 1.85 s to simulate a sharp turn. The time of 1.85 s is selected based on how long it takes for the particles to settle after sloshing. In the LVL, this lateral acceleration can only act on the walls of the hopper car and is not applied to the cargo. Inclusion of the FSI as an additional load case is meant to correct for this by determining the distributed loads on the hopper car resulting from the cargo acting upon it. Details of the evaluation of this load case by means of explicit FEM is given below. Since many of the loads in the LVL are not applied in the FSI, a load factor of 1.5 is used rather than 1.25 in order to ensure conservative results.

Explicit finite element model

The structure of the baseline hopper car is modeled using RADIOSS²⁵ based on a simplified hopper car geometry employing 42,875 quadratic and triangular elements, and 930 beam elements. Due to the costly nature of explicit dynamic analysis, a trade-off between simulation accuracy and solve time must be accepted. The number of elements is chosen such that critical areas of the car’s stress gradients are smooth and continuous while allowing for a large enough time step for reasonable solution time. For the cargo, 145,655 SPH elements are used to model granular material sloshing in the car-body. The Drücker-Prager material law was used by Ushakov et al. to describe the motion of granular cargo under inertial loading.¹⁸ Therefore, the Drücker-Prager material law is implemented in the explicit FEM to simulate granular motion of the cargo.²⁶ Properties of the cargo are shown in Table 1 along with the properties for the car structural steel according to the findings of Attar.²⁷

Table 1.

Material parameters for the steel and SPH cargo.

Material	Parameter	Value	Units
Steel	Elastic modulus, E	210	GPa
	Poisson ratio, ν	0.3	—
	Density, ρ_s	7800	kg/m³
Cargo	Initial density, ρ_cargo	850	kg/m³
	Cargo elasticity, E_cargo	26.813	MPa
	Poisson ratio, ν_cargo	0.25	—
	First material parameter, A₀	3338.5 × 10⁹	Pa²
	Second material parameter, A₁	2075700	Pa
	Third material parameter, A₂	3338.5	—
	Bulk unloading modulus, B_UNL	144	GPa
	Static angle of repose, ϕ	25°	—
	Soil cohesion parameter, c_cargo	1.5 × 10⁶	Pa

The simulation is performed in two stages over 1.85 s. In stage one, the simulation is run with kinetic energy relaxation in order to allow the SPH particles to settle in a rested position while eliminating artificially high deformations due to the settling process. In stage two, a 0.3 g lateral acceleration with a load factor of 1.5 is applied to simulate the sharp turn.

It is shown from Figure 2 that in the first 0.5 s of the simulation, the kinetic energy rapidly spiked and was subsequently relaxed as part of the settling phase. In the sloshing phase from 0.9 s onward, it is shown that the total and internal energies are equal, while kinetic energy is effectively zero, implying that all the energy in the system is translated into strain energy in the structure. This complete transfer of energy demonstrates that the simulation is successful and valid. The results of the FSI are shown in Figure 3.

Figure 2.

Time history plot of simulation energy in the explicit FEM. The energy plot shows that all kinetic energy is transferred to the structure in the form of strain energy, thus verifying the success of the model.

Figure 3.

Results of the explicit FSI model.

Implicit finite element model

From Figure 4, it is shown that after the settling phase, contact is initiated between the SPH cargo and the walls. During the sloshing, the lack of damping results in the excitation of the bending mode of the hopper car wall. This in turn causes the car wall to break contact with the SPH momentarily, then re-initiates contact as it reconnects. While this results in artificial fluctuation, the fluctuation is constant and continuous. Furthermore, application of damping would result in the reduction of the SPH contact load. Therefore, the response at the time of peak contact force is said to provide a conservative estimation of the FSI loading. This time is seen from Figure 4 to be at 1.5 s, thus the response at this time is transferred to an equivalent static model. Note that the peak at 1.5 s is the second of three peaks. This peak in specific is chosen because it is at this point where the SPH cargo is fully in contact with the wall, as opposed to the first peak, where the SPH has just begun to make contact. In taking this peak, it is ensured that the full weight of the SPH cargo is considered.

Figure 4.

Time history plot of total interface contact force of the explicit FEM.

The equivalent model uses the dynamic displacements and accelerations of the time dependent explicit model in the following equation

F_{eq} = a_{d} (M \overset{..}{x_{d}} + K x_{d})

(6)

where F_eq is the equivalent static load vector for each node, M and K are mass and stiffness matrices, respectively, and

{\overset{..}{x}}_{d}

and x_d are nodal displacements and accelerations experienced by the structure, respectively. From Figure 4, the time of maximum loading occurs at 1.5 s, implying that the time of critical loading occurs at this point and is what should be used for characterization of loading. By applying an additional dynamics safety factor of a_d = 1.25, interface forces for all times in the sloshing portion of the analysis predict a higher load than all times without mass scaling. This indicates that any time step in the sloshing phase of the mass scaled simulation with the applied safety likely to provide a conservative load case for an equivalent static model.

The final solve time is found to be 18 s on an Intel(R) Core(TM) i7 CPU with a clock speed of 2.67 GHz. In the implicit equivalent static model, all other load cases mentioned previously are employed in separate sub-cases.

Optimization framework

A multi-stage approach employing a hybridization of various optimization algorithms and strategies is proposed. The first stage considers a topology optimization of the global structure, while the second stage considers a MSDO for the component level.

Stage I—Topology optimization

The first stage in the presented method employs a topology optimization at the mesoscale whereby optimal beam location is determined. Topology optimization is performed using the Solid Isometric Material with Penalization (SIMP) method by Bendsøe.²⁸ A process flow diagram of the procedure is shown in Figure 5. The topology optimization aims to find critical load paths for the placement of the support beam structures. In each iteration, topology optimization for minimization of compliance with a volume fraction constraint is performed thereby generating an element density distribution. This plot is subsequently leveraged to either add or remove beams in critical and noncritical loading regions, respectively. A threshold relative element density of ρ_i ≥ 0.67 is used to filter topology results.

Figure 5.

Topology Optimization process flow diagram.

Upon conclusion of topology optimization, size optimization is performed using Altair OptiStruct’s the Global Response Surface Method (GRSM) for optimization.²⁹ The walls and beams of the structure are divided such that a minimum number of design variables provides a robust optimization while ensuring meaningful convergence. The variables are: a cross-sectional height, H, and width, B, are considered, along with a thickness perpendicular to each axis lying within the cross-sectional plane, T₁ and T₂. For each wall division, a plate thickness of T_p is also considered.

The optimization is constrained such that the von Mises stress in all elements does not exceed 423.3 MPa. Crippling stress, a local buckling phenomenon occurring within individual beams due to compression is also constrained in the sizing procedure. For a square cross-sectional tube, the crippling stress is given as³⁰

σ_{c c} = B {(\frac{g t^{2}}{A})}^{m} \sqrt{E^{m} σ_{c y}^{2 - m}}

(7)

where B = 0.58, m = 0.8, and g = 12 for a square cross section, while E and σ_cy are the elastic modulus and yield stress, respectively. The cross-sectional area, A, is recalculated each iteration according to the design variables. The thickness, t, is the smallest of the two thickness variables T₁ and T₂. Displacement constraints are applied per Table 2 to ensure the resulting structure is equally or more rigid than the original. Detailed results of the GRSM are found in the work of Al-Sukhon.²⁴

Table 2.

Constraints for the post-topology optimization sizing.

Constraint	Upper bound	Lower bound
Shell von Mises stress (MPa)	423.3	—
Beam von Mises stress (MPa)	423.3	—
Beam minimum stress (MPa)	—	Equation (7)
Displacement LVL (cm)	3.1	−3.1
Displacement DTL (cm)	5.0	−5.0
Displacement ECL (cm)	8.1	−8.1
Displacement FSI (cm)	6.5	−6.5
Natural frequency (Hz)	—	12.0

Stage II—Multiscale design optimization

Though traditional topology and sizing methods produce structures of minimal weight rapidly for traditional materials, it was demonstrated earlier in this paper that further savings can be achieved through implementation of sandwich panels in the walls of the hopper car. In the work of Al-Sukhon, a traditional multi-objective optimization methodology is performed to identify the design variables for a hexagonal honeycomb sandwich panel of maximum strength-to-weight ratio, with a section of the hopper car’s the sandwich panel wall modeled using basic statics. The optimization procedure produces a Pareto front which gives the ideal sandwich panel variables. However, it is determined that once these variables are implemented into the FEA of the hopper car, sandwich failure occurs. This is because any analytical and traditional optimization methodology used cannot adequately capture the loading on the hopper car, nor provide results specific to its requirements.²⁴ For this reason, the second stage of the optimization is performed on a critical region of the hopper car using a MSDO with an FEA of the hopper car embedded directly. Use of an MSDO with the FEA in each function evaluation presents two challenges. The first is that there exists no gradient to evaluate the cost function, thus stochastic methods must be employed. In stochastic methods, presence of infeasible design regions in the cost function result in many cost function evaluations since these solvers rely on a prior function evaluation to arrive at an optimal solution. The second challenge is that many of the constraints imposed on the design problem are related to stress in the structure. This means that the constraints can only be evaluated a posteriori; the FEA and cost function must be evaluated to acquire stress results to determine whether the design variables are feasible. Due to this limitation, constraints cannot be directly applied to the system, and must instead be integrated into the cost function in some manner. A penalty method for optimization employing the Kreisselmeier–Steinhausser (KS) equation is used. Though many methods exist for the synthesis of penalty-based cost functions, the KS equations have demonstrated success in prior work, and are therefore used in this paper.^31–33 Through this approach, the constraints are calculated as a secondary term in a penalized cost function, Φ(x), according to the equation

Φ (x) = f (x) + \sum_{m \in ℕ} W_{m} h_{m} (K S (x))

(8)

where m is number of constraints, f(x) is base cost function, and W_m is a penalty coefficient so that the m^th exterior penalty function,

h (KS (x))

, is appropriately proportionate to the cost function. The exterior penalty is then calculated as^31–33

h_{m} (x) = \max {(0, (1 - KS (x)))}^{2}

(9)

where KS(x) is the Kreisselmeier–Steinhausser equation given as

K S (x) = g_{\max} + \frac{1}{ρ} \ln (\sum_{i = 1}^{N} e^{(ρ g_{i} (x) - g_{\max})})

(10)

where g_i is constraint violation of the i^th of N elements or nodes (depending on the constraint type) in a design space, with constraint violations of greatest magnitude being given by g_max. The constraint violation is calculated as

g_{m} (x) = \frac{| C_{calc} (x) |}{C_{allow} (x)}

(11)

where

C_{calc} (x)

and

C_{allow} (x)

indicate the calculated and allowable value under consideration, respectively.

Now, an optimization problem is produced that considers constraints without directly constraining the design space, but has infeasible points penalized by equation (10).

The penalized mass function

The penalized mass function has been produced in terms of base mass in addition to all constraints applied as penalties in the form

\begin{array}{l} Φ_{mass} (x) = M (x) + W_{stress} (\sum_{m = 1}^{32} h_{m} (K S (g_{m} (x)))) \\ + W_{disp} (\sum_{m = 33}^{36} h_{m} (K S (g_{m} (x)))) \\ + W_{modal} h_{37} (g_{37} (x)) + W_{geom} h_{38} (g_{38} (x)) \end{array}

(12)

where Φ(x)_mass is the penalized mass function, M(x) is the base mass, h_m is the m^th constraint penalty function calculated according to equation (9), which in turn is a function of the KS equations. For the KS equations, constraint violation factors of g_m are determined according to the equations in Table 3. The base mass function is calculated from the FEM, which models the sandwich panel as a three-layer composite with two outer faces and an inner core. The three-core sandwich design variables, namely, L_c, t_c, and H_c, are used to determine the sandwich core properties of density and stiffness employing the homogenization equations developed by Burton and Noor.³⁴ Such equations of homogenized properties define sandwich core stiffness in both the flexural and shear planes in a conservative manner. Thus, the homogenized sandwich core stiffness is considered conservative in both in bending and in shear, and any displacements from shear deformation are harmonized throughout the optimization process.

Table 3.

Constraints for the penalized mass function. Note that for upper bound limits, the constraint is C_calc/C_allow, while the inverse is true for lower bound constraints. All constraints use absolute values.

Constraint	Number	Equation	Description
Face yielding	m = 1, …, 4	${(g_{f Y})}_{m} = \frac{σ_{f}}{σ_{Y}}$	Face yielding constraint for all four load cases. Enforced using element von Mises stress in the faces at the corners.
Intracellular buckling	m = 5, …, 8	${(g_{i c b})}_{m} = \frac{σ_{3}}{σ_{i c b}}$	Intracellular buckling constraint for all four load cases. Enforced using minor principal stress in the faces at the corners and equation (2b).
Face wrinkling	m = 9, …, 12	${(g_{w r})}_{m} = \frac{σ_{3}}{σ_{w r}}$	Face wrinkling constraint for all four load cases. Enforced using minor principal stress in the faces at the corners and equation (2c).
Compressive crimping	m = 13, …, 16	${(g_{c r})}_{m} = \frac{σ_{3}}{σ_{c r}}$	Shear crimping constraint due to excessive compression. Enforced using minor principal stress in the faces at the corners and equation (2d)
Shear crimping	m = 17, …, 20	${(g_{scr})}_{m} = \frac{τ_{\max}}{τ_{s c r}}$	Shear crimping constraint due to excessive shear. Enforced using shear in the faces and cores at the corners and equation (2e)
Shell yielding	m = 21, …, 24	${(g_{shells})}_{m} = \frac{σ_{V M m a x}}{σ_{Y}}$	Constraint to ensure no failure occurs in nearby shell elements using the von Mises failure criterion. Enforced using stress in default Z-offsets in shell element corners.
Beam yielding	m = 25, …, 28	${(g_{beams})}_{m} = \frac{σ_{V M m a x}}{σ_{Y}}$	Constraint to ensure no failure due to excessive stress in beams using the von Mises failure criterion. Enforced using maximum beam stress in beam elements at all points.
Crippling	m = 29, …, 32	${(g_{c c})}_{m} = \frac{σ_{\min}}{σ_{c c}}$	Constraint to ensure crippling does not occur in beams. Enforced using minimum beam stress in beam elements at all points and equation (7).
Displacement	m = 33, …, 36	${(g_{disp})}_{m} = \frac{d_{\max}}{d_{lim}}$	Displacement constraint to preserve rigidity in all load cases. Enforced using maximum norm of nodal displacement vector and displacements in ESM.
Natural frequency	m = 37	$g_{modal} = \frac{f_{lim}}{f_{n a t}}$	Lower limit on natural frequencies in the structure per MSRP regulations.²³
Thin faces	m = 38	$g_{geom} = \frac{t_{f} + H_{c}}{5.77 t_{f}}$	Constraint to ensure thin face criteria given by Allen is kept.²⁰

The penalized compliance function

The penalized compliance function is written as base compliance given by FEM in addition to all constraints applied as penalties in the form

\begin{array}{l} Φ_{comp} (x) = \sum_{j = 1}^{4} C_{j} (x) + W_{geom} h_{38} (g_{38} (x)) \\ + W_{mfrac} h_{mfrac} (\frac{M (x)}{M_{baseline}}) \end{array}

(13)

where C_j is the compliance for the j^th load case as a function of the design variables, M(x) is the calculated mass, M_baseline is the baseline mass, and W is once again a weighting value assigned to appropriately penalize infeasible design points, where its subscript indicates the affiliated constraint. Because all constraints are global, the KS equations have not been employed, and penalties are applied using equation (9) exclusively as a function of g. Relevant constraints can be referenced from Table 3. Note that where possible, element data such as stress has been evaluated at the element corners. This is because the stress at the element center is an extrapolated stress from all elemental corners. Therefore, the element corner data is always more conservative than that of the centroidal data and is therefore used for all stress evaluations. Where compressive stresses are of concern, minimum principal stress, denoted σ₃, is evaluated.

The procedure shown in Figure 6 is followed. First, a sub-model is defined and validated to reduce cost function evaluation time. After an appropriate sub-model is selected, two cost functions are defined using equation (8): one for a penalized mass, and another for penalized compliance. To gain a better understanding of the design problem and appropriately apply penalty weightings, response surfaces are created using each penalized cost function with constraints enforced using the Kreisselmeier–Steinhausser equations. Next, an individual optimization is run for minimal penalized mass and compliance. Finally, a Pareto front is found for the design space to determine the ideal mass-to-compliance ratio for the panel. Note that it is not always necessary to execute this solution in full. It is certainly possible to enter straight into the Pareto front generation. However, it is important to note that the Pareto optimization program requires a many function evaluations to produce a meaningful front. The motivation for each intermediary stage is that it leads to an increased understanding of the response of the optimization problem and the locations of the local optima. By finding the local optima and using them as seed points for each subsequent step, optimization search times are drastically reduced, leading to a faster solution.

Figure 6.

Overview of the Panel-By-Panel MSDO procedure.

Case study and results

Stage I—Topology optimization of a global hopper car structure

Topology optimization is performed by minimizing compliance of the structure, and with a volume fraction limit of one quarter. Results for each iteration of the topology optimization are shown in Figure 7, with initial models for each iteration on the left and element density contour plots on the right. From the element density contour plot shown in Figure 7(f), it is apparent that most of the wall structure is of low relative density, indicating that the beams are providing most of the rigidity in the design shown in Figure 7(e).

Figure 7.

Element density contour plots from the different topology optimization iterations in ascending order. Shown on the left are the beam layouts used in each iteration. Black indicates beams which are added for the present iteration, while light gray indicates beams which were initially present or present in previous iterations.

Upon completion of the topology optimization procedure, size optimization is performed using Altair OptiStruct²⁹ for minimal mass subject to the constraints in Table 2. The mass of the half model converges to a final value of 6644.2 kg, implying a total car-body mass of 13,288.4 kg. Now that a car-body with optimal beam placement is identified, sandwich panels are readily applied in the car walls and roofs with less concern for major changes in rigidity and stresses.

Stage II—Multiscale design optimization of a hexagonal honeycomb panel wall on a hopper car

To determine the sandwich variables giving the structure of minimum weight, a MSDO is performed on the walls of the hopper car. A reduced order model using the Superelement method is produced for timely results on the design space of interest. In Figure 8, the wall panel of interest for this case study is highlighted in black. The residual sub-model is highlighted in dark gray, while the portion of the model reduced input via the Superelement method is highlighted in light gray.

Figure 8.

Depiction of the sub-model of interest. Here, the black portion is the design space, the dark gray is the residual sub-model after model reduction, and the light gray is the portion of the model reduced to direct matrix input.

An iterative approach is taken to determine the cost function weighting factors. The base cost function is plotted next to the penalized cost function for mass and compliance and the penalty factor is slowly increased. Eventually, the penalty is large enough such that the feasible region could be easily distinguishable from the infeasible region. An infeasible point is then selected, and the cost function is re-evaluated with varying weighting factors. The weighting is increased until it is determined qualitatively that the applied penalty for infeasibility meaningfully penalizes the cost, while still resulting in smooth behavior near the barriers of the feasible design space. Weightings determined from this method are shown in Table 4.

Table 4.

Weighting values for penalty factors by constraint category.

Constraint	Weighting
Stress	100
Displacement	100
Modal	200
Mass fraction	5000
Geometric	1000

Response surfaces are created by bringing into play the listed weighting factors and running the cost function in a sample set of 625 equidistant points for two of the four variables at a time, while fixing the other two. Linear interpolation between points is used to create the resultant response surface. Figure 9 shows response surfaces with respect to face thickness, t_f, and core thickness, H_c. Conversely, Figure 10 shows the response surfaces with respect to core wall width, L_c, and core foil thickness, t_c. Response surfaces achieve two primary goals. First, response surfaces allow for confirmation that weightings applied to the KS equations appropriately penalize the design. Second, response surfaces offer preliminary design intuitions and a starting point for optimization algorithms.

Figure 9.

Objective contour plots as a function of t_f and H_c to test sensitivity and determine close-to-optimum starting points. For this plot, L_c = 6.35 × 10⁻³ m and t_c = 5 × 10⁻⁴ m.

Figure 10.

Objective contour plots as a function of L_c and t_c to test sensitivity and determine close-to-optimum starting points. For this plot, t_f = 5 × 10⁻⁴ m and H_c = 6.35 × 10⁻³ m.

Response surfaces for mass and compliance plots with respect to face and core thickness fall in line with sandwich theory. As expected, a linear relationship exists between mass of the panel and both face and core thickness. The penalized mass plot, shown in Figure 9(c), shows that as face thickness is increased, the linear relationship is preserved. This indicates that for a thick face, the design remains feasible despite mass increasing. Meanwhile, by increasing core thickness and maintaining a thin face, the design enters the infeasible design space and is penalized by the KS equations due to a high susceptibility to wrinkling per equation (2b). The compliance contour plot appears to show compliance is minimized with a thick face, while core thickness has a lesser effect. This disproportionate effect of face thickness on compliance reduction occurs because the upper bounds of the face thickness are allowed as high as 6.35 mm. The penalized compliance contour plot in Figure 9(d) assigns an appropriate cost to penalize this phenomenon, as a thicker face is infeasible due to the resulting excessive mass, thus resulting in a corrected cost function for the MSDO.

In Figures 10(a) and (c), mass and penalized mass response surfaces are shown with respect to core wall width and core foil thickness. The optimal point is near the barrier of the feasible region because a large core wall width and small foil thickness can reduce mass significantly as they are the sole factors in core density but can result in high susceptibility to failure of all forms. Therefore, the optimal point should be close to the barrier where mass is the lowest without initiation of failure.

Finally, Figures 10(b) and (d) show compliance and penalized compliance with respect to core wall width and core foil thickness. As predicted by theory, compliance has virtually no sensitivity to these variables since core density and stiffness is not a factor of importance unless it approaches the density of the faces. At the top left of the plots, where wall thickness is high and wall width is low, the core effectively behaves more as a solid material, and so an improvement in compliance is observed. Based on the penalties encroached in these areas from the KS equations due to unreasonably high mass, it is demonstrated that modification of these variables for the purpose of improving compliance is inefficient.

Points of lowest cost from all penalized cost response surface plots are denoted by the asterisks. Since two of the four design variables are fixed for the creation of the response surface, the allotted material for the design of the sandwich panels is not be optimally distributed between the sheet thicknesses and the core geometry. By plotting the response surfaces, only near optimal points are found which accordingly give useful design intuitions for refined optimization.

A seed point for the penalized mass from a combination of the best t_f and H_c values found from the response surface in Figure 9(a), and the L_c and t_c values from Figure 10(a) is produced. A similar seed point is generated for penalized compliance using the points of lowest cost from Figure 9(b) and Figure 10(b). A MATLAB algorithm employing these seed points with the patternsearch subroutine is used to obtain the new local optima for each cost function, with results shown in Table 5. Now, the newly found optima along with those from the response surfaces can be used as seed points in the next step, where a multi-objective optimization is performed.

Table 5.

Seed points for the Pareto search algorithms.

Source	L_c	t_c	t_f	H_c
Source	$[m \times 10^{3}]$
Response surface	7.144	0.287	0.569	3.175
Individual optimizations	8.275	0.378	0.557	3.176
Individual optimizations	9.523	0.024	1.386	12.160

The Pareto front produced by a multi-objective optimization for minimal mass and compliance is shown in Figure 11. Note that the mass fraction penalty constraint in the penalized compliance function has been removed for the multi-objective optimization, since mass is now part of the objectives. Two algorithms are employed to produce one front each: the first uses a pattern search-based approach, while the second uses a genetic algorithm.

Figure 11.

The Pareto front resulting from the multi-objective optimization of penalized mass and total unpenalized compliance. The Utopia point has a mass of M(x) = 826 kg and a total compliance of C(x) = 18, 726 m/N.

Selection of Pareto-optimal solutions from each front a Utopia Point is carried out using definitions provided by Marler and Arora according to the theorem³⁵

F^{°} \Leftarrow \Rightarrow F_{i}^{°} = \min {F_{i} (x) | x \in X} \forall i \in {1,2, \dots, k}

(14)

where F^° denotes the Utopia point, F_i° denotes the i^th of k objectives in the multi-objective optimization with corresponding design variable x in the Pareto feasible design space,

X

Baseline structures are produced employing a sizing optimization for minimal mass on the hopper car structure at the end of Stage I. Solid steel is used instead of honeycomb sandwich panels, with constraints applied per Table 2. Further details on the generation of the baseline design can be found in the previous work of Al-Sukhon.²⁴ In this case, a feasible design space has been defined such that a minimum reduction in mass of 8% as compared to the baseline panel mass of 983.34 kg, so at least a 4% reduction in emissions is produced.^1,2 The feasible design space is therefore limited to the left of the corresponding dashed line in Figure 11. For reference, the baseline mass is shown as well. Design points on each Pareto front of nearest proximity to the Utopia point are selected as Pareto-optimal designs and listed in Table 6.

Table 6.

Pareto-optimal designs from each Pareto front in Figure 11 for each search algorithm.

Search algorithm	L_c	t_c	t_f	H_c	m(x)	C(x)	Function evaluations
Search algorithm	$[m×1 0^{3}]$				[kg]	[m/N]	Function evaluations
Pattern search	9.5250	0.4874	1.0716	5.9531	903.58	18,726	2358
Genetic algorithm	7.3909	0.2999	0.9956	6.7820	893.24	18,825	16,523

Discussion

The Pareto-optimal points have been determined for the design of a characteristic panel. For a full design, it is possible to repeat the optimization in more panels for a further reduction in mass. For the sake of brevity, only a single panel is considered in this case study. Though the genetic algorithm search provides a lower result for mass, it required eight times as many function evaluations and produced half the number of pareto front points. Therefore, results from the pattern search are used for the final result. To demonstrate the viability of the results, the Pareto-optimal design from Table 6 employing the pattern search is extrapolated over the full structure and re-analyzed.

As a sample result, Figure 12 shows the von Mises stresses in the sandwich panel faces for all load cases considered. A full set of structural results are listed in Table 7. Note that while the structure experiences yielding in one local area upon results extrapolation, further sub-division of the optimization into smaller panel sizes is certain to resolve the failure since the maximum stress barely exceeds the allowable limit. The lowest natural frequency from the modal analysis is 12.08 Hz, which is above the allowable lower limit of 12.0 Hz.

Figure 12.

Element stress contour plots for maximum von Mises face stresses in element centers. Stress averaging is applied for post-processing, therefore element centers are used.

Table 7.

Stress results from extrapolation of optimized side panel to all wall and roof components of the structure.

Component	Stress type	Observed stress [MPa]				Allowable stress [MPa]
Component	Stress type	LVL	DTL	ECL	FSI	Allowable stress [MPa]
Faces	Max VM	169.9	373.7	367.0	396.6	423.3
Faces	Min P3	−173.7	−265.6	−322.8	−329.9	−516.6
Core	Shear XZ, YZ	37.9	131.6	29.7	126.2	174.56
Shells	Max VM	414.5	410.7	415.1	457.6	423.3
Beams	Max VM	173.4	491.6	429.5	416.7	423.3
Beams	Min norm	−89.25	−131.1	−288.8	−216.7	Variable

The baseline structure has a half model mass of 8518.8 kg, while the structure at the end of Stage I has a mass of 7110.0 kg. The structure produced by the presented method employing hexagonal honeycomb sandwich panels has produced a mass of 5731.4 kg without compromising on overall rigidity of the structure. This implies that an average mass savings of 32.72% can be anticipated by employing the novel multiscale design optimization procedure, while an average mass reduction of 13.74% can be expected over the use of conventional materials without compromising the rigidity of the structure.

Conclusion

In this paper, a novel multi-stage and multiscale design optimization methodology is introduced and successfully employed in the weight minimization of a covered hopper car. An analysis of the literature reveals that hexagonal honeycomb sandwich panels are demonstrated to offer superior strength-to-weight ratios under the expected loading experienced by a hopper car, thus making them the ideal structural candidate for optimization. A two-stage optimization procedure is then established, with the first stage determining the ideal frame structure by means of topology optimization. Meanwhile, in the second stage, a novel MSDO is introduced whereby the design variables of a characteristic honeycomb sandwich panel are optimized at the microscale and homogenized to the mesoscale to predict structural integrity. A Pareto front is produced by the MSDO, for which the Utopia point is found to determine the Pareto-optimal point. Assuming the car-body is half the mass of the hopper car, a mass savings of 16.36% can be expected, suggesting a reduction in greenhouse gas emissions of 13.09% based on available literature. Though implementation of this methodology is costly as compared to traditional methodologies, results show that for the case of a hopper car, it is able to produce results which are more reliable than those traditional methods, which can produce optimal design variables that result in failure.²⁴

Footnotes

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors acknowledge the financial support provided by Transport Canada Clean Rail Academic Grant (Grant Number 300515).

ORCID iD

Ayman Al-Sukhon

Appendix

The derivation of the wall load index is presented here as an example to show how the load indices plotted in Figure 1 are determined. Per the AAR MSRP, the pressure at the base of a sidewall on a covered hopper car is given as²³ (15)

P = W h_{w} (\frac{1 - \sin (ϕ)}{1 + \sin (ϕ)})

where W is the specific weight of the cargo assumed as 850 kg/m³ throughout the paper, h_w is the wall height assumed as 2.5 m, and ϕ is the static angle of repose assumed as 25°.

The AAR MSRP mandates that the load in equation (15) be applied as a linearly distributed load decreasing to zero at the top. For a simply supported beam under a linearly distributed load, the maximum bending moment using basic statics is (16)

M_{\max} = \frac{w_{o} L^{2}}{9 \sqrt{3}}

while the maximum shear load is (17)

V_{\max} = \frac{w_{o} L}{3}

where w_o is the load per unit width at the base of loading and L is the panel length. Now, substituting in the equations for M_max and V_max into equation (3), the load index for a simply supported beam under a linearly distributed load is given as (18)

Π = \frac{V_{\max^{2}}}{E M_{\max}} = \frac{w_{o}}{\sqrt{3} E}

Considering that the pressure at the base of the wall is effectively the load per unit length per unit with of the wall, equation (15) into equation (18) to obtain the wall load index as (19)

{(Π)}_{wall} = \frac{W h_{w}}{\sqrt{3} E} (\frac{1 - \sin (ϕ)}{1 + \sin (ϕ)})

References

Helms

Lambrecht

Höpfner

. Energy Savings by Light-Weighting. Technical report, IFEU-Institute, 2003.

Helms

Kräck

. Energy Savings by Light-Weighting - 2016 Update. Technical report, IFEU-Institute, 2016.

Rochard

Schmid

. Benefits of lower-mass trains for high speed rail operations. Proc Inst Civ Eng - Transp 2004; 157(1): 51–64. DOI:10.1680/tran.2004.157.1.51.

Lee

Kim

Sun

, et al. The next generation material for lightweight railway car body structures: Magnesium alloys. Proc Inst Mech Eng F J Rail Rapid Transit 2016; 232(1): 25–42. DOI:10.1177/0954409716646140.

Seo

Kim

Cho

. Development of a hybrid composite bodyshell for tilting trains. Proc Inst Mech Eng F J Rail Rapid Transit 2008; 222(1): 1–13. DOI:10.1243/09544097jrrt96.

Cho

Koo

Jung

. A lightweight design approach for an EMU carbody using a material selection method and size optimization. J Mech Sci Technol 2016; 30(2): 673–681. DOI:10.1007/s12206-016-0123-8.

Zinno

Fusco

Prota

, et al. Multiscale approach for the design of composite sandwich structures for train application. Compos Struct 2010; 92(9): 2208–2219. DOI:10.1016/j.compstruct.2009.08.044.

Catapano

Montemurro

. A multi-scale approach for the optimum design of sandwich plates with honeycomb core. Part II: the optimisation strategy. Compos Struct 2014; 118: 677–690. DOI:10.1016/j.compstruct.2014.07.058.

Wennberg

Stichel

Wennhage

. Substitution of corrugated sheets in a railway vehicle’s body structure by a multiple-requirement based selection process. Proc Inst Mech Eng F J Rail Rapid Transit 2012; 228(2): 143–157. DOI:10.1177/0954409712467139.

10.

Wennberg

Stichel

. Multi-functional design of a composite high-speed train body structure. Struct Multidiscipl Optim 2014; 50(3): 475–488. DOI: 10.1007/s00158-014-1056-4.

11.

Molavitabrizi

Laliberte

. Methodology for multiscale design and optimization of lattice core sandwich structures for lightweight hopper railcars. Proc Inst Mech Eng C: J Mech Eng Sci 2020; 234(21): 4224–4238. DOI: 10.1177/0954406220920694.

12.

Lee

Jung

Jang

, et al. Structural-optimization-based design process for the body of a railway vehicle made from extruded aluminum panels. Proc Inst Mech Eng F J Rail Rapid Transit 2015; 230(4): 1283–1296. DOI:10.1177/0954409715593971.

13.

Harte

McNamara

Roddy

. A multilevel approach to the optimisation of a composite light rail vehicle bodyshell. Compos Struct 2004; 63(3–4): 447–453. DOI: 10.1016/s0263-8223(03)00193-4.

14.

Kruggel-Emden

Simsek

Rickelt

, et al. Review and extension of normal force models for the discrete element method. Powder Technol 2007; 171(3): 157–173. DOI: 10.1016/j.powtec.2006.10.004. URL https://www.sciencedirect.com/science/article/pii/S0032591006004360.

15.

Di Maio

Di Renzo

. Analytical solution for the problem of frictional-elastic collisions of spherical particles using the linear model. Chem Eng Sci 2004; 59(16): 3461–3475. DOI: 10.1016/j.ces.2004.05.014. URL https://www.sciencedirect.com/science/article/pii/S0009250904003148.

16.

Govender

Wilke

Kok

. Blaze-demgpu: modular high performance dem framework for the gpu architecture. SoftwareX 2016; 5: 62–66. DOI: 10.1016/j.softx.2016.04.004. URL https://www.sciencedirect.com/science/article/pii/S235271101630005X.

17.

Fang

, et al. Quasi-real-time simulation of rotating drum using discrete element method with parallel GPU computing. Particuology 2011; 9(4): 446–450. DOI: 10.1016/j.partic.2011.01.003. URL https://www.sciencedirect.com/science/article/pii/S1674200111000654.

18.

Ushakov

Safonov

Sergeichev

, et al. Design and optimization of a vacuum infusion technological process for hopper car fabrication using polymeric composite materials. J Mach Manuf Reliab 2015; 44(3): 276–282. DOI:10.3103/S105261881503022X.

19.

Introduction to Explicit Analysis with Altair Radioss. Technical report, Altair Inc., 2019.

20.

Allen

. Analysis and Design of Structural Sandwich Panels. Elsevier, 1969. DOI:10.1016/c2013-0-02134-2.

21.

Petras

Sutcliffe

. Failure mode maps for honeycomb sandwich panels. Compos Struct 1999; 44(4): 237–252. DOI:10.1016/s0263-8223(98)00123-8.

22.

Rathbun

Zok

Evans

. Strength optimization of metallic sandwich panels subject to bending. Int J Sol Struct 2005; 42(26): 6643–6661. DOI:10.1016/j.ijsolstr.2005.06.044.

23.

Association of American Railroads , Publications Department Transportation Technology Center, Inc. P.O. Box 11130 55500 DOT Road Pueblo, CO 81001. Section C Part II Design, Fabrication, and Construction of Freight Cars, M-1001, 2015 ed.

24.

Al-Sukhon

. Multiscale Design Optimization of Hopper Cars Employing Functionally Graded Honeycomb Sandwich Composites. Master’s Thesis DOI:10.22215/etd/2020-14178.

25.

Altair Engineering Inc. Radioss . URL https://www.altair.com/.

26.

Drucker

Prager

. Soil mechanics and plastic aanalysis or limit design. Q Appl Math 1952; 10(2): 157–165. URL http://www.jstor.org/stable/43633942.

27.

Attar

. Hopper Railcars: A Preliminary Study on the Use of Fibre Reinforced Polymers in Structural Components. Master’s Thesis, Carleton University, 2015. DOI:10.22215/etd/2016-11339.

28.

Bendsøe

. Optimal shape design as a material distribution problem. Struct Optim 1989; 1(4): 193–202. DOI:10.1007/bf01650949.

29.

Altair Engineering Inc. Altair optistruct . URL https://www.altair.com/.

30.

Niu

MCY

. Airframe Stress Analysis and Sizing. Adaso/Adastra Engineering Center, 2011.

31.

Kreisselmeier

Steinhauser

. Systematic control design by optimizing a vector performance index. IFAC Proc Volumes 1979; 12(7): 113–117. DOI:10.1016/s1474-6670(17)65584-8.

32.

Wrenn

. An Indirect Method for Numerical Optimization Using the Kreisselmeier-Steinhauser function. NASA Contractor Report, National Aeronautics and Space Administration, Langley Research Center, Office of Management, Scientific and Technical Information Division, 1989. URL https://books.google.ca/books?id=v6xT-Mcr8ocC.

33.

Love

Kapania

. Structural optimization of a novel flying wing supersonic aircraft configuration. In AIAA Scitech 2020 Forum. American Institute of Aeronautics and Astronautics. DOI:10.2514/6.2020-0166.

34.

Burton

Noor

. Assessment of continuum models for sandwich panel honeycomb cores. Comput Methods Appl Mech Eng 1997; 145(3–4): 341–360. DOI:10.1016/s0045-7825(96)01196-6.

35.

Marler

Arora

. Survey of multi-objective optimization methods for engineering. Struct Multidiscipl Optim 2004; 26(6): 369–395. DOI:10.1007/s00158-003-0368-6.

Design optimization of hopper cars employing functionally graded honeycomb sandwich panels