This thesis presents a new Barely Implicit Correction (BIC) algorithm combined with a modified flux-corrected transport (FCT) algorithm for the simulation of three-dimensional (3D), low-Mach-number flows and then proceeds to apply it to the study of vortex breakdown undergoing heat addition and heat extraction. This new algorithm is based on the original, introduced by a prior work in 1987, which was a solution procedure including an explicit predictor step to solve the convective portion of the Navier-Stokes equations and an implicit corrector step to remove the acoustic limit on the integration time-step. The explicit predictor uses the flux-corrected transport (FCT) algorithm while the implicit corrector solves an elliptic equation for a pressure correction to equilibrate acoustic waves. This thesis introduces a procedure for stabilizing and implementing FCT for 3D flows and extends BIC for 3D with physical diffusion processes. A new filter is introduced to further stabilize the algorithm and the solution procedure is clarified for the inclusion of the diffusion fluxes. The new BIC-FCT algorithm is examined in four test problems with successively increased difficulty. The test problems culminate with calculations of vortex breakdown in 3D swirling flows. All the test problems demonstrate that the algorithm is able to predict accurate and robust solutions using time steps varying from near the explicit stability limit to tens and hundreds of times larger. Excellent agreement is also obtained when compared with results from other algorithms. The algorithm is then used to study how vortex breakdown is affected when heat is extracted from or added to the flow. Two heat release rates are applied to a flow with a bubble mode of breakdown upstream and double-helix mode downstream. The simulations show that heat release causes the double-helix structure to become narrower. With more heat release, the double-helix mode transitions to a columnar vortex. In addition, a lower heat extraction rate causes the columnar vortex to first transition to a spiral mode and then to a double-helix mode. With a higher heat extraction rate, the columnar vortex transitions to a double-helix mode, bypassing the spiral mode. Further investigation show that the density gradient formed by heat addition and extraction is the dominant effect in the transitions. The transition is promoted by changes in viscosity due to temperature changes from heat addition and extraction. The new algorithm presented in this thesis provides a new way to calculate low-Mach-number flows. Such vortex breakdown simulations with heat changes serve as a base for understanding the dynamics of a precessing vortex core in swirl combustors and other vortex flows with changes in heat input.