You may want to use perturbation theory. This only gives you an approximate answer, but allows for analytic treatment. Your force is considered a small perturbation to the Keplerian elliptic orbit and the resulting equations of motion are expanded in powers of $K$. For linear perturbation theory, only terms linear in $K$ are retained. This simply leads to integrating the perturbation along the unperturbed original orbit. Writing your force as a vector, the perturbing acceleration is
$$
\boldsymbol{a} = K \frac{GM}{r^2c^2}v_r\boldsymbol{v}_t
$$
with $v_r=\boldsymbol{v}{\cdot}\hat{\boldsymbol{r}}$ the radial velocity
($\boldsymbol{v}\equiv\dot{\boldsymbol{r}}$) and
$\boldsymbol{v}_t=(\boldsymbol{v}-\hat{\boldsymbol{r}}(\boldsymbol{v}{\cdot}\hat{\boldsymbol{r}}))$ the rotational component of velocity (the full velocity minus the radial velocity). Here, the dot above denotes a time derivative and a hat the unit vector.
Now, it depends what you mean with 'effect'. Let's work out the changes of the orbital semimajor axis $a$, eccentricity $e$, and direction of periapse.
To summarise the results below: semi-major axis and eccentricity are unchanged, but the direction of periapse rotates in the plane of the orbit at rate
$$
\omega=\Omega \frac{v_c^2}{c^2}
\frac{K}{1-e^2},
$$
where $\Omega$ is the orbital frequency and $v_c=\Omega a$ with $a$ the semi-major axis. Note that (for $K=3$) this agrees with the general relativity (GR) precession rate at order $v_c^2/c^2$ (given by Einstein 1915 but not mentioned in the original question).
change of semimajor axis
From the relation $a=-GM/2E$ (with $E=\frac{1}{2}\boldsymbol{v}^2-GMr^{-1}$ the orbital energy) we have for the change of $a$ due to an external (non-Keplerian) acceleration
$$
\dot{a}=\frac{2a^2}{GM}\boldsymbol{v}{\cdot}\boldsymbol{a}.
$$
Inserting $\boldsymbol{a}$ (note that $\boldsymbol{v}{\cdot}\boldsymbol{v}_t=h^2/r^2$ with angular momentum vector $\boldsymbol{h}\equiv\boldsymbol{r}\wedge\boldsymbol{v}$), we get
$$
\dot{a}=\frac{2a^2Kh^2}{c^2}\frac{v_r}{r^4}.
$$
Since the orbit average $\langle v_r f(r)\rangle=0$ for any function $f$ (see below), $\langle\dot{a}\rangle=0$.
change of eccentricity
From $\boldsymbol{h}^2=(1-e^2)GMa$, we find
$$
e\dot{e}=-\frac{\boldsymbol{h}{\cdot}\dot{\boldsymbol{h}}}{GMa}+\frac{h^2\dot{a}}{2GMa^2}.
$$
We already know that $\langle\dot{a}\rangle=0$, so only need to consider the first term. Thus,
$$
e\dot{e}=-\frac{(\boldsymbol{r}\wedge\boldsymbol{v}){\cdot}(\boldsymbol{r}\wedge\boldsymbol{a})}{GMa}
=-\frac{r^2\;\boldsymbol{v}{\cdot}\boldsymbol{a}}{GMa}
=-\frac{Kh^2}{ac^2}\frac{v_r}{r^2},
$$
where I have used the identity
$(\boldsymbol{a}\wedge\boldsymbol{b}){\cdot}(\boldsymbol{c}\wedge\boldsymbol{d})
=\boldsymbol{a}{\cdot}\boldsymbol{c}\;\boldsymbol{b}{\cdot}\boldsymbol{d}-
\boldsymbol{a}{\cdot}\boldsymbol{d}\;\boldsymbol{b}{\cdot}\boldsymbol{c}$ and the fact $\boldsymbol{r}{\cdot}\boldsymbol{a}_p=0$.
Again $\langle v_r/r^2\rangle=0$ and hence $\langle\dot{e}\rangle=0$.
change of the direction of periapse
The eccentricity vector
$
\boldsymbol{e}\equiv\boldsymbol{v}\wedge\boldsymbol{h}/GM - \hat{\boldsymbol{r}}
$
points (from the centre of gravity) in the direction of periapse, has magnitude $e$, and is conserved under the Keplerian motion (validate all that as an exercise!). From this definition we find its instantaneous change due to external acceleration
$$
\dot{\boldsymbol{e}}=
\frac{\boldsymbol{a}\wedge(\boldsymbol{r}\wedge\boldsymbol{v})
+\boldsymbol{v}\wedge(\boldsymbol{r}\wedge\boldsymbol{a})}{GM}
=\frac{2(\boldsymbol{v}{\cdot}\boldsymbol{a})\boldsymbol{r}
-(\boldsymbol{r}{\cdot}\boldsymbol{v})\boldsymbol{a}}{GM}
=\frac{2K}{c^2}\frac{h^2v_r\boldsymbol{r}}{r^4}
-\frac{K}{c^2}\frac{v_r^2\boldsymbol{v}_t}{r}
$$
where I have used the identity
$\boldsymbol{a}\wedge(\boldsymbol{b}\wedge\boldsymbol{c})=(\boldsymbol{a}{\cdot}\boldsymbol{c})\boldsymbol{b}-(\boldsymbol{a}{\cdot}\boldsymbol{b})\boldsymbol{c}$
and the fact $\boldsymbol{r}{\cdot}\boldsymbol{a}=0$. The orbit averages of these expression are considered in the appendix below. If we finally put everything together, we get
$
\dot{\boldsymbol{e}}=\boldsymbol{\omega}\wedge\boldsymbol{e}
$
with [corrected again]
$$
\boldsymbol{\omega}=\Omega K \frac{v_c^2}{c^2}
(1-e^2)^{-1}\, \hat{\boldsymbol{h}}.
$$
This is a rotation of periapse in the plane of the orbit with angular frequency $\omega=|\boldsymbol{\omega}|$. In particular $\langle e\dot{e}\rangle=\langle\boldsymbol{e}{\cdot}\dot{\boldsymbol{e}}\rangle=0$ in agreement with our previous finding.
Don't forget that due to our usage of first-order perturbation theory these results are only strictly true in the limit $K(v_c/c)^2\to0$. At second-order perturbation theory, however, both $a$ and/or $e$ may change. In your numerical experiments, you should find that the orbit-averaged changes of $a$ and $e$ are either zero or scale stronger than linear with perturbation amplitude $K$.
disclaimer No guarantee that the algebra is correct. Check it!
Appendix: orbit averages
Orbit averages of $v_rf(r)$ with an abitrary (but integrable) function $f(r)$ can be
directly calculated for any type of periodic orbit. Let $F(r)$ be the antiderivative of $f(r)$, i.e. $F'\!=f$, then the orbit average is:
$$
\langle v_r f(r)\rangle = \frac{1}{T}\int_0^T v_r(t)\,f\!\left(r(t)\right) \mathrm{d}t
= \frac{1}{T} \left[F\left(r(t)\right)\right]_0^T = 0
$$
with $T$ the orbital period.
For the orbit averages required in $\langle\dot{\boldsymbol{e}}\rangle$, we must dig a bit deeper. For a Keplerian elliptic orbit
$$
\boldsymbol{r}=a\left((\cos\eta-e)\hat{\boldsymbol{e}}+\sqrt{1-e^2}\sin\eta\,\hat{\boldsymbol{k}}\right)\qquad\text{and}\qquad
r=a(1-e\cos\eta)
$$
with eccentricity vector $\boldsymbol{e}$ and $\hat{\boldsymbol{k}}\equiv\hat{\boldsymbol{h}}\wedge\hat{\boldsymbol{e}}$ a vector perpendicular to $\boldsymbol{e}$ and $\boldsymbol{h}$. Here, $\eta$ is the eccentric anomaly, which is related to the mean anomaly $\ell$ via
$
\ell=\eta-e\sin\eta,
$
such that $\mathrm{d}\ell=(1-e\cos\eta)\mathrm{d}\eta$ and an orbit average becomes
$$
\langle\cdot\rangle = (2\pi)^{-1}\int_0^{2\pi}\cdot\;\mathrm{d}\ell = (2\pi)^{-1}\int_0^{2\pi}\cdot\;(1-e\cos\eta)\mathrm{d}\eta.
$$
Taking the time derivative (note that $\dot{\ell}=\Omega=\sqrt{GM/a^3}$ the orbital frequency) of $\boldsymbol{r}$, we find for the instantaneous (unperturbed) orbital velocity
$$
\boldsymbol{v}=v_c\frac{\sqrt{1-e^2}\cos\eta\,\hat{\boldsymbol{k}}-\sin\eta\,\hat{\boldsymbol{e}}}{1-e\cos\eta}
$$
where I have introduced $v_c\equiv\Omega a=\sqrt{GM/a}$, the speed of the circular orbit with semimajor axis $a$. From this, we find the radial velocity $v_r=\hat{\boldsymbol{r}}{\cdot}\boldsymbol{v}=v_c e\sin\eta(1-e\cos\eta)^{-1}$
and the rotational velocity
$$
\boldsymbol{v}_t = v_c\frac{\sqrt{1-e^2}(\cos\eta-e)\,\hat{\boldsymbol{k}}-(1-e^2)\sin\eta\,\hat{\boldsymbol{e}}}{(1-e\cos\eta)^2}.
$$
With these, we have [corrected again]
$$
\left\langle \frac{h^2v_r\boldsymbol{r}}{r^4}\right\rangle =
\Omega v_c^2\,\hat{\boldsymbol{k}}\,
\frac{e(1-e^2)^{3/2}}{2\pi}\int_0^{2\pi}\frac{\sin^2\!\eta}{(1-e\cos\eta)^4}\mathrm{d}\eta
=\frac{\Omega v_c^2e}{2(1-e^2)}\hat{\boldsymbol{k}}
\\
\left\langle \frac{v_r^2\boldsymbol{v}_t}{r}\right\rangle = \Omega v_c^2\,
\hat{\boldsymbol{k}}\,
\frac{e^2(1-e^2)^{1/2}}{2\pi}\int_0^{2\pi}\frac{\sin^2\!\eta(\cos\eta-e)}{(1-e\cos\eta)^4}\mathrm{d}\eta=0,
$$
in particular, the components in direction $\hat{\boldsymbol{e}}$ average to zero. Thus [corrected again]
$$\left\langle 2\frac{h^2v_r\boldsymbol{r}}{r^4}-\frac{v_r^2\boldsymbol{v}_t}{r}\right\rangle
=\frac{\Omega v_c^2e\,\hat{\boldsymbol{k}}}{(1-e^2)}
$$
