The aim of the current work is to simulate the laminar flow over square cylinder by developing numerical code using C++ programing language. The findings of this simulation are streamlines and isotherm pattern and global quantities like; lift force coefficient, Drag force coefficient, recirculation length, and Nusselt number. Flow over square cylinder is an important research topic in many engineering fields due to its numerous applications such as building aerodynamics, cooling of electronic component and compact heat exchangers. The thesis is divided into 6 chapters. The objectives and scope of the work is explained in chapter1. In chapter 2, the governing equation and the assumption used to simplify these equations are discussed. Discretization of the domain and the governing equation will be carried out in chapter 3. In chapter 4, the algorithm for pressure velocity coupling (fractional step Method), linear solver used for solving system of linear equations and different boundary conditions are discussed. The developed numerical code is verified in chapter 5 using Method of Manufactured Solution (MMS) and comparing code solution to high accurate Benchmark solution of the Driven cavity and differential heated square cavity. Finally in chapter 6, the simulation of laminar flow over square cylinder is explained, in which the effect of Reynolds number and effect of blockage ratio are investigated. Also tandem arrangement for two square cylinder and effect of adding gravity term using Boussinesq approximation for this arrangement are studied. It is concluded that, flow patterns, drag, lift coefficient, and strouhal number are affected by changing Reynolds number and blockage number. In case of tandem arrangement, downstream cylinder has different drag, lift, and Nusselt number from upstream cylinder. Moreover adding buoyancy affect different calculated parameter especially lift coefficient.