Determining the distribution of sums of random variables is a popular topic in statistics, but even more so in finance and actuarial science.
So when you are looking for literature you can also check the literature in theses fields.
There are two cases that you can consider:
First, specific margins coupled with specific dependence structures, and second a general framework.
An example for the first point is the multivariate normal (or more generally multivariate elliptical distributions), where the distribution of sums can be computed. Here it is important to note that this is only possible for a specific combination of margins and dependence structures. Finding such particular combinations and calculating the resulting distribution of $X+Y$ is (or was) a popular topic in finance and actuarial science.
If you want to go in this direction with your Rayleigh margins, you would have to try out different families of copulas to see in which cases you can solve the resulting double integral analytically (possibly taking inspiration from the literature on how to do this).
As a result you would have a model $F_Z$ for $Z=X+Y$ that will depend on the parameters of $X$, $Y$ and the chosen copula $C$, and you can then go on to estimate these parameters from you data.
This approach is nice when it works, but since it is not always possible to find closed-form solutions a numerical alternative is the following approach.
Concerning the second point, determining the distribution of $Z=X+Y$ in general for $(X,Y)$ with copula $C$ is discussed in:
Gijbels, I. and Herrmann, K. (2014). On the distribution of sums of random variables with copula-induced dependence. Insurance: Mathematics and Economics, 59, pp.27-44.
(https://www.sciencedirect.com/science/article/pii/S0167668714000961)
Your integral representation is indeed the starting point of the paper, but you can actually transform the variables into $[0,1]^2$ and get rid of one of the integrals leaving you with a one-dimensional integral only.
This integral can be computed efficiently and you can then again estimate the parameters of $F_Z$ based on your data.
Now coming back to your specific questions:
- Simulation from copulas is covered by the answer here (and you can use the copula package in R to run the simulations): How to sample from a copula?
- How to choose the copula $C$ is more difficult. If you have access to a bivariate sample $(x_i,y_i)^{i=1}_{n}$ from $(X,Y)$ you can choose a parametric model and estimate the parameters, or you can use a non-parametric approximation (like the empirical copula) to get an estimated copula $\hat{C}$.