what you calculated is the purely electric energy of an ideal (infinitely small) dipole. It is more usually written as $U = -\vec E\cdot \vec p$ with $\vec E$ the external field and $\vec p$ the dipole moment. This is useful as it can be applied directly to rigid dipoles (like some molecular dipoles). For example, you can recover using the principle of virtual work the force ($\vec F = \vec \nabla \vec E\cdot \vec p$) and torque ($\vec M=\vec p \times \vec E$) applied on the dipole by an external field. For example, at a fixed location, this means the dipole wants to align with the field. You can also calculate thermodynamic properties leading to temperature dependent susceptibility, the electric analogue of Curie paramagnetism, or some Van der Waals forces (Keesome forces).
However, for a dielectric, you need to include the internal energy which results from the dipole's deformation due to the external field. In the case of a linear dielectric dipole for example, this becomes $U = -\frac{1}{2}\vec E\cdot \vec p$. The $1/2$ comes from the quadratic form of energy, and note the sign change. Similarily, you can compute the force, and note that if the suscepibility is positive (as in most cases), this means that the dipole is attracted to regions of intense field. You can also calculate new Van der Waals forces (Debye force).
Hope it helps and tell me if you find some mistakes.
Edit:
Just as a mathematical clarification, a quicker derivation of your formula is by using the voltage $V$. Writing $q,-q$ the charges of the points at position $\vec r_1,\vec r_2$, you get $U = q(V(\vec r_1)-V(\vec r_2)) = -\vec E \cdot \vec p$ with $\vec p = q(\vec r_1-\vec r_2)$ and using a linear approximation.
For the linear dielectric case, a simple way to calculate the energy is to notice that the expression of the force $\vec F = \vec p\cdot \vec \nabla \vec E$ is always true (by the same argument as above this time using the field instead of the voltage), and since $\vec p$ is proportional to $\vec E$, you get $\vec F = -\vec \nabla(-\frac{1}{2}\vec p\cdot \vec E)$.
You can also calculate explicitly the elastic energy and add it to the previously calculated purely electric energy to get the same result, you'll have to introduce an additional elastic constant.
Edit:
Sorry, I went off-topic. The internal energy of a dipole diverges in the ideal case, similarly to the self energy of a point charge. You could have just used directly the electrostatic energy to get $U = -\frac{q^2}{4\pi\epsilon_0d}$ (you assume the force constant which is why you get the wrong sign). I don't know any applications of this, since you typically don't want to consider self energies, even when they are finite, they cancel out in energy differences being constant. Furthermore, for the link with dielectrics, you look at the response of induced charges due to the free charges so it is not the internal energy you are interested in. The only context that I can think of where it might be relevant is in particle physics where you calculate radiative correction, but for the relevant accuracy, a full quantum mechanical formalism is necessary.
