Intuitive explanation of Duhamel's principle
Don't hurt your brain too much in trying to understand the physical interpretation. Start from the mathematics! Let me do the ODE version first:
Let $X: \mathbb{R}\to \mathbb{R}^N$ be a vector valued function of one (time) variable. Let $A:\mathbb{R}^N \to\mathbb{R}^N$ be a linear transformation that is independent of the variable $t$. Let $R_t$ denote the solution operator: that is $ X(t) = R_tX_0$ where $X_0$ is a fixed vector solves the equation $$ \dot{X} = AX $$ with initial data $X(0) = X_0$. Now let $Y(t)$ be another vector valued function. Consider the expression
$$ X(t) = \int_0^t R_{t-s}Y(s) \mathrm{d}s $$
Computing it explicitly we have that
$$ \left(\frac{d}{dt} - A\right)X(t) = \int_0^t\left(\frac{d}{dt} - A\right)\left[R_{t-s}Y(s)\right]~ \mathrm{d}s + R_{t-s}Y(s)\Big|_{s = t} $$
The second term comes from the fundamental theorem of calculus when the $\frac{d}{dt}$ derivative hits the integral sign.
By definition the term under the integral sign evaluates to 0, since $Z(t) = R_{t-s}Y(s)$ solves the homogeneous equation with initial data $Z(s) = Y(s)$. So we are left with
$$ (\frac{d}{dt} - A)X(t) = R_0Y(t) = Y(t) $$
i.e. that $X(t)$ solves the inhomogeneous problem.
For the PDE version, you just replace $\mathbb{R}^N$ with a Banach or Hilbert space, and the computation formally carries through in exactly the same way.
Coming back to the physical interpretation: you see that Wikipedia is being a bit imprecise: it is not the adding up of solutions that matters really; what drives Duhamel's principle is the fact that you are adding only a limited number of solutions.
Let me clarify: imagine you have a family of initial data $Y(s)$. And you write down the expression
$$ X(t) = \int_{a}^{b} R_{t-s}Y(s) \mathrm{d}s $$
which would be what we do if we were to just add (integrate) the contributions from all of the linear waves coming from the "inhomogeneity", by linearity it is clear that $X(t)$ will still solve the homogeneous equation, since it is a fixed (as in the limits of the integral) sum of may solutions.
The magic of Duhamel's principle is in that the upper-limit of the integral is time! That is, we defined
$$ X(t) = \int_0^t R_{t-s}Y(s) \mathrm{d}s $$
As you see from the derivation above, it is this upper-limit which, when acted on using the fundamental theorem of calculus, give you the inhomogeneous term. So what is the physical interpretation then? That the upper-limits also changes represents the fact that the solution at $X(t + \Delta t)$ consists of the forward time evolution of the solution at $X(t)$ plus a new contribution from the data in $(t,t+\Delta t)$ which was not included in the computation in $X(t)$. This "adding a new contribution" is precisely what we imagine the inhomogeneous term as, that is, a source term!