Reliable shape optimization of structures subjected to transient dynamic loading using genetic algorithms

This paper presents a reliable method of solution of two dimensional shape optimization problems subjected to transient dynamic loads using Genetic Algorithms. Boundary curves undergoing shape changes have been represented by B-splines. Automatic mesh generation and adaptive finite element analysis modules are integrated with Genetic algorithm code to carry out the shape optimization. Both space and time discretization errors are evaluated and appropriate finite element mesh and time step values as obtained iteratively are adopted for accurate dynamic response. Two demonstration problems have been solved, which show convergence to the optimal solution with number of generations. The boundary curve undergoing shape optimization shows smooth shape changes. The combinations of automatic mesh generator with proper boundary definition capabilities, analysis tool with error estimation and Genetic algorithm as optimization engine have been observed to behave as a satisfactory shape optimization environment to deal with real engineering problems.


Shape optimization
Structural shape optimization has long been an area of high research. Generally, shape optimization addresses the problem of determining the best shape for a given set of load conditions. In conventional shape optimization, an initial structural layout is assumed and the boundary curves/surfaces are altered. In order to provide a reliable and effective shape optimization tool, three main aspects are to be addressed. First, a geometry representation system capable to deal with the changeable geometry of the model during optimization process. Secondly, a structural analysis technique adequate to the boundary changes, which allows an easy re-meshing and its accurate boundary displacement and stress solutions. Finally, a robust optimization tool that finds out the optimum or quasi-optimum solution.

Geometric modeling in CAD/CAM environment
The geometrical model is where the design variables are easily imposed and allows explicit integration with other design tools such as CAD or CAM systems. Initially, many authors such as Zienkiewicz and Campbell [1], Ramakrishnan and Francavilla [2] among others did not use geometric modeling. Instead, they defined nodal coordinates of the discrete finite element model as design variables. This approach required a large number of design variables and produced jagged edges shapes. In order to overcome this problem a large number of constraints were added, which complicated the optimization process. Further, the lack of an associated geometric model did not allow the integration with powerful design tools like CAD or CAM systems. Since then, several methods were used to overcome initial drawbacks. Several researchers [3][4][5][6] used mesh parameterization methods to define geometric and analytical models without the complications described above. In this approach a set of key points or master nodes were used to define the geometric entities of the mesh e.g., lines, surfaces, volumes. Parametric mappings were employed to map these geometric entities to the finite element nodes. However, parameterization methods are difficult to use because the designer has to define the model in terms of master nodes rather than dimensions, which is not an easy task when dealing with complex models and they also suffer from possible mesh degradation for large design changes [7].
Belengundu and Rajan [8] proposed a different approach, which was the natural design variable method to solve shape optimization problems. Kodiyalan et al. [9], Chen and Tortorelli [10] used a method, where a solid model of the structure to be optimized was used rather than a finite element model.
The approach as adopted in our present study is based on parameterization techniques. In order to overcome its drawbacks, an automatic mesh generator has been used. The mesh degradation is avoided by defining constraints over the topology of the mesh and by the use of side constraints. On the other hand, the geometric modeling in CAD system can be carried out using Lagrange polynomials, Bezier curves, Cubic splines and B-splines curves. Among these approaches, the B-spline curve properties can be easily controlled by using its parameters. Other useful characteristics of using B-splines are that the curve lies in a convex hull of the vertices, individual sample points affect only the local part of the curve and they allow imposition of C 1 continuity. Further, a polynomial of certain degree requires a fixed number of sample points. However, with B-spline functions additional sample points can be introduced without increasing the order of the smoothening curve.

