I think the question is assuming that each individual woman has an 82% chance of getting married, independently of what other women will do.
We aren't choosing a 2- or 3-woman sample, we are merely checking the marital status of 20 women and checking if there happen to be 2 or 3 who are unmarried.
EDIT: Another way of looking at the problem:
Let's say we have a ball pit filled with 1 million balls. 820,000 are blue and 180,000 are red. Therefore, if I pick a ball at random, I have an 82% chance of it being blue and a 18% chance of it being red.
Now, what if draw a blue ball, throw that ball away, and decide I want to draw another one? It's true that the probability distribution has changed, since there are now 819,999 blue balls and 180,000 red balls, with 999,999 total balls. But for simplicity's sake, we can assume the probability distribution it hasn't changed very much (only by ~$10^{-6}$ in fact), so keeping our 82%/18% distribution is still going to be mostly accurate.
If I draw a small number of samples relative to the total number of balls (~20 samples relative to 1 million), the distribution is approximately binomial.
So on a mathematical level, you are correct: the distribution does change when you sample without replacement, but I think the problem wants you to make a simplifying assumption.
In general, if $p_1, p_2, \ldots, p_m \in (0,1)$ are IID realizations from some probability distribution with mean $\mu_p$ and standard deviation $\sigma_p$, and $n_1, n_2, \ldots, n_m \in \mathbb Z^+$ are IID realizations of from another probability distribution with mean $\mu_n$ and standard deviation $\sigma_n$, and for each $i = 1, 2, \ldots, m$, we have random variables $$X_i \sim \operatorname{Binomial}(n_i, p_i),$$ and we are interested in the distribution of $S = \sum_{i=1}^m X_i$, then we have by linearity of expectation $$\operatorname{E}[S] = \sum_{i=1}^m \operatorname{E}[X_i].$$ In turn, for each $X_i$, we have by the law of total expectation $$\operatorname{E}[X_i] = \operatorname{E}[\operatorname{E}[X_i \mid (n_i \cap p_i)]] = \operatorname{E}[n_i p_i] = \operatorname{E}[n_i]\operatorname{E}[p_i] = \mu_n \mu_p;$$ thus $$\operatorname{E}[S] = m\mu_n \mu_p.$$ This assumes that $n_i$ and $p_i$ are independent for each $i$ (from which it follows that each $X_i$ is independent). The variance calculation is done in a similar fashion; $$\operatorname{Var}[S] \overset{\text{ind}}{=} \sum_{i=1}^m \operatorname{Var}[X_i],$$ whence by the law of total variance
$$\begin{align*}
\operatorname{Var}[X_i]
&= \operatorname{Var}[\operatorname{E}[X_i \mid (n_i \cap p_i)]] + \operatorname{E}[\operatorname{Var}[X_i \mid (n_i \cap p_i)]] \\
&= \operatorname{Var}[n_i p_i] + \operatorname{E}[n_i p_i (1-p_i)] \\
&= (\sigma_n^2 \sigma_p^2 + \sigma_n^2 \mu_p^2 + \sigma_p^2 \mu_n^2) + \mu_n \operatorname{E}[p_i(1-p_i)] \\
&= (\sigma_n^2 \sigma_p^2 + \sigma_n^2 \mu_p^2 + \sigma_p^2 \mu_n^2) + \mu_n (\mu_p - (\sigma_p^2 + \mu_p^2)).
\end{align*}$$
To understand the variance of $n_i p_i$, note that for two independent random variables $A$, $B$, with means and standard deviations $\mu_A, \sigma_A, \mu_B, \sigma_B$, respectively,
$$\begin{align*}\operatorname{Var}[AB]
&= \operatorname{E}[(AB)^2] - \operatorname{E}[AB]^2 \\
&= \operatorname{E}[A^2 B^2] - \operatorname{E}[A]^2 \operatorname{E}[B]^2 \\
&= \operatorname{E}[A^2]\operatorname{E}[B^2] - \mu_A^2 \mu_B^2 \\
&= (\operatorname{Var}[A] + \operatorname{E}[A]^2)(\operatorname{Var}[B] + \operatorname{E}[B]^2) - \mu_A^2 \mu_B^2 \\
&= (\sigma_A^2 + \mu_A^2)(\sigma_B^2 + \mu_B^2) - \mu_A^2 \mu_B^2 \\
&= \sigma_A^2 \sigma_B^2 + \sigma_A^2 \mu_B^2 + \sigma_B^2 \mu_A^2. \end{align*}$$
Note that my computation of the variance differs from yours. I have substantiated my results by simulating $m = 10^6$ observations from $X_i$ where $n_i \sim \operatorname{Poisson}(\lambda)$ and $p_i \sim \operatorname{Beta}(a,b)$, for $\lambda = 11$ and $(a,b) = (7,3)$. This should result in $\operatorname{Var}[X_i] = 1001/100$; your results do not match. I should also point out that the reason that your computation does not work is because the total variance of each $X_i$ is not merely due to the expectation of the conditional variance of $X_i$ given $n_i$ and $p_i$; the other term in the law of total variance must also be included, which captures the variability of the conditional expectation of $X_i$. In other words, there is variation in $X_i$ coming from the binomial variance even when $n_i$ and $p_i$ are fixed, but there is also additional variation in the location of $X_i$ arising from the fact that $n_i$ and $p_i$ are not fixed.
Best Answer
TLDR: The sum of two $n$-sided dice is not binomially distributed.
A discrete random variable, $X$, has a binomial distribution, $X\sim Bin(n,p)$ when $Pr(X=x) = \begin{cases}\binom{n}{x}p^x(1-p)^{n-x}&\text{for}~x\in\{0,1,2,\dots,n\}\\ 0 & \text{otherwise}\end{cases}$
For $X$ the sum of two $n$-sided dice however, $Pr(X=x) = \begin{cases} \frac{n - |x-(n+1)|}{n^2} & \text{for}~x\in\{2,3,\dots,2n\}\\ 0 & \text{otherwise}\end{cases}$
Notice that since $n$ will be a fixed number, $Pr(X=x)$ is linear on the intervals $[2,n+1]$ and again linear on the intervals $[n+1,2n]$. This is in direct contrast to the binomial distribution scenario where $Pr(X=x)$ is definitely not linear (as it has terms like $p^x$ and $\binom{n}{x}$ appearing in the formula).
As mentioned in the comments above, as $n$ grows large, the histogram for the sum of two $n$-sided dice approaches the shape of a triangle.
This becomes even more apparent as $n$ gets even larger. Here is the start of the histogram for $n\approx 30$ (its a lot of effort to complete, but you get the idea).
On the other hand, the binomial distribution appears with the all-familiar "bell-shaped" curve.
As such, these are two very different distributions and should not be confused.