Let $\mathbb{T}^n:(\mathbb{R}/\mathbb{Z})^n$ be $n$-torus.
For $p \in (1,\infty)$ and $u \in W^{2,p}(\mathbb{T}^n)$, I am aware that there exists a constant $C(n,p)>0$ such that \begin{equation} \lVert D^2 u \rVert_{L^p(\mathbb{T}^n)} \leq C(n,p)\Bigl( \lVert \Delta u \rVert_{L^p(\mathbb{T}^n)} +\lVert u \rVert_{L^p(\mathbb{T}^n)} \Bigr) \end{equation} where $(D^2 u)_{ij}=\partial_i \partial_j u \in L^p(\mathbb{T}^n)$ for $i,j=1,\cdots,n$.
This seems to have the name Calderón-Zigmund duality.
But, I cannot really find a reference with explicit formula of the above kind..
Could anyone please suggest me one?
Results of this type (controlling the $L^p$ norm of $\mathrm{D}^2u$ by $\Delta u$) are generally referred to as Calderón-Zygmund estimates. This terminology comes from how in the full-space, where we can solve the problem $\Delta u = f$ explicitly as a convolution $ u = \Gamma \ast f$ where $\Gamma$ is the Newtonian potential. Then formally $\mathrm{D}^2u = \mathrm{D}^2\Gamma \ast f$, and the work of Calderón and Zygmund allows one to obtain $L^p$ estimates for convolution operators of the type $f \mapsto K \ast f$.
Note this is not a duality result (though duality arguments are often present at some stages of the proof). There is also no formula, but rather it will use technical arguments from harmonic analysis in one way or another.
Results of this type are 'well-known' within the community, but detailed proofs adapted to specific situations are often difficult to track down. Typically one establishes these as follows:
Establish the estimate in the full space, for instance to show $$ \lVert \mathrm{D}^2\Gamma \ast f \rVert_{L^p(\mathbb R^n)} \leq C \lVert f \rVert_{L^p(\mathbb R^n)}. $$
For general domains $\Omega$, establish local estimates by cutting off $u$ away from the region of interest. That is choose a cutoff $\eta \in C^{\infty}_c(\Omega)$ and note that $$ \Delta(\eta u) = \eta \underbrace{\Delta u}_{=f} + 2\nabla u \cdot \nabla \eta + u \Delta \eta. $$ One argues the estimate in Step 1 to $\eta u$; here one has terms involving $u$ and $\mathrm{D}u$ on the right-hand side, and one needs an interpolation argument to estimate the intermediate derivative (e.g. $\lVert \mathrm{D}u\rVert_{L^p} \lesssim \lVert u \rVert_{L^p} + \lVert \mathrm{D}^2u\rVert_{L^p}$). This also works for the torus, viewing $u$ as a function on $\mathbb R^n$ via a periodic extension.
Patch these local estimates to infer a global estimate. For domains $\Omega$ one would also need local estimates near the boundary, but for the torus these do not appear.
The above is the classical way to prove these estimates, but there are also other approaches based on maximal functions or Fourier analysis (the latter may be more natural for the torus). All approaches involve technical machinery in some form or another to establish $L^p$ estimates.
Details can be found in the following references. Note that both are for domains; for elliptic regularity it's common that many results are not explicitly proven for the manifold case, but one can derive it from the corresponding flat case (as I sketched above).
Gilbarg, David; Trudinger, Neil S., Elliptic partial differential equations of second order, Grundlehren der mathematischen Wissenschaften. 224. Berlin-Heidelberg-New York: Springer-Verlag. X, 401 p. (1977). ZBL0361.35003.`
Giusti, Enrico, Direct methods in the calculus of variations, Singapore: World Scientific. vii, 403 p. (2003). ZBL1028.49001.
The latter reference only covers the range $2 \leq p < \infty$, but a suitable duality argument can extend it to the full range. Many texts on elliptic PDE will also cover this topic, though different texts may cover different settings and present different proofs.