Heuristic optimization
The heuristic combinatorial search techniques, especially the Genetic algorithms (GAs), have been the matter of recent interest in the field of structural optimization. GAs are robust search and optimization algorithms, with good balance between efficiency and efficacy, necessary to survive in many different environments. Among other considerations, GAs do not need further additional information than objective function values or fitness function. In the field of structural optimization, discrete optimization problems arise quite naturally. Since, most of the traditional optimization methods treat the design variables as continuous, they are very much inadequate in the presence of discrete design variables. On the contrary, GAs start their search with a coding of parameters set which makes the search space discrete and therefore naturally suitable for discrete optimization problems [11,12].
Annicchiarico and Cerrolaza [13] solved 2-D shape optimization problems using GAs. Boundary β-spline curves with an automatic mesh generation were used to obtain the analytical model and finally, the finite element method was used to estimate the structural responses. They showed that the combination of GAs along with geometric modeling having β-spline boundary generation capability was highly advantageous for solution of shape optimization problems. Annicchiarico and Cerrolaza [14] also presented an integrated approach for 3-D structural shape optimization using β-spline surface representation, genetic algorithms and mesh parameterization. Holst and Pulliam [15] presented a method for aerodynamic shape optimization using GA with real number encoding, where it was observed that GA with real number encoding to be more efficient than with binary coding in terms of CPU time. Obayashi et al. [16] used multi-objective genetic algorithm for multidisciplinary design of transonic wing platform.

Optimization applied to dynamics problems
The work on structural optimization under transient dynamic loading is very scanty. Feng et al. [17] considered sizing optimization problem of elastic structures with fixed geometry under dynamic loads. They used FEM, modal analysis and a generalized steepest descent method in developing the computational algorithm. Structural weight is minimized subject to constraints on displacement, stress, structural frequency and member size. Using conventional mathematical programming, shape optimization under transient dynamic loading is rather complex due to the involved design sensitivity calculations. The sensitivity analysis with stress/displacement constraints can become computationally very expensive and erroneous if proper care is not taken in choosing the proper finite element mesh and time step size in the analysis [18,19]. Brach [20] considered the problem of minimizing the maximum dynamic deflection at the center of simply supported beams of constant total mass. Yau [21] also dealt the problem of optimal design of a simply supported beam of given total mass, minimizing the upper bound on its dynamic response. Icerman [22] and Mroz [23] presented optimum designs of structures that are excited by a single harmonic load. However, both of these treatments are limited to harmonic excitation and cannot be used for general transient problems. Masad [24] solved shape optimization problem of a nonuniform beam such that its fundamental natural frequency is a maximum.
Thus, after reviewing the existing literature, it is observed that very few works exist on structural shape optimization under transient dynamic load and utilizing heuristic combinatorial search techniques, like genetic algorithms. Further, the use of error estimation and adaptivity for controlled and reliable analysis have been used for structural shape optimization under only static loads by a few researchers [25][26][27]. At the present state of knowledge, the work recorded on adaptive finite element analysis under dynamic loads is also limited. Wiberg and Li [28] proposed a post-processed type of a posteriori estimates in space and also in time while using direct integration for dynamic response evaluation. It updated the spatial mesh and time step so that the discretization errors were controlled. This process is not only time consuming but also gives rise to new errors as nodal values need to be interpolated from the previous mesh to the newly generated mesh whenever a mesh is changed. Joo and Wilson [29] arrived at the final mesh using Ritz vector as basis of transformation. In their investigation, the authors made use of modal participation and amplification factor and obtained error estimates based on Babuska's criterion using amplified Ritz modes. Cook and Avrashi [30] discussed the procedure for estimating the discretization error of the finite element modeling as applied to the calculation of natural frequency of vibrations. Meshes were obtained corresponding to each mode. Dutta [31] proposed a measure for discretization error both in space and time, which is a logical extension of Zienkiewicz and Zhu [32] error criterion for space discretization error and Zieniewicz and Xie [33] error criterion for time discretization by involving time integration to consider the variation of response with time.
In the present paper, shape optimization of structures under transient dynamic loads using genetic algorithms have been carried out. Boundary B-spline curves with an automatic mesh generator has been used to obtain the analytical model. Finite element method has been used to obtain the transient dynamic response of the structure. Moreover, the reliability of the analysis is ensured by incorporating error estimation in both space and time following the strategy proposed by Dutta [31]. Appropriate finite element meshes and time steps have been utilized for user specified accuracy limit. Two demonstration problems, which have been solved show convergence to the optimal solution with different genetic parameters. The boundary curve undergoing shape optimization shows smooth shape changes as well. Such a controlled shape optimization of structure under transient dynamic loading based on GA, adaptive finite element analysis and B-splines can be observed to be highly reliable and has not been reported anywhere in the literature.

