Moment of Inertia, why $r^2$and not $r$? So my engineering mechanics book includes a brief discussion on area moments of inertia.
Unfortunately, the ensuing chapter is predominately computational in nature.  I don't have a thorough grasp of where the equations for this come from. The equation in question is expressed as:
$$I_x = \int_A y^2dA \space \mathrm{ and } \space I_y = \int_A x^2dA$$
Because of the simple equation for torque $\tau=Fd$, it is easy to imagine that the angular momentum of an object being rotated about an axis is greater the further it's mass is from the axis.
Simply divide any continous object up into a set of differential elements and violá, you have an integral that can be vaguely matched to the angular momentum.
The problem, however, is I can't possibly derive this.  In fact, it makes more sense to me to integrate over simply $x$ or $y$ rather than $x^2$ or $y^2$.
Due to a logistical nightmare, I have no access to a physics book to look this up.  I would be grateful to anyone who could explain (ie, explain the motivation of) this to me!
 A: I think you are confussed between mass moment of inertia and area moment of inertia.
The first is an equivalent of mass in angular direction and is defined as $\int_V{r^2\rho dV}$. An angular equivalent of $F=ma$ is:
$$\tau=I\alpha$$ where $\tau$ is torque (angular equivalent of force, with units $[Nm]$), $I$ is mass moment of inertia (angular equivalent of mass, with units $[kgm^2]$) and $\alpha$ is angular acceleration (angular equivalent of linear acceleration, with units $[rad/s^2]$).
The second is an indication how well a beam can withstand a torque (how much it will bend due to it). There are many equations which can relate to this and depend on how a (cantilever) beam is constrained. But this has the definition of $\int_A{rdA}$ and has the unit $[m^4]$.
Edit:
After seeing a comment it is clear that you mean the second. And I will try to explain it a bit better. Consider a section of a beam on which a constant bending moment is applied:

This section will be compressed on one end and stretched on the other end and somewhere in the middle (depending on the shape of the section) there will be no stresses/deformation. In equation form it looks like this:
$$
\sigma(y)=\frac{My}{I},
$$
where $\sigma$ is the stress on the beam at a distance $y$ from the beams neutral axis, $M$ is the applied bending moment.
A way to derive this is that the integral of the stresses times a small cross-section area times the the distance from the neutral axis should give the bending moment:
$$
M=\int_A{\sigma(y)ydA}
$$
When substituting the equation for $\sigma$ into this integral and divide both sides by $\frac{M}{I}$ you get:
$$
I=\int_A{y^2dA}
$$
I am not sure if this makes it more clear, since this might seem as an circular reasoning, since I assume a certain equation for $\sigma$. But I hope that you can see that the stresses will be proportional to the bending moment when assuming small and elastic deformations ($\sigma\propto M$). And that the stresses will also be linearly depend on the distance away from the neutral axis, since $\sigma dA=dF$ and $Fr=M$, so $\sigma rdA=dM$.
A: Here's a motivation for where the inertia tensor $I=(I_{ij})$ (and by extension moments of inertia) comes from.  It's a quantity that's analogous to mass for rotational motion in the sense that the kinetic energy of a rotating object is essentially proportional to the inertia tensor times the square of the body's angular velocity.  More precisely
\begin{align}
  T(t) = \frac{1}{2} \boldsymbol \omega(t)^tI(t)\boldsymbol \omega(t). \tag{1}
\end{align}
where $\boldsymbol \omega(t)$ is the instantaneous angular velocity of the body.  Compare this, for example, to the expression for the kinetic energy of a particle of mass $m$ moving with speed $v$;
\begin{align}
  T = \frac{1}{2}mv^2.
\end{align} 
To prove expression $(1)$, start with a rigid body consisting of points $\mathbf x_i$ undergoing pure rotation.  There exists a time-dependent rotation $R(t)$ that generates the motion of all points in the rigid body;
\begin{align}
  \mathbf x_i(t) = R(t)\mathbf x_i(0) \tag{2}
\end{align}
The kinetic energy of the body is the sum of the kinetic energies of the individual particles;
\begin{align}
  T(t) 
&= \frac{1}{2}\sum_i m_i \dot{\mathbf x}_i(t)\cdot\dot{\mathbf x}_i(t) \\
&= \frac{1}{2}\sum_i m_i \big(\dot R(t)\mathbf x_i(0)\big)\cdot\big(\dot R(t)\mathbf x_i(0)\big) \\
&=\frac{1}{2}\sum_i m_i \big(\dot R(t)R(t)^t\mathbf x_i(t)\big)\cdot\big(\dot R(t)R(t)^t\mathbf x_i(t)\big) \\
\end{align}
where in the last equality I used the fact that $R^tR = I$ for rotations so that eq. $(2)$ gives $\mathbf x_i(0) = R(t)^t \mathbf x_i(t)$.  Now, we note that
\begin{align}
  \dot R(t)R(t)^t\mathbf x_i(t) = \boldsymbol\omega(t)\times\mathbf x_i(t)
\end{align}
where $\boldsymbol\omega$ is the angular velocity vector of the body.  See the following for a detailed derivation of this fact:
Angular Velocity expressed via Euler Angles
So putting this together, we have
\begin{align}
  T(t) 
