Implementation of a Variable Fractional Delay with Lagrange Interpolation using Farrow Structure
The task at hand is a pretty common one; most symbol timing synchronizers need to do exactly what you describe.

Same goes for many types of resamplers – they need to reconstruct varying fractional inter-sample "snapshots" from the input sample stream to construct a non-rationally related rate output.

So that's where I'd start looking.

In fact, the farrow structure idea takes you pretty far!

Imagine you could have a large bank of linear-phase filters. A linear-phase filter has a constant group delay.

By modulating said filters, you can change the phase slope and hence the group delay. Now, you just pick the filter that fits your needs.

In a polyphase filterbank, such as implemented by the GNU Radio "PFB Arbitrary Resampler" block, exactly that is done: You construct a whole bunch of filters, and for each sample, pick the one that's closest to what you need (or the two closest, and linearly interpolate between these).

Because it's used as a resampler, for every non-1.0 resampling ratio, the next sample will be choosen from a different filter branch, with the appropriate delay.

Now, since you, most of the time, just want to delay, you seldom need the "resampling" capability of that resampler. (Resampling is exactly what happens when you change the delay linearly)

So, instead, you'd have a single filter (if possible, use a filter that you need to have, anyway), and you'd modulate it.

Think about this in frequency domain: assume you had a linear, continuous Fourier transform that you apply to your signal and you multiply that with both the FT of your filter impulse response and $e^{j2\pi\tau f}$ before transforming back to time domain. To adjust delay, you adjust $\tau$. Tadaaah! Variable fractional delay.

Problem certainly is that you're in DSP, so you don't get that linear Fourier transform to implement linear convolution; you get the DFT, but multiplication in discrete frequency domain is equivalent to circular convolution in time domain.

But that's been solved already – fast convolution allows you to apply long filters with little effort (and low numerical problems – the FFT is exceptionally stable, usually even more stable than the direct convolution!). Now, if you go ahead and implement that e.g. via the overlap-save method, all you need to do is make sure that the phase of the complex exponential doesn't "jump" between FFT segments. You can then change (even gracefully!) the $\tau$ between individual convolutions.