The equations you propose remind me of CFD. Mathematica (or any other similar system) is not well suited to solving CFD or the equations you propose; rather, such tools are best used for equations that have an analytical solution or a simple numerical one. Complex equations with even more complex boundary conditions (such as is the case with CFD) have led to a wide variety of approaches, each created to suit a limited subset of the total possible configurations. For example, the FDM might work quite well for flow through a reservoir, while the FEM might work much better for flow over a complex shape such as a Porsche. I doubt this approach will yield satisfactory results.

why you are not applying the newly developed Coiflet wavelet Homotopy analysis method to solve these PDEs. for details please read the work of Professor Liao and Yu Qiang Article Analysis of mixed convection flow in an inclined lid-driven ...