Could you visualize solution without numerical artefacts by adding options PlotPoints -> 100 and Chop?

@AlexTrounev. MWE updated as you have proposed.

Row 1 and 2 are very nice, but Row 3 pictures look wrong. Please, use option PlotPoints -> 100 without Chop for plt3DfluxLC and plt3Dflux.

@AlexTrounev. I removed Chop@ but didn't notice any difference. Then I added "MaxRecursion->4" and again I didn't notice any change. The updated row 3 is now under the old row 3. Please note that radial flux oscillates around 0 near zero of velosity u.

Try DensityPlot[ Norm[ flux[Sqrt[vx^2 + vy^2], vx/Sqrt[vx^2 + vy^2], {\[Mu], \[Eta]}]] /. F -> fs, {vx, umin, umax}, {vy, umin, umax} , AxesLabel -> {"\!\(\*SubscriptBox[\(v\), \(z\)]\)", "\!\(\*SubscriptBox[\(v\), \(\[UpTee]\)]\)"}, PlotLabel -> lbl, PlotLegends -> Automatic , ColorFunction -> Hue , ImageSize -> iSize, PlotRange -> Automatic, PlotPoints -> 100];

@AlexTrounev: Done, Note that PlotRange was changed to Automatic.

Please, use exactly what I proposed {vx, umin, umax}, {vy, umin, umax}

Let us continue this discussion in chat.

To me it's not exactly clear the what the problem is and especially what you expect.

@user21: I expect that radial flux near zero velocity, at u -> umin tends to zero. From physics of the problem this would mean that the solution tends to Maxwellian distribution function. I expected that boundary condition NeumannValue[0, (0 <= \[Xi] <= 1) && (u == umin)] provides zeroing of the redial flux. But this is not the case. Best of all it is seen from the right picture in the last row. This picture shows that solution is isotropic in the vicinity of u=umin as expected [we see no dependence on [Xi]]. However it shows thar radial flux does not tend to zero.

I plotted the radial flow versus speed at a fixed angle equal to the angle of the source. It can be seen from the graph that the flux decreases linearly from the maximum at the source to zero at zero velocity. However, near zero velocity, a numerical instability appears in the form of an oscillator. The amplitude and frequency of oscillations decreased after changing the variable u -> w=u^2. Thus, the question of how to correctly set the BC must be recognized as obsolete. The boundary condition works as expected. But the question arises how to prevent instability.

I still do not understand anything, but have you tried to refine the mesh in that area to at least contain the oscillations?

@user21. Yes. here are mesh for some of recent runs: NDSolveFEMQuadElement[<3611>] NDSolveFEMQuadElement[<18515>] NDSolveFEMElementMesh[{{4.*10^-6,36.},{0.,1.}},{NDSolveFEMQu‌​adElement[<18515>]}] Perhaps, it is worth to refine mesh locally near u=0? but I don't find how to do that,

@user21 I am trying to invoke ToGradedMesh function introduced in v.13.

I implemented non-uniform mesh: