Average order of the reciprocal of sum of divisors 
Question: What is the order of $$S(x)=\sum_{n\le x}\frac{1}{\sigma(n)},$$ where $\sigma(n)=\sum_{d|n}d$?

Numerical experiment suggests $S(x)\sim C\ln x$ for some constant $C$. How can we find $C$?
Some Quick Results
$$S(x)<\sum_{n\le x}\frac1n\sim\ln x$$
since $\frac{1}{\sigma(n)}<\frac{1}n$.
Assuming RH, we know $$\frac{1}{\sigma(n)}>\frac{1}{e^\gamma n\ln\ln n}$$ for sufficient large $n$, which means $S(x)$ is bounded asymptotically by $\int^x\frac{dx}{x\ln\ln x}\sim\frac{\ln x}{\ln\ln x}$
 A: A standard method in the mean value of multiplicative arithmetic function is the use of Dirichlet series. 
Let
\begin{align}
G(s)&=\prod_p \left(1+\frac{p-1}{p^s(p^2-1)}+\frac{p-1}{p^{2s}(p^3-1)}+\cdots\right)\left(1-\frac1{p^{s+1}}\right)\\&=\prod_p\left(1+\sum_{n=1}^{\infty} \frac{p-1}{p^{ns}(p^{n+1}-1)}\right)\left(1-\frac1{p^{s+1}}\right)\\
&=\prod_p\left(1-\sum_{n=1}^{\infty} \frac{(p-1)^2}{p^{ns}(p^{n+1}-p)(p^{n+1}-1)}\right)=\sum_{n=1}^{\infty} \frac{g(n)}{n^s}
\end{align}
which converges absolutely for $\Re(s)>-1$.
Then we have 
$$
\frac1{\sigma(n)}=\sum_{d|n} g(d) \frac nd = \sum_{dk=n} g(d)\frac 1k.
$$
Then
\begin{align}
\sum_{n\leq x}\frac1{\sigma(n)}&=\sum_{d\leq x}g(d)\sum_{k\leq \frac xd} \frac1k\\&=\sum_{d\leq x}g(d)\left( \log\frac xd + \gamma+O\left(\frac dx\right)\right)\\&=\sum_{d\leq x} g(d)\left(\log x - \log d + \gamma\right) + O(x^{-1+\epsilon})\\
&=G(\log x + \gamma)+H+O(x^{-1+\epsilon}),
\end{align}
where 
$$
G=G(0), \ H=G'(0).
$$
To compute $G(0)$, we put $s=0$ in the Euler product, 
$$
G(0)=\prod_p \left(1+\sum_{n=1}^{\infty} \frac{p-1}{p^{n+1}-1}\right)\left(1-\frac1p\right).
$$
To compute $G'(0)$, we use the logarithmic differentiation of the Euler product,
$$
\frac{G'(0)}{G(0)}=\sum_p \left(\frac{-\sum
\limits_{n=1}^{\infty} \frac{n(p-1)}{p^{n+1}-1}\log p}{1+\sum\limits_{n=1}^{\infty} \frac{p-1}{p^{n+1}-1}} +\frac{\log p}{p-1}\right).
$$
Running the SAGE code, 
s=1
for p in primes(1,10000):
    b=0
    for n in range(2,40):
        b+=1/(p^n-1)
    s=s*(1+(p-1)*b)*(1-1/p)
print numerical_approx(s)

gives an approximation of $G$, 0.672744912303011
t=0 
for p in primes(1,10000):
    b=0
    c=0
    for n in range(2,40):
        b+=1/(p^n-1)
        c+=1/(p^n-1)*(n-1)
    t+=(-(p-1)*c)/(1+(p-1)*b)*log(p)+(log(p)/(p-1))
print numerical_approx(t)

gives an approximation of $G'(0)/G(0)$,  0.507239459457213 
Then $G\gamma + H$ is approximately 
.672744912303011(.577215664901533+.507239459457213)
=0.72956166753
A: Not a real answer, just for fun:
I played a little with Mathematica. Here is how I defined $S(x)$ in it:
sfun[1] = 1

sfun[n_] := sfun[n] = sfun[n - 1] + N[1/DivisorSum[n, # &], 100] 

And here is the plot of $S(x)$ for $x\in[2,400000]$:
 
If you create a logarithmic plot with $\ln x$ on the horizontal axis, the plot becomes almost a straight line:

I tried to interpolate it and got the following expression:
$$S(x)\approx 0.731115 + 0.672623 \ln x\tag{1}$$ 
...which gives some fairly accurate predictions. For example, the exact value for $S(200000)$ is 8.94111 while my approximation predicts 8.94119.
Then I tried some extrapolation. The accurate value for $S(800000)$ is 9.87372 while (1) predicts 9.87365. So (1) seems to be fairly accurate.
