Sums of this sort are well known to the experts - but since none of them have answered so far, let me try. Denoting the fractional part of a real $x$ by ${x}$, \{x\}$, you can write your sum as
$$ S = n^2 \sum_{d=1}^n \frac1{d^2} - 2n \sum_{d=1}^n \frac1d\left\{\frac nd\right\} + \sum_{d=1}^n \left\{\frac nd\right\}^2 = \frac{\pi^2}{6}\, n^2 + O(n\log n). $$
If you need more precision, you have to find the main term of the sum $\sum_d d^{-1}\{d^{-1}n\}$. The standard technique here, I believe, would be to use the Fourier expansion of the fractional part function, but you'd better contact experts for details, to avoid re-inventing the wheel.
Here is a different kind of answer, depending on what you are after. Your sum counts the number of triples $d,x,y\in[1,n]$ with $xd,yd\le n$. Since there are $\sum_{k=1}^n \tau(k)$ such triples with $x=y$, splitting the sum into two part according to whether $x\ge y$ or $y\ge x$, we can write it as
$$ S = 2 \sum_{dx\le n} x - \sum_{k=1}^n \tau(k). $$
Letting $k=dx$, we get
$$ S = 2 \sum_{k=1}^n \sigma(k) - \sum_{k=1}^n \tau(k), $$
where $\sigma$ is the sum-of-divisors function. This gives you ``an identity containing arithmetic functions'', as you requested.