Implementing GA for optimization process
In order to use GAs in optimization problems, some parameters of interest in the system to be optimized have to be chosen. These parameters are called design variables. In the present study, the y-coordinates of position vectors of the points on the B-spline curve control vertices (y p ) are taken as design variables. These design variables are represented by some set of strings coded in binary or other codes, which corresponds to the chromosomes of living things. In the present case, an individual design is represented by a binary string of appropriate length incorporating, generally by simple concatenation, the values of all design variables. Design = y p1 , y p2 , . . . , y pn (1) Chromosome = 10011100 . . . 001110 These strings form the initial population. The variable y p is bounded between upper (y u p ) and lower limits (y l p ). The decimal value of the design variable can be computed from where q is the string length of binary coded design variable and b is the binary digit (1 or 0). In the present study, fifteen bits have been taken to code each of the design variables.
Once the population has been defined, a fitness function or objective function that measures the behavior of each individual in its environment has to be defined. This function provides a direct indication of the performance of each individual to solve the optimization problem subjected to the imposed constraints from the environment. With the population ranked according to the fitness, a group of chromosomes are selected from the population. There exist several methods to select parents.
Following schemes can be used:

1) Stochastic sampling with replacement (SSWR). 2) Remainder stochastic sampling without replacement (RSSWR) 3) Tournament selection (TS)
Selection methods can be classified roughly into two groups: fitness-proportionate and rank-based selection. Fitness-proportionate methods select individuals probabilistically depending on the ratio of their fitness and the average fitness of the population. Both Stochastic sampling with replacement (SSWR), which is a fancy name of roulette wheel selection and Remainder stochastic sampling without replacement (RSSWR) are included in this category. The scheme like RSSWR has been proposed to reduce the stochastic error associated with spinning the roulette wheel numerous times. Tournament selection comes under the category of rank-based selection. The method works by first picking s individuals (with or without replacement) from the population and then selecting the best of the chosen s individuals. If performed without replacement in a systematic way, this selection scheme can assign exactly s copies of the best individual to the mating pool at every generation. The control of this selection pressure has made tournament selection popular. However, with high selection pressure, a rapid convergence of the population towards solution may be observed, which may not readily evolve to acceptable solutions. In the present study, a binary tournament selection with s = 2 is used, which is generally used in most GA application.
Using the crossover and mutation operators, the selected chromosomes are then reproduced. The crossover operation consists in taking two selected chromosomes as parents. Then, they are either crossed by using a certain probability value in order to create two new chromosomes (children) or they are directly included into new population. For the binary string in the present study, single point crossover has been used. Mutation is a randomly applied change, which is incorporated to a single gene to simulate copying errors in real organisms. This change is applied with a probability defined by the mutation rate. This operation is performed with the help of a random number in the range of 0 to 1. If the random number is less than the probability of mutation, then the bit under consideration will be switched (i.e. 0 to 1 or 1 to 0).
The values of probability of crossover and mutation have been chosen from those suggested by Goldberg [34]. High probability of cross over and low probability of mutation (inversely proportional to population) are considered. Crossover operator ensures the exchange of characteristics in the mating pool among strings selected based on their fitness and unless the probability is high, the chance of good intermixing of genes may be low. On the other hand, a mutation is the occasional random alteration of a binary digit. When used sparingly (a low probability is needed) with the reproduction and crossover operators, premature loss of important genetic material at a particular position can be avoided. Termination criteria has been set on the basis of number of generation and convergence to the optimum result has been ascertained on the basis of the nature of distribution of best fit solution, constraints and all the associated variables.
The optimization scheme can be summarized as follows: 1) Select Parents 2) Apply genetic operators in order to create the next generation.
3) Evaluate objective function and constraints. 4) If the optimum has not been reached yet, repeat 1-4 until the optimum is found or the process will be stopped by the end criteria.
The bulk of Genetic Algorithms processing power is due to the simple transformation carried out by selection and crossover operators. Mutation plays a secondary role in the operation of genetic algorithm and it is needed because even though selection and crossover effectively search and recombine extant notions, occasionally they may become overzealous and they can lose some potential useful genetic material (1's or 0's at a particular locations). Thus, the mutation operator protects against such irrecoverable premature loss of important notions.

Adaptive finite element analysis
A controlled finite element analysis of structure subjected to transient dynamic load has been carried out for reliable evaluation of response. This is ensured by carrying out error estimation for both space and time discretization. Appropriate mesh size and time step, ∆t values are selected iteratively for a prescribed accuracy limit.

