Consider a small volume element
$$\{(s_1,s_2,s_3)\in \mathbb{R}^3: x_1 \leq s_1 \leq x_1 + \Delta x_1,x_2 \leq s_2 \leq x_2 + \Delta x_2,x_3 \leq s_3 \leq x_3 + \Delta x_3\}.$$
The surface (stress) force acting on that element in the $i$ direction consists of contributions from each of three pairs of faces where coordinates $s_1$, $s_2$ and $s_3$ are constant, respectively.
Consider one pair of where $s_1$ is fixed: $s_1 = x_1$ and $s_1 = x_1 + \Delta x_1$. The net force in the $i$ direction is
$$F_{i,1}\approx \sum_{j=1}^{3}\sigma_{ij}n_j|_{s_1=x_1 + \Delta x_1}\Delta x_2\Delta x_3+\sum_{j=1}^{3}\sigma_{ij}n_j|_{s_1=x_1}\Delta x_2\Delta x_3$$
This is approximate because we evaluate $\sigma_{ij}(s_1,s_2,s_3)$ at representative coordinates $s_2= x_2$ and $s_3=x_3$ . In general, we assume that the stress tensor is a continuously differentiable function of spatial position.
Here we use the fact that the surface force (in terms of the stress tensor and normal vector) is the product $\mathbb{\sigma}\cdot\mathbb{n}.$ On the face $s_1 = x_1 + \Delta x_1$ we have $\mathbb{n} = (1,0,0)$, and on the face $s_1 = x_1$ we have $\mathbb{n} = (-1,0,0).$
Hence,
$$F_{i,1}\approx \sigma_{i1}|_{s_1=x_1 + \Delta x_1}\Delta x_2\Delta x_3-\sigma_{i1}|_{s_1=x_1}\Delta x_2\Delta x_3$$
The contribution to the force per unit volume is
$$\frac{F_{i,1}}{\Delta x_1 \Delta x_2 \Delta x_3}\approx \frac{\sigma_{i1}|_{s_1=x_1 + \Delta x_1}-\sigma_{i1}|_{s_1=x_1}}{\Delta x_1}.$$
Now consider the limit as the element shrinks to the point $(x_1,x_2,x_3)$ and $\Delta x_k \rightarrow 0.$
Then the contribution to the net surface force per unit volume at that point is
$$\hat{F_{i,1}}= \lim_{\Delta x_1 \rightarrow 0}\frac{\sigma_{i1}|_{s_1=x_1 + \Delta x_1}-\sigma_{i1}|_{s_1=x_1}}{\Delta x_1}=\frac{\partial \sigma_{i1}}{\partial x_1}.$$
Applying the same derivation to the other pairs of faces we obtain the total contribution to the net surface force per unit volume (in the $i$ direction) as
$$\hat{F_{i}}=\sum_{j=1}^{3}\frac{\partial \sigma_{ij}}{\partial x_j}$$
You are conflating incompressible with inviscid. Here we have flow of an incompressible viscous fluid and the term $\nu \nabla^2 \mathbf{v}$ may not be neglected.
You are asked to find "a solution" under the assumptions that the flow is steady and uniform in the x-direction. This implies fully-developed unidirectional flow where the only non-vanishing component of velocity is $u$ (i.e., the component in the x-direction).
For an incompressible fluid, the equation of continuity gives
$$\nabla \cdot \mathbf{v} = \frac{\partial u}{\partial x} = 0$$
Thus, $u$ is a function only of the y-coordinate and the Navier-Stokes equations reduce to
$$0 = -\frac{1}{\rho}\frac{\partial p}{\partial x} +g \sin \alpha + \nu \frac{d^2u}{dy^2}, \\0 = -\frac{1}{\rho}\frac{\partial p}{\partial y} -g \cos \alpha $$
Since the flow is gravity-driven, we can neglect the x-component of the pressure pressure gradient $\frac{\partial p}{\partial x}$ in the first equation and find $u$ by applying the two boundary conditions to the solution of the second-order differential equation
$$\nu \frac{d^2 u}{dy^2} = -g \sin \alpha$$
The second equation allows for solution of the pressure as a function of the y-coordinate.
See if you can finish now.
Best Answer
We can see how the conjugate arises by deriving the Blasius formula, where to be perfectly clear, the squared derivative of the potential $w$ should appear -- not the squared modulus.
In inviscid flow, the fluid velocity is tangent to the boundary $\partial B$ which is represented as a closed contour in the complex plane. Let $\alpha$ be the angle between the tangent and the real axis. The normal vector directed into the body is then, $\mathbb{n} = -\sin \alpha + i\cos \alpha = i e^{i \alpha}.$
The components of velocity are related to the complex potential by
$$\frac{dw}{dz} = -u_x + i u_y = \sqrt{u_x^2 + u_y^2}(- \cos \alpha + i \sin \alpha) = -\sqrt{u_x^2 + u_y^2}e^{-i \alpha}.$$
The pressure is given by Bernoulli's equation
$$p = p_0 - \frac{\rho}{2}\left( u_x^2 + u_y^2\right) = -\frac{\rho}{2}\left( u_x^2 + u_y^2\right),$$
where we take the constant $p_0=0$, since there is no net contribution when integrated around a closed contour.
Hence,
$$p = -\frac{\rho}{2}\left( \frac{dw}{dz}\right)^2e^{2i\alpha}.$$
Note that the above expression involves the squared derivative of the potential, not the squared modulus of the derivative.
The complex representation of the force on the body is obtained as the integral with respect to arclength $s$ around the contour.
$$F = F_x + i F_y = \int_{\partial B} p\, n \, ds = \int_{\partial B} p\, ie^{i \alpha} \, ds. $$
We take the conjugate in order to get a cancellation of the complex exponential when we later substitute for pressure:
$$\overline{F} = -i\int_{\partial B} p\, e^{-i \alpha} \, ds.$$
Note that the pressure is real-valued and is not affected upon taking the conjugate.
The differential of the complex variable is related to arclength differential by $dz = e^{i\alpha} ds.$
Hence,
$$\overline{F} = -i\int_{\partial B} p\, e^{-2i \alpha} \, dz.$$
Substituting for pressure , we get
$$\overline{F} = F_x - iF_y = -i\int_{\partial B} \frac{\rho}{2}\left( \frac{dw}{dz}\right)^2e^{2i\alpha}\, e^{-2i \alpha} \, dz \\ = -\frac{i\rho}{2}\int_{\partial B} \left( \frac{dw}{dz}\right)^2 \, dz .$$