I suggest that it doesn't make much sense to say that the planets orbit the barycenter of the solar system.
Beware that you are going very much against the grain of the best models of the solar system in writing that. All three of the leading ephemeris models (JPL's Development Ephemeris, the Russian Institute for Applied Astronomy's Ephemerides of the Planets and Moon, and the IMCCE's INPOP) use a barycentric rather than heliocentric coordinate system to model the motions of the bodies that comprise the solar system.
That said, the following plot shows the distances between Venus and the Sun (red) and Venus and the solar system barycenter (black) from January 1970 to December 2014. The horizontal (time) axis is in days from 12 Noon TT, 1 January 2000.
I used Venus rather than the Earth for two reasons: The Earth's orbit is a bit more complex than is Venus' orbit because of the Moon, and because Venus' orbit is the closest to circular of all the planet,s so variations should be small. The red curve, the distance between Venus and the Sun, is very close to a small amplitude sinusoidal. The black curve exhibits beats and other nastiness. A much simpler picture results when Venus is modeled as orbiting the Sun as opposed to the solar system barycenter.
So why do solar system modelers inevitably use a barycentric rather than heliocentric coordinate system? First, I'll digress a bit and look at the question "Do the planets orbit the Sun or the solar system barycenter?" This is posed as an either-or question, implying that only one point of view is valid. Both points of view are equally valid. All frames of reference are equally valid. Do the math right, and even a Nix-centered, Nix-fixed coordinate system would work just fine, at least for a little while. One problem with a Nix-centered, Nix-fixed frame is that Nix's rotational motion is chaotic (see https://www.youtube.com/watch?v=zwSFC-aPEG0). Another is that it just doesn't make much sense, physically. But it can be done!
Solar system modelers have moved beyond using perturbed ellipses to describe the behaviors of the planets. Numerical integration provides more precise and faster predictions, so long as a frame with simple dynamics is used. A barycentric frame is preferred because that is the frame in which the equations of motion take on their simplest form.
Short answer: Use the radial acceleration component $a_\mathrm r$ of the total $\bf a$ by the central force and then deduce the SHM equation using the properties of circular orbital motion (effective potential energy curves can also be used) that $\dot r= 0$ and keeping in mind the fact that SHM would only occur when the circular orbit is stable.
Here the approach is elaborately discussed:
The radial component of the central force $\bf F$ is given by
$$\mathbf F\cdot \mathbf e_\mathrm r~=~ m\left[\ddot r - r(\dot \theta)^2\right]\;.\tag 1 $$
Due to the small perturbation, the radius of the circular orbit changed from $r_0$ to $r(t) ~=~r_0 + \rho(t)~~_; ~~~\rho(0)\ll r_0\;.$
Using $\mathit l= mr^2\dot \theta=\textrm{const}$, $(1)$ becomes
$$\ddot\rho(t)- \frac{l^2}{m^2(r_0+ \rho(t))^3}~=~ \frac{F_\mathrm r(r_0+ \rho(t))}{m}~\tag {1-a} $$
Now, \begin{align}(r_0+ \rho(t))^{-3}&=r_0^{-3}\left(1+ \frac{\rho}{r_0}\right)^{-3}\\ &\approx r_0^{-3}\left(1-3\frac{\rho}{r_0}\right)\\ F_\mathrm r(r_0+\rho)&= F_\mathrm r(r_0) +\frac{\mathrm dF_\mathrm r}{\mathrm dr}\bigg|_{r=r_0}~\rho+ \ldots \\ &\approx F_\mathrm r(r_0) +\frac{\mathrm dF_\mathrm r}{\mathrm dr}\bigg|_{r=r_0}~\rho\end{align}
So,
$$\ddot\rho- \frac{l^2}{m^2r_0^3}\left(1-3\frac{\rho}{r_0}\right)= \frac{F_\mathrm r(r_0)}m +\frac{\mathrm dF_\mathrm r}{\mathrm dr}\bigg|_{r=r_0}~\cdot \frac{\rho}m\tag{1-b} $$
From the orbit-equation$^\ddagger$, we get
$$\mathit l^2=- mr_0^3 F_\mathrm r(r_0)_; $$
So,
$$\ddot \rho+ \left[\frac{3\mathit l^2}{m^2r_0^4}- \frac{1}{m}\frac{\mathrm dF_\mathrm r}{\mathrm dt}\bigg |_{r=r_o}\right]\rho~=~0\;.$$
This can be written as
$$\ddot \rho + \omega^2 \rho ~=~0\;.\tag 2$$
Now, this is a second-order differential equation for simple harmonic oscillator with frequency $\omega\;,$
where $$\omega^2 ~=~ \frac{3\mathit l^2}{m^2r_0^4}- \frac{1}{m}\frac{\mathrm dF_\mathrm r}{\mathrm dt}\bigg |_{r=r_o}\;.$$
Unless, the circular orbit is unstable, a simple harmonic radial oscillation about $r=r_0$ would ensue (i.e. SHM would occur if $\omega^2\gt 0\;.$)
$\ddagger$ We would derive the orbit-equation from $(1)$ expressing $r$ and its derivatives in terms of $\theta\;.$
\begin{align}\dot r &= \frac{\mathrm dr}{\mathrm d\theta}\frac{\mathrm d\theta}{\mathrm dt}\\ \ddot r&=\frac{\mathrm d}{\mathrm dt}\left(\frac{\mathrm dr}{\mathrm d\theta}\dot \theta\right)\\&=\frac{\mathrm d^2r}{\mathrm d\theta^2}~\dot\theta ^2 + \frac{\mathrm dr}{\mathrm d\theta}~\ddot \theta \end{align}
Let
$$\mathbf H \stackrel{\text{def}}{=}
\mathbf r\times \dot{\mathbf r} =\frac{\mathbf L}m \;.$$
Therefore
\begin{align}\dot \theta &= \frac{|\mathbf H|}{r^2}\\ \implies \ddot \theta &=-\frac{2|\mathbf H|}{r^3} \, \dot r\;.\end{align}
Putting all these in $(1),$ we get
\begin{align}\frac{\mathrm d^2r}{\mathrm d\theta^2}~\dot\theta ^2 + \frac{\mathrm dr}{\mathrm d\theta}~\ddot \theta- r(\dot \theta)^2& = \frac{F_\mathrm r}{m} \\ \implies\frac{\mathrm d^2r}{\mathrm d\theta^2}~\left(\frac{|\mathbf H|}{r^2}\right)^2 +\frac{\mathrm dr}{\mathrm d\theta}~\left(-\frac{2|\mathbf H|}{r^3}\dot r\right) - r ~\left(\frac{|\mathbf H|}{r^2}\right)^2 &= \frac{F_\mathrm r}{m}\\ \implies \left(\frac{|\mathbf H|}{r^2}\right)^2\left[\frac{\mathrm d^2r}{\mathrm d\theta^2} - 2\frac{\left(\frac{\mathrm dr}{\mathrm d\theta}\right)^2}{r}-r\right]&=\frac{\mathrm F(r)}{m}\\ \implies \frac{\mathrm d^2r}{\mathrm d\theta^2} - 2\frac{\left(\frac{\mathrm dr}{\mathrm d\theta}\right)^2}{r}-r&= \frac{r^4 F_\mathrm r}{|\mathbf H|^2 m}\;.\tag{i}\end{align}
$(\rm i)$ is the orbital equation.
Now, applying the facts that for circular orbits, $\dot r,~\ddot r~=~0\,,$ we get the value of $$|\mathbf H|~=~\sqrt{ -\frac{r^3F_\mathrm r}{m}}\;,$$ where $r~=~r_0\;.$
Best Answer
As Guillermo Angeris correctly pointed out, this is essentially a numerical roundoff problem, not a physical situation.
As a physical example, there are sungrazing comets that get very close to to Sun, yet they maintain their original elliptical (or hyperbolic) orbit, without the orbit precessing a full third of a circle as you seem to be seeing.
Computationally, there are a few interesting issues. As Kyle pointed out in a comment, many integration schemes are indeed unreliable in that roundoff error (which is always present in floating-point computations) can accumulate in a runaway feedback. Indeed I often advise using leapfrog methods over Euler (used by Box2D) or even Runge-Kutta (see for instance What is the correct way of integrating in astronomy simulations? over at the Computational Science Stackexchange).
However, I suspect your problem is even simpler, in the sense that even an unstable numerical scheme should work for one or two orbits. Given that everything is going wrong in just one pass, it seems that your timesteps are simply too large. A brief glimpse at the Box2D documentation suggests you don't change the timestep mid-simulation, so I presume you are just using a good value to simulate the whole process in reasonable time. The problem is that when gravitating bodies get close in their orbit, they move quickly, sometimes very quickly. The way the code works is it updates each object's position and velocity at each timestep, where the new velocity is determined by the force. As far as I can tell, this is done in line 206 of
b2Island.cpp(v. 2.2.1):Without looking at your code, I am guessing you simply calculate the gravitational force the body should feel at that moment, and have the simulation chug away. The problem is this moves the orbiting object in a straight line for the next timestep, and that straight line takes it too far away from the gravitating mass for that mass to properly curve its orbit into a closed ellipse. The quick schematic below shows the blue object moving to the tip of the red arrow, rather than staying on the path.
Physically, your timestep should be smaller than any timescale you encounter in the problem. Now for an orbit conserving angular momentum, the product of the orbiting body's mass, tangential velocity $v$, and distance from the other object $r$ should be constant: $v \sim 1/r$. At the same time, the acceleration $a$ it feels is given by Newton's law of gravity: $a \sim 1/r^2$. So one natural timescale in this problem is $$ t \sim \frac{v}{a} \sim \frac{1/r}{1/r^2} \sim r $$ (omitting dimensional constants), which goes to show that if your timescale is just barely small enough and then you tweak the orbit so as to half the periapsis distance (distance of closest approach), then you would expect to need timesteps at least twice as small in order to preserve the integrity of the simulation.