Space discretization error
Finite element stresses are calculated as where D and B are usual elasticity and strain displacement matrices respectively. The approximate solution σ containing discretization errors differs from the accurate solution σ * and the difference is the error e σ Thus, e σ = σ * − σ A recent and elegant technique for the evaluation of accurate stress is done using superconvergent patch recovery method by Zienkiewicz and Zhu [35]. In the recovery process, it is assumed that the accurate nodal values σ * belong to a polynomial expansion σ * e of the same complete order p as that present in the basis function N and which is valid over an element patch surrounding the particular assembly node considered. Such a 'patch' represents a union of elements containing this vertex node. This polynomial expansion will be used for each component of σ * e and we get σ * e = Pa (5) where P contains the appropriate polynomial terms and a is a set of unknown parameters. For a general plate and shell problem, the polynomial can be expressed as a function of x, y and z. However, since only plate problems are solved in this paper, x and y are consider only as the parameters in the polynomial expression. Thus it can be written as and The determination of the unknown parameters a of the expression given in Eq. (5) is best made by ensuring a least square fit of this to the set of superconvergent points existing in the patch considered. To do this we minimize where (x i , y i ) are the co-ordinates of a group of sampling points, n = mk is the total number sampling points and k is the number of the sampling points on each element m j (m j = 1, 2, . . . , m) of the element patch. The minimization condition of F (a) implies that a satisfies This can be solved in matrix form as Once the parameter a are determined the recovered nodal vales of σ * are simply calculated by insertion of appropriate co-ordinate into the expression for σ * e . The stresses in the nodes inside the patch are evaluated using a from Eq. (10). The stresses at boundary nodes can be determined using interior patches.
The mathematical basis and broad details of the algorithm are described in [31]. The complete procedure of discretization error measure is given below in an algorithmic form.
Calculate smooth stress σ * (j) Energy norm for error e σ at element level Energy norm for FEM solution at element level } for the whole structure (ii) Overall Domain Discretization Error where e m = η u and − η is the permissible domain discretization error Time averaged error at element level, If η η and Where φ = 1/p for no singularity and φ = 1/λ for element close to singularity (p is the order of defining polynomial and λ is the strength of singularity) Endif } Go to step 1

Time discretization error
In order to control the time discretization error in the time marching scheme, some means of estimating the errors should be introduced so that the steps can be suitably adjusted. The adaptive time stepping is aimed at maintaining largest possible step size while keeping the accuracy within the prescribed limit. In this work, the error measure proposed by Zienkiewicz and Xie [33] is adopted with some modification. The total time domain T is divided into a finite number (n T ) of sub-domains. A uniform value of ∆t is maintained in a sub-domain and the time discretization errors are calculated. Time integration of errors is carried out over individual sub-domain interval to arrive at a new value of step size for that sub-domain.
The local error in time is estimated as wherez n+1 is acceleration vector at (n+1) th time step,z n is acceleration vector at n th time step and β = 1/4. The L 2 norm of the local error is Usually, it is inconvenient to specify an absolute error tolerance. To measure the relative error, it may be defined as where ( u ) max is the maximum value of the corresponding norm of the displacement solution recorded during the computation. Time integration of error over sub-domain leads to where T i is the time interval of sub-domain 'i'. In this section, the error η Ri in a sub-domain 'i' is used to calculate new step size for the next iteration so that the error is roughly equal to the prescribed tolerance. With this tolerance given as η t , an upper limit (γ 1 η t ) and lower limit γ 2 η t is also specified, where the parameter γ 1 and γ 2 are in the ranges of 0 γ 1 1 and γ 2 1.
When the error η Ri exceeds the upper limit the step size needs to be reduced and the new step size in a sub-domain 'i' may be predicted as Similarly, if the error η Ri is smaller than lower limit, then the step size may be increased using the above expression. Hence the overall strategy can be written as below.
Starting with a coarse mesh and a value of ∆t for every time sub-domain, overall domain discretization error and discretization errors at element level are calculated and time discretization errors are calculated for every time sub-domain. Mesh is then adaptively refined based on the space discretization error values at element level and new values of ∆t are obtained for every time sub-domain. Iteration is carried out until the mesh is an optimal one and time discretization errors in all the sub-domains are within the prescribed limit.

