Since you say you're interested in computing the expectation, here's a way to calculate the moments numerically.

Rearranging a little and defining $Y = X - \log c$, the expectation you're after, which I'll call $I$, is

$$I =  c^{-n} \cdot  \mathbb{E} \left[ \frac{ 1 }{ \left( 1 + e^{Y} \right)^{n}} \right].$$

Massaging the term inside the integral still further, this can be rewritten

$$ I = c^{-n} \cdot \mathbb{E} \frac{ e^{-nY/2} }{ \left( e^{Y/2} + e^{-Y/2} \right)^{n} } = c^{-n} \cdot \mathbb{E} \left[ e^{-nY/2} \frac{ 1 }{ \left( 2 \cosh (Y/2) \right)^{n} } \right] \ .$$

The point of doing this is that $1/(2\cosh Y/2)^{n}$ has a nice Fourier transform. Specifically, 
$$\frac{1}{(2 \cosh (Y/2))^{n}} = \int_{-\infty}^{\infty} e^{iYz} \mathscr{F}_{n}(z) \, dz \ , $$
where
$$ \mathscr{F}_{n}(z) = \frac{1}{2 \pi} \frac{ \Gamma(n/2+iz)\Gamma(n/2-iz)}{ \Gamma(n) } $$
are rapidly decaying, even functions of $z$. (For integer $n$, they can be expressed more simply: for example, $\mathscr{F}_{1}(z) = \frac{1}{2} \operatorname{sech} \pi z$.)

So
$$I = c^{-n} \mathbb{E} \int_{-\infty}^{\infty} e^{-n Y/2} e^{iYz} \mathscr{F}_{n}(z) \, dz \ . $$

Writing out the characteristic function of the Normal random variable $Y$ explicitly, and rearranging a bit, this becomes
> $$I = c^{-n/2} \int_{-\infty}^{\infty} \cos \left[ \left( \log c + \frac{1}{2} \sigma^{2} n \right) z \right] e^{ \frac{1}{2} \sigma^{2} \left( \frac{n^{2}}{4} - z^{2} \right) } \mathscr{F}_{n}(z) \, dz \ , $$

which is easy to evaluate numerically. Here's what the moments look like for $\sigma=1$ and $c$ equal to 0.7, 1, and 1.3:

[![enter image description here][1]][1]

This approach is also useful if $X$ and $Y$ are not Normal, so long as you know their characteristic functions.


  [1]: https://i.sstatic.net/DIGLa.png