I need to evaluate $$ \frac{\partial^2W}{\partial F^2}(I)(F,F)=\sum_{i,j,k,l=1}^3\frac{\partial^2W}{\partial F_{ij}F_{kl}}(I)F_{ij}F_{kl} $$ for $W=\det\sqrt{F^TF}$. Here $F\in\mathbb{R}^{3\times3}$.
From First derivative we know that $$ \frac{\partial W}{\partial F}=WF^{-T}. $$ To evaluate the sum above, I write $$ \frac{\partial^2W}{\partial F^2}=W\left[\left(F^{-T}\right)^2+\frac{\partial F^{-T}}{\partial F}\right]․ $$
Here I got stack.
Denote the gradient of $W$ as $$\eqalign{ G &= \frac{\partial W}{\partial F}=F^{-T}W \cr \frac{1}{W}\,G &= F^{-T} \cr }$$ and find its differential and gradient $$\eqalign{ dG &= F^{-T}\,dW + W\,dF^{-T} \cr &= F^{-T}(G:dF) - WF^{-T}\,dF^{T}\,F^{-T} \cr &= F^{-T}(G:dF) - WF^{-T}\,{\mathbb E}\,F^{-1}:dF^{T} \cr &= \frac{1}{W}\,\Big(G\star G - (G\,{\mathbb E}\,G^T):{\mathbb B}\Big):dF \cr\cr \frac{\partial G}{\partial F} &= \frac{1}{W}\,\Big(G\star G - (G\,{\mathbb E}\,G^T):{\mathbb B}\Big) \cr }$$ where $(\star)$ denotes the dyadic product, $(:)$ denotes the Frobenius product, and $({\mathbb E, B})$ are fourth-order isotropic tensors with components $$\eqalign{ {\mathbb E}_{ijkl} &= \delta_{ik}\,\delta_{jl} \cr {\mathbb B}_{ijkl} &= \delta_{il}\,\delta_{jk} \cr\cr }$$ So, as expected, the second derivative of $W$ is a fourth-order tensor $$\eqalign{ \frac{\partial^2 W}{\partial F^2} &= \frac{\partial G}{\partial F} \cr\cr }$$
For the sake of completeness, I'll mention that there is one more isotropic tensor ${\mathbb A}=I\star I$ with components $$\eqalign{ {\mathbb A}_{ijkl} &= \delta_{ij}\,\delta_{kl} \cr }$$ We didn't need it for this problem, but it does illustrate how the dyadic product works.
Anyway, the isotropic tensors have properties which are useful for working with matrix expressions $$\eqalign{ P:{\mathbb E} &= {\mathbb E}:P = P \cr P:{\mathbb B} &= {\mathbb B}:P = P^T \cr P:{\mathbb A} &= {\mathbb A}:P = I(I:P) = P\,\operatorname{tr}(P) \cr P\cdot Q\cdot R &= (P\cdot {\mathbb E}\cdot R^T):Q \cr P\cdot Q^T\cdot R &= (P\cdot {\mathbb E}\cdot R^T):{\mathbb B}:Q \cr \cr }$$ If you are uncomfortable with tensors, you can apply vectorization $$\eqalign{ g &= \operatorname{vec}(G) \cr f &= \operatorname{vec}(F) \cr }$$ and then calculate $$\frac{\partial g}{\partial f}$$ which will be an ordinary matrix.
Update
Now I understand the whole question. You wish to evaluate the above Hessian at $F=I$ and then fully contract it with the tensor $(F\star F)$.
At $F=I$, the gradient reduces to $G=WI$ and Hessian reduces to $$\eqalign{ {\mathbb H} &= \frac{\partial G}{\partial F} = W\,\Big(I\star I - {\mathbb E}:{\mathbb B}\Big) \cr &= W\,({\mathbb A}-{\mathbb B}) \cr\cr }$$ Now it can be contracted to a scalar value $$\eqalign{ F:{\mathbb H}:F &= W\,\Big(F:({\mathbb A}-{\mathbb B}):F\Big) \cr &= W\,\Big(\operatorname{tr}(F)^2 - \operatorname{tr}(F^2)\Big) \cr }$$
