This might be because of the expression used for lift computation being slightly different from the popular form of the Kutta-Joukowski formula (rho V gamma).
If you look up textbooks like Low Speed Aerodynamics by Katz, J. and Plotkin, A., you'll notice that in panel methods the panel pressure difference is first computed and then the panel lift vector.
For a 2D flat plate/panel of chord dx, inclined at alpha with freestream velocity V, the panel pressure difference is typically computed using the local tangential velocity (V cos(alpha)), giving,
panel pressure difference, dP = rho V cos(alpha) Gamma / dx
The panel normal force then turns out to be
panel normal force, dN = dP dx
and panel lift, dL = dN cos(alpha) = rho V Gamma cos^2(alpha)
This approach results in an additional cos^2(alpha) term along with the linear part. It may be intuitive in some sense to think about the limiting case when alpha=90 deg and the linear theory still produces 'unnatural' lift. On the other hand, having the cos^2(alpha) term will cause the lift to fall to zero. It is a consequence of computing lift from the normal force (from unsteady Bernoulli theorem to be precise).
If you divide your 'nonlinear' curve by cos^2(alpha) and you get back a straight line, it's almost definitely due to this reason. The attached figure shows the discrepancy for a 2d flat plate. You will see differences depending on your wing's aspect ratio and local angle of attack.
Disclaimer: I'm not sure if VSPAero does exactly this, since the force calculation was recently switched over to computing using Kutta Joukowski. Maybe there's something related that causes the cos^2 term to come up. However, I've seen the cos^2(alpha) term show up in a few other inviscid solvers too.