Engineers often insist on using locally conservative methods such as finite volume, conservative finite difference, or discontinuous Galerkin methods for solving PDEs.
What can go wrong when using a method that is not locally conservative?
Okay, so local conservation is important for hyperbolic PDEs, what about elliptic PDEs?
In the solution of nonlinear hyperbolic PDEs, discontinuities ("shocks") appear even when the initial condition is smooth. In the presence of discontinuities, the notion of solution can only be defined in the weak sense. The numerical velocity of a shock depends on the correct Rankine-Hugoniot conditions being imposed, which in turn depends on numerically satisfying the integral conservation law locally. The Lax-Wendroff theorem guarantees that a convergent numerical method will converge to a weak solution of the hyperbolic conservation law only if the method is conservative.
Not only do you need to use a conservative method, in fact you need to use a method that conserves the right quantities. There's a nice example that explains this in LeVeque's "Finite Volume Methods for Hyperbolic Problems", Section 11.12 and Section 12.9. If you discretize Burgers' equation
$$u_t + 1/2 (u^2)_x = 0$$
via the consistent discretization
$$U^{n+1}_i = U^n_i - \frac{\Delta t}{\Delta x} U^n_i (U^n_i-U^n_{i-1})$$
you will observe that shocks move at the wrong speed, no matter how much you refine the grid. That is, the numerical solution will not converge to the true solution. If you instead use the conservative discretization
$$U^{n+1}_i = U^n_i - \frac{\Delta t}{2\Delta x} ( (U^n_i)^2-(U^n_{i-1})^2)$$
based on flux-differencing, shocks will move at the correct speed (which is the average of the states to the left and the right of the shock, for this equation). This example is illustrated in this IPython notebook I wrote.
For linear hyperbolic PDEs, and for other types of PDEs which typically have smooth solutions, local conservation is not a necessary ingredient for convergence. However, it may be important for other reasons (e.g., if the total mass is a quantity of interest).
I think one answer to your question is that certain communities simply always used conservative schemes and so it has become part of "the way it's done". One may argue whether that's the best way to do it, but that's about as fruitful as asking the British to drive on the right because it would simply be more convenient to have only on standard side.
That said, I do see cases where it's useful. Think, for example, of two-phase porous media flow. This problem is commonly posed in the following way: $$ u + K \nabla p = 0 \\ \nabla \cdot u = f \\ \partial_t S + u \cdot \nabla S = q. $$ Here, part of the problem is solving the mixed Laplace that makes up the first two equations, a task traditionally done using Raviart-Thomas elements. They are often chosen because of "the importance of ensuring mass conservation", and in a sense I can understand that: if you end up with a velocity field that's not mass conservative, you'll get a saturation equation that does not conserve the overall mass of the transported fluid. Of course one can argue that that wouldn't be so bad because it would all be the same in the limit $h\rightarrow 0$, but the insistence on ensuring that this property holds even for finite mesh sizes does make some sense.
Many times, the equations to be solved represent a physical conservation law. For example, the Euler equations for fluid dynamics are representations of conservation of mass, momentum, and energy. Given that the underlying reality that we are modeling is conservative, it is advantageous to choose methods that are also conservative
You can also see something similar with electromagnetic fields. Maxwell's laws include the divergence-free condition for the magnetic field, but that equation is not always used for the evolution of the fields. A method that conserves this condition (for example: constrained transport) helps match the physics of reality.
Edit: @hardmath pointed out that I forgot to address the "what could go wrong" part of the question (Thanks!). The question specifically refers to engineers, but I'll provide a few examples from my own field (astrophysics) and hope that they help illustrate the ideas enough to generalize to what might go wrong in an engineering application.
(1) When simulating a supernova, you have fluid dynamics linked to a nuclear reaction network (and other physics, but we'll ignore that). Many nuclear reactions depend strongly on the temperature, which (to a first-order approximation) is some measure of the energy. If you fail to conserve energy, your temperature will be either too high (in which case your reactions run much too fast and you introduce far more energy and you get a runaway that shouldn't exist) or too low (in which case your reactions run much too slow and you can't power a supernova).
(2) When simulating binary stars, you need to recast the momentum equation to be conservation of angular momentum. If you fail to conserve angular momentum, then your stars cannot orbit each other correctly. If they gain extra angular momentum, they separate and stop interacting correctly. If the lose angular momentum, they crash into each other. Similar issues occur when simulating stellar disks. Conservation of (linear) momentum is desirable, because the laws of physics conserve linear momentum, but sometimes you have to abandon linear momentum and conserve angular momentum because that is more important to the problem at hand.
I have to admit, despite citing the divergence-free condition of magnetic fields, I'm not as knowledgeable there. Failure to maintain the divergence-free condition can generate magnetic monopoles (which we have no evidence of currently), but I don't have any good examples offhand of issues that might cause in a simulation.
Today I come across a thesis The EMAC Scheme for Navier-Stokes Simulations, and Application to Flow Past Bluff Bodies and notice Section 1.2 of it answers OP's question, at least partially. The relevant parts are:
It is widely believed in the computational fluid dynamics (CFD) community that the more physics is built into the discretization, the more accurate and stable the discrete solutions are, especially over longer time intervals. N. Phillips in 1959 [42] constructed an example for the barotropic nonlinear vorticity equation (using a finite-difference scheme), where the long-time integration of the convection terms results in a failure of numerical simulations for any time step. In [4] Arakawa showed that one can avoid instability issues with integration over long time if kinetic energy and enstrophy (in 2D) are conserved by a discretization scheme. … . In 2004, Liu and Wang developed that conserves helicity and energy for three-dimensional flows. In [35], they present an energy and helicity-preserving scheme for axisymmetric flows. They also show that their dual conservation scheme eliminates the need for large nonphysical numerical viscosity. …
… It has been known for decades in CFD, that the more physical quantities are conserved by a finite element scheme, the more accurate the prediction, especially over the long time intervals. Thus the solutions provided by a more physically accurate scheme are also more physically relevant. If one could afford a fully resolved mesh and infinitely small time step, all commonly used finite element schemes are believed to provide the same numerical solutions. However, in practice one cannot afford a fully resolved mesh in 3D-simulations, especially for time-dependent problems. For example in chapter 2 we need 50-60 thousand time steps, where each time step requires solving a sparse linear system with 4 million unknowns. This required 2-3 weeks of computational time with highly parallel code on 5 nodes with 24 cores each.
