setupBPT[terms_List]:=With[{ud=Union[Flatten[UnitDimensions/@terms,1][[All,1]]]},If[MemberQ[ud,_UnitDimensions],Throw[None]];{Length@terms-Length[ud],ud}]
UDToExp[dimensions_List]:=Times@@(dimensions/.List[d_String,n_]:>Power[d,n])
Nondimensionalize[expr_,vars_List]:=With[{exp=expr/.{Power[q_Quantity,n_]:>UDToExp[UnitDimensions[q]]^n,qu_Quantity:>UDToExp[UnitDimensions[qu]]}},With[{res=Solve[Thread[((List@@PowerExpand[exp])/.Power[unit_String,n_]:>n)==0],vars]},If[Flatten[res]==={},Throw[None],With[{nondim=expr/.q_Quantity:>QuantityMagnitude[q]},First[nondim/.res]]]]]
Buckinghamπ[terms_List,const:{__Integer}]:=Catch[With[{res=setupBPT[terms]},If[Length[Union[const]]=!=First[res],Throw[None]];If[MemberQ[const,x_/;x>Length[terms]],Throw[None]];With[{c=#},Module[{temp,vars={}},Check[temp=If[Not[MemberQ[const,#]],With[{u=Unique[]},AppendTo[vars,u];terms[[#]]^u],If[c=!=#,1,terms[[#]]]]&/@Range[Length[terms]],Throw[None],{Part::partw}];Nondimensionalize[Times@@temp,vars]]]]&/@const]

BPTMatrixLabels={"pressure","dynamic viscosity","velocity","characteristic length","mass density"};
P=Quantity[p,"Newtons"/"Meters"^2];
mu=Quantity[μ,("Pascals"*"Seconds")];
V=Quantity[v,"Meters"/"Seconds"];
L=Quantity[d,"Meters"];
Rho=Quantity[ρ,"Kilograms"/"Meters"^3];
terms={P,mu,V,L,Rho};

MatrixForm[Outer[Buckinghamπ[terms,{#1,#2}]&,Range[5],Range[5]],TableHeadings->{BPTMatrixLabels,BPTMatrixLabels}]