You're not getting your facts right at all.
How do we know from this $\langle W \rangle = \int_{-\infty}^{\infty} \bar{\Psi}\left(-\frac{\hbar^2}{2m} \frac{d^2}{dx^2} + W_p \right) \Psi dx$ or this $\hat{H} = -\frac{\hbar^2}{2m}\frac{d^2}{dx^2} + W_p$ that we have an eigenfunctiuion and eigenvalue.
Answer: we don't.
All I know about operator $\bar{H}$ so far is this equation where $\langle W \rangle$ is an energy expected value:
\begin{align}
\langle W \rangle = \int_{-\infty}^{\infty} \bar{\Psi}\left(-\frac{\hbar^2}{2m} \frac{d^2}{dx^2} + W_p \right) \Psi dx
\end{align}
No, you don't.
Here's the mathematical side of what an eigenfunction and eigenvalue is:
Given a linear transformation $T : V \to V$, where $V$ is an infinite dimensional Hilbert or Banach space, then a scalar $\lambda$ is an eigenvalue if and only if there is some non-zero vector $v$ such that $T(v) = \lambda v$.
Here's the physics side (i.e. QM):
We postulate that the state of a system is described by some abstract vector (called a ket) $|\Psi\rangle$ that belongs to some abstract Hilbert space $\mathcal{H}$.
Next we postulate that this state evolves in time by some Hermitian operator $H$, which we call the Hamiltonian, via the Schrodinger equation. What is $H$? you guess and compare to experimental results (that's what physics is anyway).
Next we postulate for any measurable quantity, there exists some Hermitian operator $O$, and we further postulate that the average of many measurements of $O$ is given by $ \langle O \rangle = \langle \Psi | O | \Psi \rangle$.
Connection to wavefunctions: we pick the Hilbert space $L^2(\mathbb{R}^3)$ to work in, so $\Psi(x) = \langle x | \Psi \rangle$, and $\langle O \rangle = \int_{-\infty}^{\infty} \Psi^*(x) O(x) \Psi(x) dx$.
Ok, that's the end. The form of $H$ doesn't follow from the energy expected value.
Wait! I haven't even talked about eigenvalues and eigenfunctions. This is a useless post!
Answer: well you don't have to. But it is useful to find the eigenvalues and eigenfunctions of $H$, because the eigenfunctions of $H$ form a basis of the Hilbert space, and certain expressions become diagonal/more easily manipulated when we do whatever calculations we want to do.
So to find the eigenvalues of $H$, we simply solve the eigenvalue equation as stated above:
Solve
\begin{align}
H | \Psi_n \rangle = E_n | \Psi_n \rangle.
\end{align}
This is in the form $T(v) = \lambda v$.
So as Alfred Centauri says, we simply want to find the eigenfunctions of $H$. A more subtle question would be, how do we know they exist? The answer lies in spectral theory and Sturm-Liouville theory but nevermind for now, as physicists we assume they always exist.
So your additional question:
$\hat{a} \psi$ is an eigenfunction of operator$\hat{H}$ with
eigenvalue $(W-\hbar \omega)$.
Well.... that just follows straightaway. You said you already proved that $H a^\dagger \psi = (W - \hbar \omega) a^\dagger \psi$. So here $T$ = $H$, $a^\dagger \psi = v$, and $\lambda = (W - \hbar \omega)$. which is an eigenvalue equation $T(v) = \lambda v$. Thus, $a^\dagger \psi$ is an eigenfunction of $H$ with eigenvalue $(W-\hbar \omega)$.
There is nothing to "derive" for the kinetic energy operator. By definition, classical kinetic energy is $\frac{p^2}{2m}$, and so $\hat{E}_\text{kin} = \frac{\hat{p}^2}{2m}$ quantumly. It's not exactly clear why you think this doesn't make sense mathematically, but it does: In words, it says "apply the momentum operator twice, then divide the result by $2m$".
Note that $\langle p^2\rangle \neq \langle p\rangle^2$, the difference is precisely what the standard deviation is defined as and what is usually called the "uncertainty" $\Delta p$ in most physics texts.
Best Answer
For the specific case of the particle-in-a-box problem, the momentum operator is much trickier to handle than what you're allowing for, but that's not the issue at the heart of your current confusion.
If you know that two operators $A$ and $B$ commute, i.e. $[A,B]=0$, then there's two common ways to understand the consequences, one of which is right and one of which is wrong:
You are guaranteed the existence of at least one shared eigenbasis of both operators; but
You are not guaranteed that all eigenbases of $A$ will be eigenbases of $B$, or vice versa.
(For other occurrences of the misconception that the second point does happen, see e.g. this answer or this one.)
Thus, the 'contradiction' you observe is also present for the case of the free particle, without the box: the wavefunction $\psi(x)=\cos(kx)$ is an eigenfunction of $\hat H=\frac{1}{2m}\hat P^2$, but it is not an eigenfunction of $\hat P$. That is, the fact that the hamiltonian's eigenfunctions over a compact interval are not momentum eigenfunctions is not surprising and it is not specific to that configuration.
But, that said, you can still ask something like
and it's a reasonable question. Here what happens is that the subtleties with the compact-interval momentum operator kick in, at a much deeper level than just the commutation: the result, in full, reads
and it breaks because $\hat P$-in-a-box is not a self-adjoint operator: it is symmetric, but it has domain problems that prevent it from being self-adjoint. The consequences of this are as deep as they get: there simply isn't a momentum eigenbasis in this Hilbert space. You can extend the momentum operator to make it self-adjoint; this extension is not unique, but there is a reasonable choice (setting $\alpha=0$ in V. Moretti's answer) that comes close to making $\hat H$ the square of the extended $\hat P_\alpha$, but ultimately this can't work, as they have different domains. (More specifically, the domain of $\hat H$ is contained in the domain of $\hat P_\alpha$, but the eigenfunctions of $\hat P_\alpha$ do not fall into that subspace.)