Micropolar fluids are fluids with microstructures. They
belong to a class of fluids with a nonsymmetric stress tensor.
Micropolar fluids consist of rigid, randomly oriented
or spherical particles with their own spins and microrotations,
suspended in a viscous medium.
Physical examples of micropolar fluids can be seen in
ferrofluids, blood flows, bubbly liquids, liquid
crystals, and so on, all of them containing intrinsic polarities.
The following (including notations) is based on textbook
G. Lukaszewicz, Micropolar Fluids: Theory and Applications, Birkhauser, Boston, 1999
http://books.google.ru/books?id=T3l9cGfR9o8C
We start with Cauchy momentum equation
$$
\rho \frac{D\vec{v}}{Dt} = \rho \vec{f} + \nabla \cdot \hat{T},
$$
where $\vec{f}$ is body force and $\hat{T}$ is stress tensor.
If we assume that fluid is polar, that is it has its own internal angular momentum (independent of the motion of fluid as a whole) then we need an additional
equation expressing conservation of angular momentum (for nonpolar fluid conservation of angular momentum is a consequence of Cauchy equation):
$$
\rho \frac{D}{Dt}(\vec{l} + \vec{x}\times\vec{v}) = \rho \vec{x} \times \vec{f} + \rho\vec{g} + \nabla \cdot (\vec{x} \times \hat{T} + \hat{C}).
$$
Here $\vec{l}$ is an intrinsic (internal) angular momentum per unit mass, $\vec{g}$ is body torque and $\hat{C}$ is a new object called couple stress tensor
(this equation could be seen as its definition).
Now in order to close this system of equation we need to express the stress tensor and couple stress tensor through the characteristics of fluid dynamics.
For this we need to make certain assumptions: the absence of preferred direction and position, reduction to hydrostatic pressure in the absence of deformations and
linear dependence on the velocity spatial derivatives $v_{i,j}$ (deformation). For polar fluid we also define the vector field $\vec{\omega}$ -- microrotation
which represents the angular velocity of rotation of particles of the fluid. We further assume that the fluid is isotropic and $\vec{l} = I \,\vec{\omega}$ with $I$ a scalar called the microinertia coefficient. Couple stress tensor should be a linear function of the spatial derivatives of the microrotation field: $\omega_{i,j}$. All this assumptions allow us to specify the general form for stress tensor:
$$
T_{ij}=(- p +\lambda v_{k,k})\,\delta_{ij}+\mu\,(v_{i,j}+v_{j,i})+\mu_r\,(v_{j,i}-v_{i,j}) - 2 \mu_r\, \epsilon_{mij}\omega_m,
$$
and couple stress tensor:
$$
C_{ij} = c_0\, \omega_{k,k} \delta_{ij} + c_d\, (\omega_{i,j}+\omega_{j,i}) + c_a \,(\omega_{j,i}-\omega_{i,j}).
$$
Note: we have three parameters $c_0$, $c_d$ and $c_a$ (called coefficients of angular
viscosities) because there are three irreducible representations for the action of $SO(3)$ group on rank 2 tensor: scalar (times $\delta_{ij}$), traceless symmetric tensor, antisymmetric tensor. The same for the triplet $(\lambda,\mu, \mu_r)$ (which are called second viscosity coefficient, dynamic
Newtonian viscosity and dynamic microrotation viscosity).
Substituting this expressions for $T_{ij}$ and $C_{ij}$ into Cauchy equation and angular momentum equation we obtain the equations from the question.
The notations correspondence is: $$ \vec{\omega} \to \vec{N}^{*}, \qquad \mu_r \to k_1^{*}, \qquad I\to j^{*} \qquad c_a+c_d \to \gamma^{*}.$$
Additionally we see that the equations in the question assume $\mathrm{div} \vec{v} =0 $ (incompressible flow) and $\mathrm{div} \vec{\omega} =0 $ -- this
one is generally not true, but symmetries of the system could make it so. Also we see that there are no body torque and the body force term corresponds to Lorentz force.
Let's suppose that the boundary is the x-axis. So along the boundary, the stream function is constant. So, $$\frac{\partial^2\psi}{\partial x^2}=0$$ And from the no-slip boundary condition, $$\frac{\partial \psi}{\partial y}=0$$ This can be used to establish the 2nd order finite difference approximation to the value of the vorticity at the boundary:$$\omega=\frac{\partial^2 \psi}{\partial x^2}+\frac{\partial^2 \psi}{\partial y^2}=\frac{2(\psi(I,1)-\psi(I,0))}{(\Delta y)^2}=\omega(I,0)$$
This is used for the boundary condition on $\omega$.
ADDENDUM
First of all, the stream function is known (and constant) at the solid boundaries (because the solid boundaries are stream lines). So it doesn't have to be solved for. It is determined up to an arbitrary constant, and can thus be taken to be zero at one of the boundaries. At the other boundary, the stream function is equal to the volumetric throughput rate per unit width of channel (which is typically known). So you don't need to solve for the stream function at the solid boundaries.
If the tangential derivative $\partial \psi/\partial x=0$ at all locations along the boundary, it's second partial with respect to x must also be equal to zero. This, of course, all follows from the fact that $\psi$ is constant at the boundary.
To integrate the vorticity equation, you need a boundary condition on the vorticity (or at least a 2nd order finite difference approximation to a boundary condition). Just because $\partial \psi/\partial y=0$ does not mean the the second partial of $\psi$ with respect to y is equal to zero at the boundary; this would imply that the vorticity at the boundary is equal to zero, which we know is not correct.
The variable I in the relationships refers to the I'th x grid point. So, back to the boundary condition on vorticity: We have shown so far that, at the boundary, $$\omega=\frac{\partial^2 \psi}{\partial y^2}$$ subject to the constraint that $\partial \psi/\partial y=0$. If we represent these two conditions in 2nd order finite difference form, we obtain:
$$\omega(I,0)=\frac{\psi(I,1)-2\psi(I,0)+\psi(I,-1)}{(\Delta y)^2}$$and$$\frac{\psi(I,1)-\psi(I,-1)}{2\Delta y}=0$$If we combine these two finite difference equations, we obtain a 2nd order finite difference approximation to the value of the vorticity at the boundary:
$$\omega(I,0)=\frac{2(\psi(I,1)-\psi(I,0))}{(\Delta y)^2}$$
I've successfully used this approach to solving these equations many times.
Best Answer
I don't know a good answer to your first question (I'd be interested in a good text for that myself), but I can answer the second.
It's easier to explain if we temporarily imagine $\phi$ represents the concentration of some dye made up of little particles suspended in the fluid. The convective term (aka advective term) is transport of $\phi$ due to the fact that the fluid is moving: a single "particle" of $\phi$ will tend move around according to the velocity of the fluid around it. The diffusive term, on the other hand, represents the fact that the dye tends to spread out, regardless of the motion of the fluid, because each particle is undergoing Brownian motion. So if you were moving along at the same velocity as the fluid you would see a small spot of dye tend to become more and more blurred over time.
For quantities like energy and momentum the diffusion happens for a slightly different reason (transfer of the quantity between fluid molecules when they collide) but the principle is the same. The property is transported along with the fluid's bulk velocity (convective term) but also tends to spread out and become blurred of its own accord (diffusive term).