1. Why is force linear in the velocity?
On a conceptual level, the idea behind stokes flow is that it is a very overdamped system because the viscosity is very high. When I say "overdamped" I am trying to make an analogy with a damped harmonic oscillator. You gave the equation
$$\ddot{x}=-k^2 x,$$
but a damped harmonic oscillator has a drag term:
$$\ddot{x}=-k^2 x - \gamma \dot{x}.$$
Now the limit of stokes flow is like the limit and $\gamma \gg k$. In this case there is a lot of drag and so everything goes slowly. Since the $\ddot{x}$ term involves to powers of time in the denominator, and the time scale has become very large, the $\ddot{x}$ term ought to be very small. In fact we would expect $\dot{x}$ to be small as well. The only reason that term is significant is because $\dot{x}$ is multiplied by the very large $\gamma$.
In fact we can look the solution and see exactly how things play out. The overdamped solution the form
$$x(t) = A e^{-r_-t} + Be^{-r_+ t},$$
where the decay rates $r_\pm$ are given by
$$ r_\pm = \frac{1}{2}\left(\gamma \pm \sqrt{\gamma^2-4 k^2}\right).$$
In the limit of large $\gamma$, we get to lowest order
$$ r_+ = \gamma$$
and
$$ r_- = k^2/\gamma.$$
Therefore $r_+ \gg r_-$, and so the $r_+$ term decays very quickly and we can ignore it. We are left with
$$x(t)=A e^{-\left( k^2/\gamma \right) t}.$$
With this solution, we can see the $\ddot{x}$ term goes like $\dfrac{k^4}{\gamma^2}$, the $\dot{x}$ term goes like $k^2$, and the $x$ term also goes like $k^2$. Therefore the intertial term goes to zero, and the external force acting on the object (in this case the spring force) is equal to the drag force.
Now in the navier stokes equation, the inertial term looks like
$$\newcommand{\u}{\mathbf{u}}\rho\left(\dot{\u} + \left(\u \cdot \nabla\right) \u\right).$$
This term contains both a non-linearity, and a $\dot{\u}$, which provides the system with memory, so that the current state can depend on the past. But as we saw, in the overdamped case this inertial term is neglected, and the remaining terms are linear and have no time derivatives. Looking at the equation with the remaining terms:
$$\newcommand{\z}{\mathbf{0}} \z = \nabla p - \mu \nabla^2 \u,$$
we see that the present value of $p$ is linearly related to the present value of $\u$. Since the force is basically the integral of the pressure over the surface of an object, we expect the force to be linearly related to the velocity.
2. Why is the propulsion matrix symmetric
The symmetry of the propulsion matrix basically follows from the symmetry of the laplacian. I will go into detail in the following two steps. First I will show that the green's function for the stokes equation is symmetric, then I will show why this symmetry gives you the symmetry of the propulsion matrix.
Symmetry of the green's function
Before I get into why the green's function for the stokes equation is symmetric, let's consider the simpler case of electrostatics. Here we have that the charge density $\rho$ is related to the electric potential by $\rho = \nabla^2 \phi$. Now notice that $\nabla^2$ is a symmetric operator in the sense that $$\int \newcommand{\r}{\mathbf{r}}\phi_1(\r) \nabla^2 \phi_2(\r) d \r = \int \newcommand{\r}{\mathbf{r}}\phi_2(\r) \nabla^2 \phi_1(\r) d \r,$$
as can easily be seen by integrating by parts. Now since the green's function is the inverse of the operation $\nabla^2$, and the inverse of a symmetric operator is symmetric, we expect the green's function should be symmetric. In fact, this can be directly verified since we know the inverse of $\nabla^2$ is the operator $$\phi(\r) \mapsto \int d\r_0 \frac{1}{4 \pi |\r-\r_0|} \phi(\r_0).$$
We can easily see that $$ \int \phi_1(\r) \int d\r_0 \frac{1}{4 \pi|\r-\r_0|} \phi_2(\r_0) d\r = \int \phi_2(\r) \int d\r_0 \frac{1}{4 \pi|\r-\r_0|} \phi_1(\r_0) d\r.$$
In particular, plugging in $\phi_1(\r)=\delta(\r-\r_s)$ and $\phi_2(\r)=\delta(\r-\r_o)$, we find that the potential at $\r_o$ due to a source chrage at $\r_o$ is $\frac{1}{4 \pi |\r_s -\r_o|}$, which is the same as the potential at $\r_s$ from a charge at $\r_o$.
Now switching back to the stokes equation, if we add a term for an external force density $\newcommand{\f}{\mathbf{f}}\f$, we get the equation $\f=\nabla p - \mu \nabla^2 \u$. This equation looks a lot like the relation between charge density and potential from electrostatics. The difference is that there is now an extra field $p$. These extra degrees of freedom would make the problem underdetermined, but in stokes flow, we have an extra incompressibility constraint: $\nabla \cdot \u = 0$. Although this differential operator is not so simple in this case, it is still symmetric, and so the green's function is still symmetric. In fact, the green's function is given by
$$\f(\r) \mapsto \int \frac{1}{8 \pi \mu |\mathbf{\Delta}|} \left( \mathbb{I} + \hat{\Delta} \otimes \hat{\Delta} \right) \f(\r_0) dr_0, $$
where $\mathbf{\Delta} = \r - \r_0$. We see that this green's function is basically the same as the greens function for the laplacian execept for two differences. One is that since we started with a vector differential equation, the green's function is a tensor. The second difference is the addition of the $\hat{\Delta} \otimes \hat{\Delta}$ piece, which enforces incompressibility.
Symmetric green's function implies symmetric propulsion matrix
Now that we have found the green's function, we will next find the propulsion matrix for an object. The object is characterized by its two-dimensional surface $S$.
We want to consider what happens if the object is subjected to an external force $\newcommand{\F}{\mathbf{F}}\F$ and torque $\newcommand{\t}{\boldsymbol{\tau}}\t$. To see the effect of the force and torque, it is necessary only to find the surface forces they create. (It is not necessary to find the internal stress distribution in the object.) This is because momentum is conserved locally, so we must have that for each point on the surface, the hydrodynamic momentum flux balances the external momentum flux. The hydrodynamic momentum flux is $$p \hat{n} -\mu \left(\hat{n} \cdot \nabla\right) \u,$$ where $\hat{n}$ is the normal to the surface $S$. This must be balanced by the surface force density from the stress distrubtion of the object, which is given by $\sigma \cdot \hat{n}$, where $\sigma$ is the stress in the object. For convenience of notation, I will refer to the external surface force density simply as $\f$ instead of $\sigma \cdot \hat{n}$. So the net effect of what I have said this paragraph is that any force $\F$ and torque $\t$ can be represented by a surface force density $\f$ defined on $S$, which will serve as a boundary condition for the flow field $\u$.
Given $\f$, we already know that the solution for $\u$ is
$$ \u(\r) = \int_S \frac{1}{8 \pi \mu |\mathbf{\Delta}|} \left( \mathbb{I} + \hat{\Delta} \otimes \hat{\Delta} \right) \f(\r_0) d \r_0, $$
where again $\mathbf{\Delta}=\mathbf{\r}-\mathbf{\r_0}$.
Now the propulsion matrix is not concerned with the velocity $\u$ everywhere, but only $\u$ on the surface $S$. But we may simply consider the restriction of $\u$ to $S$. If we do that, then we have a linear relationship
$$\newcommand{\M}{\mathbf{M}}\u = \M \f$$
between the two vector fields $\u$ and $\f$, where $\M$ is the symmetric linear operator defined by
$$\f(\r) \mapsto \int _S \frac{1}{8 \pi \mu |\mathbf{\Delta}|} \left( \mathbb{I} + \hat{\Delta} \otimes \hat{\Delta} \right) \f(\r_0) dr_0. $$
Now the operator $\M$ takes an $\f$ and gives a $\u$, but the propulsion matrix takes a $\u$ and gives an $\f$. Therefore we must invert $\M$. Now I won't get into the question of whether $\M$ is invertible, because it will be sufficient to take the pseudo-inverse $\M^+$. This pseudoinverse has the property that if $\u_0 = \M \f$ has a solution for $\f$, then the solution is given by $\f = \M^+ \u_0$. This will be an important property. Also, since $\M$ is symmetric, the pseudo-inverse of $\M$ must also be symmetric.
Now we are only interested in rigid motions of a body with velocity $\newcommand{\v}{\mathbf{v}}\v$ and angular velocity $\newcommand{\w}{\boldsymbol{\omega}}\w$. We now define a linear operator $\newcommand{\U}{\mathbf{U}}\U$ defined by the map
$$\U (\v,\w) = \v + \w \times \r,$$
where the right hand side is a vector field on $S$.
Then we have $\f = \M^+ \U \left( \v , \w\right)$. Here we are relying on the property that if a solution exists, then the pseudoinverse will give the solution, for we believe that there ought to be some surface force distribution $\f$ that gives rise to the rigid motion velocity field given by $\U \left( \v , \w\right)$. Now we are interested in calculating the force and torque from $\f$. The force is simply given by
$$\F=\int_S \f d\r$$
and
$$\t = \int_S \r \times \f d\r.$$
This can be written as $\left( \F, \t\right) = \newcommand{\V}{\mathbf{V}}\V \f$.
Putting everything together, we have
$$\left( \F, \t\right) = \V \M^+ \U \left( \v , \w\right).$$
Therefore the quantity $\V \M^+ \U$ is the propulsion matrix. To show the propulsion matrix is symmetric all we have to show is that $\V = \U^\dagger$. But this is the same as showing $\langle \V \f, (\v,\w)\rangle = \langle \f, \U (\v,\w)\rangle$. But the left hand side is
$$ \v \cdot \int_S \f d \r + \w \cdot \int_S \r \times \f d \r $$
and the right hand side is
$$ \int \f \cdot \left( \v + \w \times \r \right) d \r.$$
Expanding the above line and rearranging the triple product in the second term, we see the right hand side is indeed equal to the left hand side, and so $\V$ is indeed $\U^\dagger$. And so the propulsion matrix is indeed symmetric.
