Apart from the obvious reason, that knowledge is a value in itself, there are indeed practical engineering application of multivariate calculus about order $n \geq 2$.
One example is the study of stability, which clearly benefits from Taylor expansions up to second order.
There are examples from computational methods too: if the engineer is using a Finite Element code and at some point along an analysis gets a "non positive definite Hessian", he would certainly benefit from knowing about second-order expansions to understand what is going on.
Another fitting example in my view is given by Laplace's method (alternatively, saddle-point) to approximate the solution of certain integrals, very common in Statistical Mechanics and not rare at all in Engineering. The method mandatorily requires a second order expansion.
Talking about thermodynamics, how would one prove the principle of minimum energy to an engineer not accustomed to second order expansions? With multivariate Taylor expansions, it takes just a moment, by expanding the entropy $S$ as function of an internal parameter $X$ and energy $U$ to second order, as $ \frac{\partial S}{\partial X} = 0$
$$ \mathrm{d}S \approx \frac{\partial^2 S}{\partial X^2}\mathrm{d}X^2 +\frac{\partial S}{\partial U} \mathrm{d} U $$ and the minimum nature of $U$ is on display.
Moreover, referring to elasticity as in your post, the elastic moduli are related to the strain energy function $W$ via $$ \mathcal{C} = \frac{\partial^2 W }{\partial E^2} $$ where $E$ stands for a deformation tensor, so a second-order expansion of the strain energy is needed.
In $1$ D the "first-order only" engineer could avoid the problem by looking at stress-strain curves only. Indeed, in $1D$ $$ \sigma = \frac{\mathrm{d} W}{\mathrm{d} \epsilon}$$ where $\epsilon$ stands for the strain. As for small enough strains $\sigma = E \epsilon$, in this case only first-order condsiderations are needed to estimate the elastic moduli
(in fact, by letting the experiment perform a differentiation).
But in a more general case, dealing for example with a biaxial test (and this is not an academic example for modern engineers busy with polymers), it can often be easier to start from the energy $W$, perform a multivariate expansion up to second order, and get the elastic moduli matrix out of it.
The difference stems from practical reasons: in a $1D$ test one stretches in one direction, and what happens in the other directions is easy to characterise.
In a biaxial test this is far more complex: it is difficult to impose one strain value along a given direction, while the strain along the perpendicular direction is varied. In other words, it is difficult to "perform partial derivatives" experimentally.
All this could be especially useful, and releavant, when a component is loaded by a constant stress state, on top of which small oscillations are superimposed. It is then useful to re-define the constant stress-state as the "unloaded" state, and calculate "fictitious" elastic moduli.
Best Answer
Try this:
$$\sqrt{(1+k+\epsilon)^2-4k\epsilon} = \sqrt{(1+k)^2 + 2(1+k)\epsilon + \epsilon^2 -4k\epsilon} = \sqrt{(1+k)^2 + 2(1-k)\epsilon + \epsilon^2} = (1+k)\sqrt{1 + \frac{2(1-k)\epsilon + \epsilon^2}{(1+k)^2}}$$
and the first order expansion is
$$(1+k)(1 + \frac{(1-k)\epsilon}{(1+k)^2}) = (1 + k) + \frac{1-k}{1+k}\epsilon$$