Needs["NDSolve`FEM`"]Ω=Rectangle[{0, 0}, {1, 1}];op1=
(0,1)
v
[x,y] +
(1,0)
u
[x,y];op2={u[x, y], v[x, y]}. Inactive[Grad][u[x, y], {x, y}] +Inactive[Div][{p[x, y],0*p[x,y]}, {x,y}] -Pr*Inactive[Div][{
(1,0)
u
[x,y],
(0,1)
u
[x, y]}, {x, y}];op3={u[x, y], v[x, y]} . Inactive[Grad][v[x, y], {x, y}] +Inactive[Div][ {0*p[x,y],p[x, y]},{x,y}] -Pr*Inactive[Div][{
(1,0)
v
[x, y],
(0,1)
v
[x, y]}, {x, y}]-Pr*Ra*θ[x,y];op4={u[x, y], v[x, y]} . Inactive[Grad][θ[x, y], {x, y}] - Inactive[Div][{
(1,0)
θ
[x, y] ,
(0,1)
θ
[x, y]}, {x, y}];Pr=1;Ra=10^5;OP={op1==0,op2==0,op3==0,op4==0};bc={DirichletCondition[θ[x,y]==1,y==0],DirichletCondition[θ[x,y]==0.,x==1],DirichletCondition[θ[x,y]==0.,x<1],DirichletCondition[p[x,y]==0,x==0 && y==0],DirichletCondition[{u[x,y]==0.,v[x,y]==0.},True]}; NDSolveValue[{OP, bc},{u, v, p, θ},{x,0,1},{ y,0,1} ,Method  {"FiniteElement","InterpolationOrder"  {u  2, v  2, p  2, θ  1}}]