Take the three equations in AnonSubmitter85's answer. Let's call these equations (1), (2), and (3). Now form three new equations by subtracting: (1) - (2), (2) - (3), and (3) - (1). The important thing is that the three new equations will not longer have any quadratic terms; they will just be linear in $x,y,z$, and therefore fairly easy to solve.
The three linear equations represent three planes, of course. What three planes are these? Well, if you intersect any two of the spheres, you will get a circle, as you said. Doing this two-at-a-time intersection will give you three circles. The three planes mentioned above are the planes containing these three circles. So, you were somewhat on the right track.
Edit: Despite the fact that it got two upvotes, what I wrote above is not really correct. It can't possibly be correct because intersecting the three planes will give only a single point (usually), whereas the intersection of three spheres will give two points (usually).
The three planes will actually intersect in a line, not a point. That's why you get NaN if you try to solve the system of three linear equations. You should take two of the planes and intersect them to get a line. Then intersect this line with a sphere to get the two desired points. My apologies. Moral: don't believe everything you read on the internet :-)
The tedioius details are worked out in this answer.
Also, good answers here, including some with code.
There is software that will compute the intersection (or union) of two closed triangle meshes as another closed triangle mesh. In fact, I wrote a program that reliably computes arbitrary triangle mesh intersections and unions.
Once you have a closed triangle mesh of the intersection region, there is in fact a very nice formula for the volume. The formula is obtained by using the divergence theorem:
$$
\int\int\int_V (\nabla\cdot\mathbf{F}) \operatorname{d}\!V = \int\int_S (\mathbf{F}\cdot \mathbf{n}) \operatorname{d}\!A \tag{1}
$$
where $\mathbf{F}:\mathbb{R}^3\rightarrow\mathbb{R}^3$ is any $C^1$ vector field and $\mathbf{n}$ is the exterior pointing unit surface normal. Define $\mathbf{F}(x,y,z) = (x,y,z)$. Then the left side of $(1)$ is just $3V$ (three times the volume of the region enclosed by the mesh). To compute the surface integral for a single triangle $ABC$, the barycentric coordinates surface parameterization works very well. The barycentric coordinates for $A$ are $(u,v)=(1,0)$, for $B$ they are $(u,v)=(0,1)$ and for $C$ $(0,0).$ Barycentric coordinates map $(u,v)$ to $(x,y,z)$ like this:
$$
(x,y,z) = uA + vB + (1-u-v)C = u(A-C) + v(B-C) + C. \tag{2}
$$
Since $\mathbf{n}$ is orthogonal to all vectors in the plane of the triangle, $\mathbf{F}\cdot \mathbf{n}$ is very easy to compute: it's just $C\cdot\mathbf{n}$ or equivalently $A\cdot\mathbf{n}$ or $B\cdot\mathbf{n}$ since $(A-C)\cdot\mathbf{n} = (B-C)\cdot\mathbf{n} = 0.$ Using $(2)$ it is easy to compute $\mathbf{x}_u = A-C$ and $\mathbf{x}_v = B-C.$ So the surface integral is
$$
\begin{align}
\int\int_S (\mathbf{F}\cdot \mathbf{n}) \operatorname{d}\!A &= \int_0^1\int_0^{1-v} (C\cdot\mathbf{n}) |(A-C) \times (B-C)| \operatorname{d}\!u \operatorname{d}\!v \\
& = (C\cdot \mathbf{n})|(A-C) \times (B-C)| / 2 \\
& = (C\cdot ((A-C)\times (B-C))) / 2.
\end{align}
$$
The last equality above is true since $\mathbf{n} = ((A-C)\times (B-C)) / |(A-C)\times (B-C))|.$ So the formula for the volume of a closed triangle mesh is just
$$
V = \frac{1}{6}\sum_{i=1}^N (C_i\cdot ((A_i - C_i)\times (B_i - C_i))) \tag{3}
$$
where $A_i$, $B_i$, and $C_i$ are the vertices of the $i$th triangle of the mesh and $N$ is the number of triangles in the mesh.
I use the formula $(3)$ all the time. You could test it on some sphere meshes. Just make sure that the face vertex order is consistent and produces an exterior facing normal. This means that if a triangle of your mesh is parallel to the screen and the exterior pointing normal pointing directly at you, the vertices $A$, $B$, and $C$ of the triangle must have counter-clockwise orientation.
Best Answer
The following holds for a specific case where the volume enclosed by 3 spheres is contained within the triangle formed by the 3 sphere centers.
I am hesitant to put this answer up even though I feel it is correct, because the closed form I have found is actually quite ugly. The derivation is quite pretty but the resulting integral is ugly. I will discuss it a bit at the end.
I will post my thoughts anyways in the hope that someone can come along and simplify these formulas into a nicer form, or perhaps find a much simpler formula altogether!
Derivation
Let $T$ be the solid triangle defined by sides of length $a,b,c$. Define 3 disks of raius $r$: $D,E,F$ and let each vertex of $T$ be the center of a disk in the natural way. See figure included below (at the very bottom) for a graphic description.
Let the wedges be
$$W_1 := D \cap T$$ $$W_2 := E \cap T$$ $$W_3 := F \cap T$$
Let the lenses be
$$L_1 := D \cap E$$ $$L_2 := E \cap F$$ $$L_3 := F \cap D$$
and let the rounded triangle be
$$ R := D \cap E \cap F .$$
It should be noted that each $L_i$ is bisected by a side of $T$. That is
$$ \left| L_i \cap T \right| = \left| L_i \cap \bar{T} \right| $$ where $\bar{T}$ denotes the compliment of $T$ and $| \cdot |$ denotes area. So naturally, define the half-lenses
$$H_1 := T \cap L_1 $$ $$H_2 := T \cap L_2 $$ $$H_3 := T \cap L_3 $$
Let $$ f(r) = \left| R \right| $$, it follows that the volume $V$ contained in the overlap of 3 spheres with equal radii $r_s$ is equal to
$$ V = 2 \int _0 ^{m} f\left(\sqrt{r_s^2-h^2}\right) dh$$
where $m = \sqrt{r_s^2 - \left( \frac {\max\{a,b,c\}} {2} \right)^2}.$
It remains to find a closed form of $f(r)$. We can find one applying the inclusion exclusion principle as
$$ f(r) = |T| - \sum_i |W_i| + \sum_i |H_i| $$
We have $$\sum |W_i| = \frac{\pi}{2} r^2$$
since the sum of the internal angles of $T$ is $\pi.$
By Heron's elegant formula we have
$$ |T| = \sqrt{s(s-a)(s-b)(s-c)}$$ where $s = \frac{1}{2} (a+b+c).$
Lastly, we can compute the area of each half-lense with
$$|H_i| = r^2 \cos^{-1} \left( \frac{d_i}{2r} \right) - \frac{1}{4} d_i \sqrt{4r^2-d_i}$$ where $d_i = c,b,a$ for $i=1,2,3$ using the formulas for lenses but halved.
Discussion
This seems to be a closed form, as the most difficult integrals to compute, namely those corresponding to the $|H_i|$ terms are found to have a closed form as per a quick query here and here on wolfram's alpha.
But because these indefinite integrals are so ugly, I fear this closed form has little practical use unless a simplification or other formula can be found!
Though one may still be able to find a tight upper and lower bound by simplifiying/approximating these ugly integrals.