&= \frac{1}{2}\sum_im_i\big(\boldsymbol\omega(t)\times \mathbf x_i(t)\big)\cdot \big(\boldsymbol\omega(t)\times \mathbf x_i(t)\big) \\
&= \frac{1}{2}\sum_im_i \sum_{j,k}\omega_j\big(\mathbf x_i^2\delta_{jk} - (x_i)_j(x_i)_k\big)\omega_k \\
&= \frac{1}{2} \sum_{j,k}\omega_j\left[\sum_im_i\big(\mathbf x_i^2\delta_{jk} - (x_i)_j(x_i)_k\big)\right]\omega_k \\
\end{align}
Now, if we simply note that the inertia tensor is defined as the quantity whose components $I_{jk}$ are in the big brackets, then we have the desired formula.
Note, in particular, that when $j=k$, namely when we consider only the diagonal components of the inertia tensor, then we obtain the $j$th moment of inertia
\begin{align}
  I_{jj} = \sum_i m_i\big(\mathbf x_i^2 - (x_i)_j^2\big)
\end{align}
so, for example, the $x$ moment is
\begin{align}
  I_{xx} = \sum_i m_i(y_i^2+z_i^2)
\end{align}
and if the object is in the $x$-$y$ plane, then $z=0$ and we get
\begin{align}
  I_{xx} = \sum_i m_i y_i^2
\end{align}
and if the body is continuous, then sums get replaced with the appropriate integrals;
\begin{align}
  m_i\to dm, \qquad I_{xx}\to \int y^2 dm
\end{align}
A: The moment of inertia relates to resisting rotational acceleration, which involves motion. With respect to a rotational influence twisting at the axis, the radius affects the amount of this resistance two separate (although intertwined) ways. Firstly, if the radius is doubled then the mechanical advantage of the lever arm (radius) is reduced by a factor of two, it will require twice the torque at the axis to produce the same force out at the end of the lever arm. Secondly, doubling the radius doubles the length of the arc sweep when motion occurs. For any arbitrary bit of rotation at the axis, doubling the radius will double the distance along the arc swept by the end of the lever arm (radius), necessitating the tip of the lever arm to travel at twice the velocity. 
So an increased radius participates in resisting rotational acceleration two ways, accounting for the appearance of the $r^2$ rather than simply $r$.
Put another way, when you double the radius, for any given twist applied at the axis, you are trying to accelerate the mass twice as quickly with only half as much leverage.
A: The angular momentum of a single particle with mass $m$ in motion about an axis, with angular speed $\omega$, a distance $r$ from the axis, is
$L = r (m v) = m r^2 \omega$.
When we consider an extended body, the sum up the contribution ($m r^2$) from each particle in motion inside the body, and this is the moment of inertia.
More generally,
$$\begin{align}
\mathbf{L} &= \sum_p m_p \mathbf{r}_p \times \mathbf{v}_p\\
& = \sum_p m_p \mathbf{r}_i \times (\boldsymbol{\omega} \times \mathbf{r}_i)\\
& = \sum_p m_p \left( \boldsymbol{\omega} \, (\mathbf{r}_p \cdot \mathbf{r}_p)  - \mathbf{r}_p (\mathbf{r}_p \cdot \boldsymbol{\omega})\right)\\
& = \sum_p m_p \left(  (\mathbf{r}_p \cdot \mathbf{r}_p)\boldsymbol{\omega}  - \mathbf{r}_p (\mathbf{r}_p \cdot \boldsymbol{\omega})\right)\\
L_j &= \sum_p m_p \sum_k \left(  r_p^2 \delta_{jk} - (r_p)_j (r_p)_k \right) \omega_k \\
&=  \sum_k I_{jk} \omega_k \\
\end{align}$$
The term $\sum_p m_p \left(  r_p^2 \delta_{jk} - (r_p)_j (r_p)_k \right) $ is the moment of inertia tensor $I_{jk}$
Note that $I_{xx} = \sum_p m_p (r^2 - x^2) = \sum_p m_p (y^2 + z^2)$. In 2 dimensions, you just ignore the $z^2$ part, so this says that $I_{xx} = \sum_p m_p y^2$. For a continuous body of uniform density, we would get $I_{xx} = \rho \int\!dA\, y^2$.
Now suppose this body is rotating about the z axis (even though it's 2 dimensional in your case, and $z=0$, we put in a z axis to give something to rotate about; Really we are rotating in the x-y plane.) Then the angular velocity vector is $\boldsymbol{\omega} = \omega \,\hat{\mathbf{z}}$. We would write
$$
\begin{align}
L_z &= L_3\\
& = \sum_k I_{3k} \omega_k\\
& = I_{33} \omega\\
& = \rho \int\!dA\, (x^2 + y^2) \omega\\
\end{align}
$$
A: The other answers here are great, but if they don't work for you, here is one other thing you might consider...
In one dimension, we find the center of mass with an average over the masses
$$ \bar x=\frac {\sum m_i x_i}{\sum m_i} $$
We find the gyro-radius (one dimension still) with the root mean square over the masses
$$ x_{gyro}=\sqrt{\frac {\sum m_i x_i^2}{\sum m_i}} $$
The root mean square tells us how spread out something is.  And the more spread out the mass is, the harder it is to "rev up", to get spinning and once it is spinning, the harder it is to stop.  The gyro-radius gives us the distance from the center all the mass could be relocated to and have the same rotational inertia.
So we do need the square in the integral, since this gives us how spread out the mass is rather than the first power which gives us the center of mass.
