Treating the problem more geometrically: There is a well-known partition of $\mathbb{R}^3$ into nested one-sheet hyperboloids of revolution, along with their common axis of symmetry; https://en.wikipedia.org/wiki/Skew_lines has an image. Using the $z$ axis as the axis of symmetry, the equation of a hyperboloid of revolution is $\frac{x^2+y^2}{a^2} - \frac{z^2}{c^2} = 1$ with $a \neq 0, c \neq 0$.

This is a family with two parameters and almost all pairs of members intersect. To obtain a family that partitions $\mathbb{R}^3$, fix a relationship between $a$ and $c$ so that every choice of $x^2+y^2 \neq 0$ and $z^2$ admits exactly one choice of $a$ and $c$, thus reducing the family to one parameter.

Each hyperboloid can be ruled in two ways. All of the lines of a ruling have the same slope relative to the $z$ axis, which is $\pm a/c$, with $+$ for one ruling and $-$ for the other. Clearly, all of the lines in a ruling are mutually skew, and to ensure the lines in different hyperboloids are mutually skew, the chosen relationship between $a$ and $c$ must not admit two solutions with the same $a/c$. A trivial choice is to fix $c=1$.

With the relationship between $a$ and $c$ fixed, there is an additional (uncountable) choice of which ruling of each hyperboloid.

This family of solutions to the problem can be extended by noting that all of the required properties are preserved under affine transformations of $\mathbb{R}^3$, as well as preserving the component hyperboloids as hyperboloids (but not necessarily as hyperboloids of revolution).

Projective geometry allows further extension: Skewness of lines is preserved by projective transformations. (In projective geometry, parallel lines become lines whose intersection is a point at infinity, so the union of parallel lines and intersecting lines is projectively invariant, despite that neither set alone is.) The extension of the above construction fills all of projective 3-space except the "equatorial" line which is the intersection of the $z=0$ plane of symmetry with the plane at infinity, because none of the lines are parallel to the $x=y=0$ plane. So the equatorial line must be added to the partition, providing a dual of some sort to the $z$ axis.

Any projective transformation can be applied to the extended construction, followed by restriction to $\mathbb{R}^3$, leaving a partition into skew lines (and possibly one transformed line is contained in the plane at infinity). In particular a projective transformation can be applied which maps the equatorial line into a line in $\mathbb{R}^3$ while mapping the $z$ axis into a line still in $\mathbb{R}^3$. (My memory is that given two pairs of skew lines, there is a projective transformation that maps the first pair into the second pair.)

This constructs a solution that is difficult to visualize. It has two skew "axis" lines, one mapped from the $z$ axis and one mapped from the equatorial line (the line at infinity with $z=0$). The hyperboloids are mapped into quadratic surfaces (because projective transformations preserver quadratics). Each axis is surrounded by nested one-sheet hyperboloids. I suspect that both families, going outward from their axis, expand more quickly on the side away from the other axis, and the families almost partition $\mathbb{R}^3$, with the limiting surface between them being a hyperbolic paraboloid.

I'm sure this was worked out in the 19th century, and I have a vague memory I've seen an image of this structure. Does anyone know a reference?

I suspect that looking at projective 3-space using homogeneous coordinates (which is more or less $\mathbb{R}^4$ modulo scalar multiplication, intersected with a 3-plane that does not contain the origin) is isomorphic to some of the above constructions using $\mathbb{C}^2$ intersected with $\mathbb{R}^3$.