Euler-Lagrange equations for a damped mass-spring system. Is there a way to use Euler-Lagrange equations to derive the following formula?

I can derive the term not due to damping, I guess the damping force should be somehow modeled using a generalized force. Can you please show me how?
My attempt was starting by
$$
F_{d,ij} = k_d \lVert v_{ij} \rVert u_{ij}
$$
Where $v_{ij}$ is the relative velocity vector and $u_{ij}$ is the unit vector going from mass $i$ to $j$. But I'm missing how should I go from here to the generalized forces.
Thank you.
 A: *

*OP is essentially asking to find a velocity-dependent potential for coupled damped oscillators. This can be reduced to the question of finding  a velocity-dependent potential $U({\bf x},\dot{\bf x},t)$ for the 1-particle force
$$ {\bf f} ~=~-k_d{\bf x}\frac{{\bf x} \cdot\dot{\bf x}}{|{\bf x}|^2} 
~=~-k_d{\bf x}\frac{d\ln |{\bf x}|}{dt}.$$
One may show, using methods e.g. outlined in my Phys.SE answer here that such velocity-dependent potential $U({\bf x},\dot{\bf x},t)$ does not exist. So OP's model cannot be obtained from a stationary action principle (under the assumption that we shouldn't modify the kinetic sector of the action).

*However, it is still possible to describe OP's model via Lagrange equations
$$  \frac{d}{dt} \left(\frac{\partial L}{\partial \dot{\bf x}_i}\right)-\frac{\partial L}{\partial {\bf x}_i}~=~-\frac{\partial {\cal F}}{\partial \dot{\bf x}_i} $$
using the Rayleigh dissipation function$^1$
$$ {\cal F} ~=~ \frac{1}{2}\sum_{i<j}  k_d^{ij} \frac{({\bf x}_{ij} \cdot\dot{\bf x}_{ij})^2}{|{\bf x}_{ij}|^2}, \qquad {\bf x}_{ij}~:=~{\bf x}_{i}-{\bf x}_{j}. $$
The Lagrangian itself is
$$L~=~T-V, \qquad T~=~\frac{1}{2}\sum_i m_i |\dot{\bf x}_i|^2, \qquad V~=~\frac{1}{2}\sum_{i<j}k^{ij}_s (|{\bf x}_{ij}| - \ell^{ij}_0)^2. $$
See also this related Phys.SE post.
--
$^1$ In this answer, we have slightly generalized OP's model to allow the spring and dissipation constants to depend on the particle pair $(i,j)$.
A: Let us derive the undamped case. The kinetic energy $T$ and the potential energy $V$ are
$$
T = \frac{1}{2}m\left({\dot x_i}^2 + {\dot x_j}^2 \right) \qquad\text{and}\qquad V = \frac{1}{2}k_s (|x_j - x_i| - l_0)^2 \, .
$$
Therefore, the Euler-Lagrange equations $\frac{\text d}{\text d t}\frac{\partial \mathcal{L}}{\partial \dot{\boldsymbol x}} = \frac{\partial \mathcal{L}}{\partial {\boldsymbol x}}$ deduced from the Lagrangian $\mathcal{L} = T-V$ are
\begin{aligned}
m \ddot x_i &= \phantom{-}\frac{x_j - x_i}{|x_j - x_i|} k_s (|x_j - x_i| - l_0) = \phantom{-}f|_{k_d=0} \, ,\\
m \ddot x_j &= -\frac{x_j - x_i}{|x_j - x_i|} k_s (|x_j - x_i| - l_0) = -f|_{k_d=0} \, .
\end{aligned}
The dissipative term with damping constant $k_d\neq 0$ is not included in the present Euler-Lagrange equations, but it can be done upon introducing a dissipation potential (see other answers)
