$\def\ZZ{\mathbb{Z}}\def\RR{\mathbb{R}}$I'm not sure if this is an answer to the question, since it does refer to $H_{\ast}(M, \ZZ)$, but I think it sheds some interesting light on why the problem is hard.
We begin with a flawed proof attempt. Write $\delta$ for the diagonal map $M \to M \times M$, and $\pi_1$ and $\pi_2$ for the projections from $M \times M$ onto its first and second factor. Let $\rho$ be an integer $k+\ell$-cycle and let $\alpha$ and $\beta$ be a $k$-form and an $\ell$-form on $M$. Then $\int_{\rho} \alpha \wedge \beta = \int_{\delta(\sigma)} \pi_1^{\ast}(\alpha) \wedge \pi_2^{\ast}(\beta)$. Suppose $\rho$ were homologous in $M \times M$ to $\sum \sigma_i \times \tau_i$, for various cycles $\sigma_i$ and $\tau_i$ in $H_{\ast}(M)$, with $\dim \sigma_i + \dim \tau_i=k+\ell$. Then we would have
$$\int_{\rho} \alpha \wedge \beta = \sum \int_{\sigma_i \times \tau_i} \pi_1^{\ast}(\alpha) \wedge \pi_2^{\ast}(\beta) =\sum_{(\dim \sigma_i, \dim \tau_i) = (k, \ell)} \int_{\sigma_i} \alpha \int_{\tau_i} \beta.$$
This would prove the result. (The integrals over terms where $\dim \sigma_i \neq k$ would drop out. If $\dim \sigma_i<k$, then $\pi_1^{\ast}(\alpha)|_{\sigma_i \times \tau_i}=0$ so $\int_{\sigma_i \times \tau_i} \pi_1^{\ast}(\alpha) \wedge \pi_2^{\ast}(\beta)=0$, and likewise if $\dim \tau_i < \ell$.)
Unfortunately, this need not be true. Saying that $\delta(\rho)$ should always be homologous to such a class is saying that $\delta_{\ast} H_m(M)$ should lie in the image of $\bigoplus_{i+j=m} H_i(M, \ZZ) \otimes H_j(M, \ZZ) \to H_{m}(M \times M, \ZZ)$. But the map $\bigoplus_{i+j=m} H_i(M, \ZZ) \otimes H_j(M, \ZZ) \to H_{m}(M \times M, \ZZ)$ need not be surjective -- the cokernel is $\bigoplus_{i+j=m-1} \mathrm{Tor}_1(H_i(M, \ZZ), H_j(M, \ZZ))$ by the Kunneth theorem. There is no reason that $\delta_{\ast}$ needs to land in the image of $\bigoplus H_i \otimes H_j$; I give an example where it doesn't in this Mathoverflow thread,
So we have to work harder. The good news is that the full statement of Kunneth's theorem includes the statement that
$$0 \to \bigoplus_{i+j=m} H_i(M, \ZZ) \otimes H_j(M, \ZZ) \to H_{m}(M \times M, \ZZ) \to \bigoplus_{i+j=m-1} \mathrm{Tor}_1(H_i(M, \ZZ), H_j(M, \ZZ)) \to 0$$
is (noncanonically) split. Choose such a splitting
$$\iota: \bigoplus_{i+j=m-1} \mathrm{Tor}_1(H_i(M, \ZZ), H_j(M, \ZZ)) \to H_m(M, \ZZ).$$
Then we can write $\delta(\rho) = \sum \sigma_i \times \tau_i + \iota(\sum \phi_j)$ for some $\phi_j$ in Tor groups. But the good news is that Tor groups are torsion, so an integer multiple of $\phi_j$ is homologous to zero. This means that $\int_{\iota(\phi_j)} \pi_1^{\ast}(\alpha) \wedge \pi_2^{\ast}(\beta)$ must be zero. So we still have $\int_{\rho} \alpha \wedge \beta = \sum \int_{\sigma_i \times \tau_i} \pi_1^{\ast}(\alpha) \wedge \pi_2^{\ast}(\beta)$ and now the proof concludes as before. $\square$.
What I find interesting about this argument is two things:
(1) We never had to build cup product in $H^{\ast}(M, \ZZ)$, or prove its compatability with de Rham cohomology, we just tried to shove cycles around naively.
(2) Proving Kunneth with the full statement that the sequence is split is not hard. But I have only ever seen it done in a very algebraic way, not by explicitly writing a cycle on $X \times Y$ as $\sum \sigma_i \times \tau_j + \phi$ for some torsion $\phi$. Moreover, there is no canonical choice of splitting. This may explain why no one seems to be able to prove this in a natural way.
The relative de Rham cochain complex $Ω(M,S)$ is simply the homotopy fiber of the pullback map $Ω(M)→Ω(S)$ (computed as the mapping cocone).
Likewise, the relative singular cochain complex $S(M,S)$ is the homotopy fiber of the restriction map $S(M)→S(S)$ (computed as the kernel of a surjective map).
Thus, the natural weak equivalence of functors $Ω→S$ induces a natural weak equivalence (possibly a zigzag) $Ω(M,S)→S(M,S)$ of homotopy fibers, as desired.
Best Answer
After a few days of pondering I think I understood the answer: Let $M$ be a $d$-dimensional manifold as above, and let $i: N\hookrightarrow M$ be the embedding of a submanifold of dimension $\ell$. The cap product between $[N]\in H_\ell (M)$ and a $k$-form $[\omega]\in H^k(M)$ is an element in $\phi\in H_{\ell-k}(M)\cong \big(H^{\ell-k}(M)\big)^*$ (universal coefficient theorem). The element is specified by the formula:
$$ \phi(\zeta)=\int_{N} i^*\omega\wedge i^*\zeta, \quad \zeta\in H^{\ell-k}(M) $$
This is follows from the fact the isomorphism $H_m(M)\cong \big(H^m(M)\big)^*$ is precisely given by $\frown:[C]\mapsto [C]\frown\cdot$ for $[C]\in H_m(M)$, and the cap product satisfies $[C]\frown ([\alpha]\smile[\beta])=([C]\frown[\alpha])\frown[\beta]=([C]\frown[\alpha])([\beta])$ for $[\alpha]\smile[\beta]\in H^m(M)$.
Here we use de Rham (co-)homology classes, but we may also work on a (co)-chain level.