This is a question I ended up with while trying to make a program that would find an (at least locally-) minimal-energy configuration for a piecewise linear 1-dimensional object in $\mathbb{R}^3$ in relation to a practical knot theory type of question. If it helps at all for the purposes of context or intuition, the setup described here is essentially asking for the electrostatic force between two one-dimensional nonconducting rods, as seen from one endpoint of one of the rods, with uniform charge distribution.
Let $(x_0, y_0, z_0)$, $(X_0, Y_0, Z_0)$, $(x_1, y_1, z_1)$, and $(X_1, Y_1, Z_1)$ be four elements of $\mathbb{R}^3$, with the standard Euclidean norm $\lVert (x, y, z) \rVert = \sqrt{x^2 + y^2 + z^2}$. What is $$\int_{0}^{1} \!\! \int_{0}^{1} \!\! \frac{\displaystyle \big( (x_0, y_0, z_0) + (X_0 - x_0, Y_0 - y_0, Z_0 - z_0) u \big) - \big( (x_1, y_1, z_1) + (X_1 - x_1, Y_1 - y_1, Z_1 - z_1) v \big)}{\displaystyle {\left\lVert \big( (x_0, y_0, z_0) + (X_0 - x_0, Y_0 - y_0, Z_0 - z_0) u \big) - \big( (x_1, y_1, z_1) + (X_1 - x_1, Y_1 - y_1, Z_1 - z_1) v \big) \right\rVert}^3} \ \mathrm{d} v \ u \ \mathrm{d} u?$$
I'm pretty sure this should have a closed-form solution, but after a bunch of WolframAlpha, find-and-replace, and failing to apply manual partial integration correctly, I've only gotten rid of one of the integrations and the remaining one looks so messy I reckon I must have missed some important simplification somewhere, although it numerically comes out the same as the original double integral so I doubt I made a mistake in manipulation. With $a_0 = X_0 - x_0$, $b_0 = Y_0 - y_0$, $c_0 = Z_0 - z_0$, $a_1 = X_1 - x_1$, $b_1 = Y_1 - y_1$, and $c_1 = Z_1 - z_1$, a screenshot from Desmos of the result I have obtained so far is below.
EDIT: I realized that the form of the integral that I was solving in the screenshot below, as well as the one I originally asked about, simply had a $\mathrm{d} u$ term as opposed to the $u\ \mathrm{d} u$ I actually meant to be investigating, so it (as per the physical intuition) was measuring the total force over the entire line from $(x_0, y_0, z_0)$ to $(X_0, Y_0, Z_0)$ instead of just the force percieved by $(x_0, y_0, z_0)$. Just imagine it ends in $u\ \mathrm{d} u$ instead of just $du$. However, if the $\mathrm{d} u$ form is substantially easier to solve, that would be acceptable as well; ideally I would be able to test both forms, but $u\ \mathrm{d} u$ is I believe more relevant to the physical intuition.
UPDATE: After some messing around I have managed to determine that the integral in question here, for both the $\mathrm{d} u$ and $u\ \mathrm{d} u$ versions, is invariant under applying an isometric transformation to the values $(x_0, y_0, z_0)$, $(X_0, Y_0, Z_0)$, $(x_1, y_1, z_1)$, and $(X_1, Y_1, Z_1)$ and then applying the inverse map to the result of the integral, which allows me to assume WLOG that all of $x_0$, $y_0$, $z_0$, $Y_0$, and $Z_0$ are all zero, which turns the above integral into
after a $u$- (or I guess more properly a $v$- in this case) substitution of $v = X_0 u$, of course replacing the $\mathrm{d} v$ by $\frac{v}{X_0}\ \mathrm{d} v$ for the $u\ \mathrm{d} u$ form. This is a bit better, I think, but I doubt it actually helps much if at all toward getting a solution.