Solution (source code)

= Solution

In the bilinear identity choose field $a$ to be the <complex conjugate> of the field at <frequency> $\omega$, equivalently its $-\omega$ partner, and field $b$ to be the field itself. The squared <frequencies> are equal, so the <mass density> terms cancel:
$$
\int_{\mathcal D}(-\widehat f_i^*\widehat u_i+\widehat u_i^*\widehat f_i)dV=\int_{\partial\mathcal D}(t_i^*\widehat u_i-\widehat u_i^*t_i)dS.
$$
Multiply by $i\omega$ and use $\widehat v_i=-i\omega\widehat u_i$, $\widehat v_i^*=i\omega\widehat u_i^*$ to obtain
$$
\boxed{\int_{\mathcal D}(\widehat f_i^*\widehat v_i+\widehat v_i^*\widehat f_i)dV=\int_{\partial\mathcal D}(-t_i^*\widehat v_i-\widehat v_i^*t_i)dS.}
$$
For the real fields defined in this question as the sum of positive and negative <frequency> components, time averaging cancels the $e^{\pm2i\omega t}$ terms. Their mean body-force <work> <mass density> is exactly $\widehat f_i^*\widehat v_i+\widehat v_i^*\widehat f_i$, and their outward mean <elastic-wave energy flux> is $-t_i^*\widehat v_i-\widehat v_i^*t_i$. Thus \b[mean body-force <power> equals mean outward boundary <power>]. Mean stored <energy> does not accumulate in a periodic steady field.

There is no extra factor $1/2$ here: the specified real field is $\widehat f e^{-i\omega t}+\widehat f^*e^{i\omega t}$, not $\operatorname{Re}(\widehat f e^{-i\omega t})$. For the latter phasor convention, both averages instead have the usual factor $1/2$. At $\omega=0$ <velocity> and both powers vanish.