The question is:
A half cylinder with the square part on the $xy$-plane, and the length $h$ parallel to the $x$-axis. The position of the center of the square part on the $xy$-plane is $(x,y)=(0,0)$.
$S_1$ is the curved portion of the half-cylinder $z=(r^2-y^2)^{1/2}$ of length $h$.
$S_2$ and $S_3$ are the two semicircular plane end pieces.
$S_4$ is the rectangular portion of the $xy$-plane
Gauss' law:
$$\iint_S\mathbf E\cdot \mathbf{\hat{n}}\,dS=\frac{q}{\epsilon_0}$$
$\mathbf E$ is the electric field $\left(\frac{\text{Newton}}{\text{Coulomb}}\right)$.
$\mathbf{\hat{n}}$ is the unit normal vector.
$dS$ is an increment of the surface area $\left(\text{meter}^2\right)$.
$q$ is the total charge enclosed by the half-cylinder $\left(\text{Coulomb}\right)$.
$\epsilon_0$ is the permitivity of free space, a constant equal to $8.854\times10^{-12}\,\frac{\text{Coulomb}^2}{\text{Newton}\,\text{meter}^2}$.
The electrostatic field is:
$$\mathbf{E}=\lambda(x\mathbf{i}+y\mathbf{j})\;\text{,}$$
where $\lambda$ is a constant.
Use this formula to calculate the part of the total charge $q$ for the curved portion $S_1$ of the half-cylinder:
$$\iint_S\mathbf E\cdot \mathbf{\hat{n}}\,dS=\frac{q}{\epsilon_0}=\iint_R\left\{-E_x[x,y,f(x,y)]\frac{\partial f}{\partial x} -E_y[x,y,f(x,y)]\frac{\partial f}{\partial y} +E_z[x,y,f(x,y)] \right\}\,dx\,dy$$
The goal is to find the total charge $q$ enclosed by the half-cylinder, expressed in terms of $\lambda$, $r$ and $h$.
The solution should be:
$$\pi r^2\lambda h\epsilon_0$$
This is what I've tried:
First calculate Gauss' law for $S_1$:
\begin{align}
f(x,y)&=z=(r^2-y^2)^{1/2}=\sqrt{(r^2-y^2)} \\
\frac{\partial f}{\partial x}&=\frac12(r^2-y^2)^{-\frac12}\cdot 0=0 \\
\frac{\partial f}{\partial y}&=\frac12(r^2-y^2)^{-\frac12}\cdot -2y=-\frac{y}{\sqrt{(r^2-y^2)}}=-\frac yz \\
\\
\mathbf{E}&=\lambda(x\mathbf{i}+y\mathbf{j}) \\
E_x[x,y,f(x,y)]&=\lambda x \\
E_y[x,y,f(x,y)]&=\lambda y \\
E_z[x,y,f(x,y)]&=0 \\
\\
\text{length}&=h \\
\end{align}
Using the formula
$$\iint_R\left\{-E_x[x,y,f(x,y)]\frac{\partial f}{\partial x} -E_y[x,y,f(x,y)]\frac{\partial f}{\partial y} +E_z[x,y,f(x,y)] \right\}\,dx\,dy$$
we get:
\begin{align}
&\iint_R\left\{-\lambda x\cdot 0-\lambda y\cdot -\frac{y}{z} + 0\right\}\,dx\,dy \\
&=\iint_R\frac{\lambda y^2}{z}\,dx\,dy \\
&=\lambda\iint_R\frac{y^2}{\sqrt{r^2-y^2}}\,dx\,dy \\
\end{align}
Since the length is $h$ and the length is parallel to the $x$-axis:
\begin{align}
&\lambda \int_R\int_0^h\frac{y^2}{\sqrt{r^2-y^2}}\,dx\,dy \\
&=\lambda\int_R\left[\frac{y^2x}{\sqrt{r^2-y^2}}\right]_0^h\,dy \\
&=\lambda\int_R\frac{y^2h}{\sqrt{r^2-y^2}}\,dy \\
\end{align}
Subsitute:
\begin{align}
y&=r\sin\theta \\
\theta&=\arcsin\left(\frac1r y\right) \\
\frac{dy}{d\theta}&=\frac{d}{d\theta}\left(r\sin\theta\right)=r\cos\theta \\
dy&=r\cos(\theta)\,d\theta \\
\\
&\lambda\int\frac{hr^2\sin^2\theta}{\sqrt{r^2-r^2\sin^2\theta}}\cdot r\cos(\theta)\,d\theta \\
&=\lambda h\int\frac{r^3\sin^2\theta\cos\theta}{r\sqrt{1-\sin^2\theta}}\,d\theta \\
&=\lambda hr^2\int\frac{\sin^2\theta\cos\theta}{\sqrt{\cos^2\theta}}\,d\theta \\
&=\lambda hr^2\int\frac{\sin^2\theta\cos\theta}{\cos\theta}\,d\theta \\
&=\lambda hr^2\int\sin^2\theta\,d\theta \\
&=\lambda hr^2\int\frac{1-\cos2\theta}{2}\,d\theta \\
&=\frac12\lambda hr^2\int1-\cos2\theta\,d\theta \\
&=\frac12\lambda hr^2\int1\,d\theta-\int\cos2\theta\,d\theta \\
&=\frac12\lambda hr^2\left[\theta-\frac12\sin2\theta\right] \\
\\
\text{substitute back } \theta=\arcsin\left(\frac1r y\right)\text{:} \\
&=\frac12\lambda hr^2\left[\arcsin\left(\frac1r y\right)-\frac12\sin\left(2\arcsin\left(\frac1r y\right)\right)\right] \\
\text{the boundaries of }y\text{ are }-r\text{ and }r\text{:} \\
&=\frac12\lambda hr^2\left[\arcsin\left(\frac{y}{r}\right)-\frac12\sin\left(2\arcsin\left(\frac{y}{r}\right)\right)\right]_{-r}^r \\
&=\frac12\lambda hr^2\left(\frac{\pi}{2}-0\right) – \left(-\frac{\pi}{2}-0\right) \\
&=\frac12\pi\lambda hr^2
\end{align}
Calculate Gauss' law for $S_2$ and $S_3$:
The surfaces of $S_2$ and $S_3$ are equal.
Since:
$\bullet$ the position of the center of the square part on the $xy$-plane is $(x,y)=(0,0)$, the direction of the electrostatic field at both surfaces is opposite: $(\lambda x \mathbf{i})$,
$\bullet$ and the unit normal vectors are in opposite direction,
the addition of the result of Gauss' law will not be equal to 0.
The surface of each of the surfaces is $\frac12 \pi r^2$.
The electric field in the $x$-direction is $\lambda x\mathbf{i}$.
$x$ for $S_2$ = $\frac12 h$.
$x$ for $S_3$ = $-\frac12 h$.
$\mathbf{\hat{n}}$ for $S_2$ = $\mathbf{i}$.
$\mathbf{\hat{n}}$ for $S_3$ = $-\mathbf{i}$.
Therefore for $S_2$:
\begin{align}
\mathbf{E}\cdot \mathbf{\hat{n}} \times \text{surface area}&=\lambda x\mathbf{i} \cdot \mathbf{i} \times \frac12 \pi r^2 \\
&=\lambda \frac12 h\mathbf{i} \cdot \mathbf{i} \times \frac12 \pi r^2 \\
&=\frac14 \pi\lambda hr^2 \\
\end{align}
And for $S_3$:
\begin{align}
\mathbf{E}\cdot \mathbf{\hat{n}} \times \text{surface area}&=\lambda x\mathbf{i} \cdot -\mathbf{i} \times \frac12 \pi r^2 \\
&=\lambda (-\frac12 h)\mathbf{i} \cdot -\mathbf{i} \times \frac12 \pi r^2 \\
&=\frac14 \pi\lambda hr^2 \\
\end{align}
Calculate Gauss' law for $S_4$:
Since $S_4$ lies in the $xy$-plane, the electrostatic field $\mathbf{E}=\lambda(x\mathbf{i}+y\mathbf{j})\;\text{,}$ lies parallel to the surface, thus the result of Gauss' law is $0$.
The net result is:
\begin{align}
\iint_S\mathbf E\cdot \mathbf{\hat{n}}\,dS&=\frac{q}{\epsilon_0} \\
&=\frac12 \pi\lambda hr^2 + \frac14 \pi\lambda hr^2 + \frac14 \pi\lambda hr^2 +0 \\
&=\pi\lambda hr^2 \\
\end{align}
The total charge $q$ enclosed in the half-cylinder is thus:
$$q=\pi\lambda hr^2 \epsilon_0$$.
This solves the problem I was having.
Best Answer
The lower and upper boundaries for $y$ in $S_1$ should be $-r$ and $r$ instead of $0$ and $r$. I was thinking in polar coordinates instead of in Cartesian coordinates...
I'll edit the original post accordingly