I haven't done much astronomical image processing before, but as this question is unanswered I'll give it a shot - hopefully to some avail. If the problem is more specific, a code sample/image sample would probably be useful for further diagnosis, but otherwise this example may help. It discusses the process of writing a 3-channel image to separate FITS images. I would try this first, and check if this output is also garbage - it might just be an issue with how you're using the HDUList.writeto() function. The relevant code is pasted below (with a couple of edits, as the example uses the Pillow Image class at first):
Split the three channels (RGB):
r, g, b = img[:, :, 0], img[:, :, 1], img[:, :, 2]
Write out the channels as separate FITS images
red = fits.PrimaryHDU(data=r)
green = fits.PrimaryHDU(data=g)
blue = fits.PrimaryHDU(data=b)
Also out of curiosity - why is your "numpy RGB array" storing values between 0 and 512? The standard for images is usually integers between 0-255 or floats between 0 and 1; but my guess is that astronomy requires a larger bit depth? In any case, you should also make sure the dtype of your numpy array is large enough for the desired bit depth - often images are stored with dtype="uint8" which ranges from 0 to 255.