The overall optimization process
An integrated approach combining geometric modeler with mesh generation capability, reliable dynamic FE analysis module and an optimization engine (i.e. GA) is considered to simulate optimal design of structures subjected to transient loading satisfying all imposed constraints e.g., stress constraint, geometric constraints etc. GAs are usually not efficient (in terms of CPU time) as compared to deterministic methods, since they evaluate the objective function many times. However, the present study is intended for obtaining a robust and reliable optimization technique for complex problems. Thus GA is utilized as optimization tool and optimal solution is obtained with number of generations. The flow chart for the simulation process is shown in Fig. 1.
The stages involved can be summarized as below: First, the geometric model is created. Key points for generating the boundary B-spline curve, which has to undergo shape change, will be defined. Number of divisions in each direction (x and y) is specified for finite element modeling using 8-noded quadrilateral elements. Thus, the structural domain under consideration is discretized using valid FE model along with a suitable definition for the changing boundary. Appropriate boundary conditions, loading and material properties (modulus of elasticity, Poisson's ratio, density etc.) are specified. Two-dimensional dynamic analysis is carried out using Newmark-β integration method considering β and γ as 0.25 and 0.5 respectively and for an arbitrarily chosen time interval (∆t). Nodal stresses are evaluated from gauss point stresses using superconvergent patch recovery [35] corresponding to each time step. Error estimation is carried out for both space and time discretization. Appropriate mesh and ∆t values are selected iteratively for a prescribed accuracy limit.
Genetic Algorithm is used to optimize the objective function subjected to imposed constraints. The objective function is stated by Additional side constraints are also imposed, which depend on the problem to be optimized. They are used to prevent possible occurrence of any unacceptable changes in the configuration of the geometric model.
In order to incorporate all the constraints, penalty method has been used. The evaluation function, which is the minimization of the model area, penalized with the constraints is written as Where denotes the absolute value of the operand, if the operand is negative and returns a value zero, otherwise. ∆σ j is the allowed stress minus acting stress, ∆y p k is the side constraints violation in terms of control point coordinates. The parameters λ, µ are adjusted by trial and in the present study, they have been evaluated in such a way that 10% of violation in every constraint increases the original area by about 10%. If the convergence criteria are satisfied and the optimum solution has been found the simulation process is terminated otherwise a new set of values of design variables will result in the optimization iteration. New coordinates of shape controlling points are sent to the mesh generator, which automatically generates a new FE analysis model and the whole process is repeated.

Numerical examples
In order to demonstrate the robustness and power of the integrated shape optimization module, two examples are presented and discussed. The objective in both the cases is to find optimal surface area of the structures subjected to constraints on Von-Mises stress in addition to some geometric constraints. All the design variables (control vertices points) are considered within a lower and upper bound and the binary coded variables have been generated randomly. Different combinations of genetic parameters have been tried for both the examples. Each problem corresponding to each combination of genetic parameters is simulated three times from different initial population and then each result is shown as the average of the three results. The maximum, average and worst values of the objective functions are also noted.

Cantilever beam subjected to suddenly applied loading at the free end
The initial geometry, dimensions, boundary conditions, loading and prescribed design variables of a cantilever beam is shown in Fig. 2. The modulus of elasticity, Poisson's ratio and mass density of the material is taken as 2.1e11 N/m 2 , 0.3 and 7850 kg/m 3 respectively. The objective of the problem is to find out the optimal surface area provided that Von-Mises stresses anywhere in the structure must not exceed 250 MPa, which is computed over the  Step loading applied at the free end of the cantilever.  whole time domain. The boundary AB, which has to be controlled for the optimal solution is defined by a B-spline with six shape controlling points and Y coordinates of these points have been considered as design variables y 1 , y 2 , . . . , y 6 (Fig. 2). A suddenly applied step loading of constant magnitude of 75 kN acting at the free end of the cantilever beam for a duration of 0.05 sec is considered (Fig. 3). A finite element model of the cantilever beam using 8 noded isoparametric elements has been considered. The initial model is having 36 elements with an area of 1.2 m 2 (Fig. 4). A time step of magnitude 1.667e-3 sec (T 1 /15) has been chosen for time integration using Newmark-β integration method (constant average acceleration). The dynamic analysis is carried out and Von-Mises stresses on all the nodes are evaluated corresponding to each time step. Both space and time discretization errors are evaluated and have been obtained as 15.52% and 12.91% respectively. For a target accuracy of 3% and 1% for space and time discretization respectively, a finite element mesh having 108 elements as shown in Fig. 5 and a time step of 0.0005 sec (T 1 /50) has been obtained after two iterations. The maximum Von-Mises stress for the initial shape corresponding to the refined mesh (Fig. 5) has been obtained as 199 MPa. The control of error in both space and time is very important, since constraints are imposed on stresses during the course of optimization. If we take the coarse mesh (Fig. 4) and the initial coarse time step 1.667e-3 sec, the maximum Von-Mises stresses is obtained as 188 Mpa. Even if we take a finer time step of 0.0005 sec corresponding to the same coarse mesh (Fig. 4), the maximum Von-Mises stress still remains underestimated. Thus, a combination of appropriate mesh and time step size will only ensure accuracy to the finite element solution and subsequent evaluation of maximum stresses. Since the evaluation of optimal result using GA will largely depend on the values of stresses as constraint, the error in both space and time should be under control from the beginning of optimization. Over the entire time domain, the maximum Von-Mises stress value among all the nodes considered are evaluated and a stress constraint is imposed such that the maximum stress values anywhere should not exceed permissible value of 250 MPa. Additional constraint like y 1 > y 2 > y 3 > y 4 > y 5 > y 6 are also introduced from the point of view of bending stress distribution. Since the maximum stress constraint value is 250 MPa, which is more than the maximum Von-Mises stress corresponding to the initial shape, scope for reduction in surface area of the problem domain exists during optimization iterations. Table 1 contains the genetic parameters and operators used during the optimization of the cantilever. Three configurations (by varying population, probability of crossover and mutation) have been tried, which are designated as Run 1, Run 2 and Run 3. The converged optimized shape of the cantilever is shown in Fig. 6 corresponding to configuration Run1. Figure 7 shows the plot for maximum Von-Mises stresses vs. time steps corresponding to different updated geometry of the cantilever beam obtained from modified shape controlling variables during the optimization process. It may be noted that the surface area has reduced corresponding to the increase in stress values. For optimal surface area, the effective Von-Mises stress in the cantilever has reached to 250 MPa, which is the permissible limit. The final area (best optimized) is 0.734 m 2 , which indicates a reduction of 38.83% reduction. It may be noted that the average and worst optimized areas are 0.814 m 2 and 0.922 m 2 respectively. The error in space and time corresponding to the converged optimal shape has been observed to be quite close to the prescribed accuracy limit. Figure 8 shows the optimization history in terms of non-dimensional area factor (both best optimized area/initial area as well as average optimized area/initial area) and stress factor (maximum stress/permissible stress). It is observed that the non-dimensional area and stress factor are always below one, which clearly indicates the stability of the optimization process. Moreover, the stable nature of these curves around the optimal value is a clear indication of the good convergence of the genetic algorithm. Further, from Fig. 9 it can also be observed that the values of all the design variables are also stable, which indicates global convergence of the optimal solution. Figure 10 shows the plot of non-dimensional area factor (best optimized area/initial area) corresponding to different genetic parameters as detailed in Table 1. It is observed that the converged optimal values of the surface area corresponding to three different combinations of genetic parameters are more or less same, indicating the global nature of the convergence. The optimization process has been carried out up to 150 generation for convergence study. However, the converged result has been obtained at 100 generation and the total process up to 100 generation has taken about 11 hours, including IO operations in a 533 MHz Pentium-Pc based computer.

Simply supported beam subjected to transient step loading at middle
The initial geometry with dimensions, boundary conditions, loading and prescribed design variables of a simply supported beam is shown in Fig. 11. The objective of this exercise is also to find optimal surface area provided that the Von-Mises stress anywhere within the domain does not exceed 250 MPa evaluated over the prescribed time domain. External boundary AB has been modeled using B-spline (Fig. 11). Seven control vertices define the spline, which undergoes shape changes and Y coordinates of these points have been considered as design variables (y 1 , y 2 , . . . , y 7 ). 8 noded isoparametric quadrilateral elements with 2 × 2 Gauss quadrature is chosen for finite element analysis. 24 elements having an initial area of 1.2 m 2 define the finite element model as shown in Fig. 12

Generation Stress & Area factor
Best area factor Average area factor Stress factor  magnitude 20 kN for a total time of 0.01 sec. applied at the middle of the beam has been considered as shown in Fig. 13. A time step of magnitude 5e-4 sec (T 1 /10) has been chosen for time integration using Newmark-β integration method (constant average acceleration). In order to ensure the accuracy of dynamic response, both space and time discretization errors are evaluated and have been obtained as 13.62% and 16.86% respectively. For a target accuracy of 3% and 1% for space and time discretization, a finite element mesh having 72 elements as shown in Fig. 14     Step loading for simply supported beam. a time step of 0.0001 sec (T 1 /50) has been obtained after two iterations. Von-Mises stresses on all nodes for all the time steps are calculated corresponding to the refined mesh. The maximum Von-Mises stress for initial shape has been obtained as 204 MPa corresponding to the refined mesh (Fig. 14). Since the maximum permissible stress is 250 MPa, scope for reduction of surface area with consequent increase in stress value exists during optimization  process. Side constraints like y 1 < y 2 < y 3 < y 4 > y 5 > y 6 > y 7 are imposed keeping in mind the bending stress distribution in a simply supported beam for the applied loading. The final optimized geometry of the beam should be symmetric because of loading pattern. In order to ensure a symmetric optimal design, additional equality constraints like y 1 = y 7 ; y 2 = y 6 ; y 3 = y 4 have been introduced. Table 2 contains the genetic parameters and operators used during the optimization of the cantilever. Three configurations (by varying population, probability of crossover and mutation) have been tried, which are designated as Run 1, Run 2 and Run 3. The converged optimized shape of the simply supported beam is shown in Fig. 15 corresponding to configuration Run 1. Figure 16 shows the maximum Von-Mises stresses vs. time steps plot for the different updated geometry of the simply supported beam corresponding to modified shape controlling variables in the optimization process. For optimal surface area the effective Von-Mises stress in the simply supported beam has reached to 248 MPa. The final area (best optimized) is 1.072 m 2 , which means a reduction of 10.7%. It may be noted that the average and worst optimized areas are 1.106 m 2 and 1.151 m 2 respectively. The error in space and time corresponding to the converged optimal shape has been also observed to be quite close to the prescribed accuracy limit. Figure 17 shows the optimization history in terms of non-dimensional area factor (both best optimized area/initial area as well as average optimized area/initial area) and stress factor (maximum stress/permissible stress). The stability of these curves at the end of the process has been observed, which is a clear indication of the convergence of the optimization process. Figure 18 shows design variables vs. number of generations plot and steady values of the variables are observed after 60 generations, which also indicates the global convergence of the search process. Figure 19 shows the plot of non-dimensional area factor (best optimized area/initial area) corresponding to different genetic parameters as detailed in Table 2. It is observed that the converged optimal values of the surface area corresponding to three different combinations of genetic parameters are nearly similar, indicating the global nature of the convergence. The optimization process has been carried out up to 100 generation for convergence study. However, the converged result has been obtained at 60 generation and the total process up to 60 generation has taken about 7 hours, including IO operations in a 533 MHz Pentium-Pc based computer. Time steps

Stress (Pas)
Initial geometry Partially optimized geometry corresponding to 20 gen. of GA Optimized geometry corresponding to 60 gen. of GA

Conclusion
In the present study of simulation of optimal shape of structures subjected to transient dynamic loading, an integrated approach combining geometric modeling, reliable dynamic finite element analysis using implicit Newmark-β integration method and a robust optimization tool i.e. Genetic algorithm has been presented. B-spline curves have been used to model the changing boundary of the structure and an automatic mesh generator has been used to generate the analytical model from the geometric model through a set of geometric parameters. As shown in the numerical examples, use of this integrated approach leads to evaluation of reliable optimal shape satisfying all imposed constraints and also demonstrates that this integrated module is capable to deal with general shape optimization problems. The use of Genetic algorithms as optimization tool has helped to avoid the evaluation of computing intensive and involved design sensitivity calculations, which is common in shape optimization using traditional mathematical programming techniques. However, because of population approach GAs are fairly expensive in terms of computational time.