The development of a shape optimization solver using the existing Computational Fluid Dynamics (CFD) codes is taken up as topic of research in this thesis. A shape optimizer was initially developed based on Genetic Algorithm (GA) coupled with a CFD solver in an earlier work. The existing CFD solver is based on Kinetic Flux Vector Splitting and uses least squares discretization. This solver requires a cloud of points and their connectivity set, hence this CFD solver is a meshless solver. The advantage of a meshless solver is utilised in avoiding re-gridding (only connectivity regeneration is required) after each shape change by the shape optimizer. The CFD solver is within the optimization loop, hence evaluation of CFD solver after each shape change is mandatory. Although the earlier shape optimizer developed was found to be robust, but it was taking enoromous amount of time to converge to the optimum solution (details in Appendix). Hence a new evolving method, Ant Colony Optimization (ACO), is implemented to replace GA. A shape optimizer is developed coupling ACO and the meshless CFD solver. To the best of the knowledge of the present author, this is the first time when ACO is implemented for aerodynamic shape optimization problems. Hence, an exhaustive validation has become mandatory. Various test cases such as regeneration problems of (1) subsonic - supersonic nozzle with a shock in quasi - one dimensional flow (2) subsonic - supersonic nozzle in a 2-dimensional flow field (3) NACA 0012 airfoil in 2-dimensional flow and (4) NACA 4412 airfoil in 2-dimensional flow have been successfully demonstrated. A comparative study between GA and ACO at algorithm level is performed using the travelling salesman problem (TSP). A comparative study between the two shape optimizers developed, i.e., GA-CFD and ACO-CFD is carried out using regeneration test case of NACA 4412 airfoil in 2-dimensional flow. GA-CFD performs better in the initial phase of optimization and ACO-CFD performs better in the later stage. We have combined both the approaches to develop a hybrid GA-ACO-CFD solver such that the advantages of both GA-CFD and ACO-CFD are retained with the hybrid method. This hybrid approach has 2 stages, namely, (Stage 1) initial optimum search by GA-CFD (coarse search), the best members from the optimized solution from GA-CFD are segregated to form the input for the fine search by ACO-CFD and (Stage 2) final optimum search by ACO-CFD (fine search). It is observed that this hybrid method performs better than either GA-CFD or ACO- CFD, i.e., hybrid method attains better optimum in less number of CFD calls. This hybrid method is applied to the following test cases: (1) regeneration of subsonic-supersonic nozzle with shock in quasi 1-D flow and (2) regeneration of NACA 4412 airfoil in 2-dimensional flow. Two applications on shape optimization, namely, (1) shape optimization of a body in strongly rotating viscous flow and (2) shape optimization of a body in supersonic flow such that it enhances separation of binary species, have been successfully demonstrated using the hybrid GA-ACO-CFD method. A KFVS based binary diffusion solver was developed and validated for this purpose. This hybrid method is now in a state where industrial shape optimization applications can be handled confidently.