Your lack of fit p-value is what you use to detect curvature, but you cannot determine the source if you have > 1 continuous factors. You cannot actually model any quadratic effects with this design. The 2-level factorial designs will only allow you to fit a hyperplane surface to the data. If the center points are significantly biased as a group away from the hyperplane, your lack of fit test will have a small p-value.

As @Peter_Bartell alluded to, you have a huge run budget. Using center points in a factorial design (full or fractional) is typical of screening scenario where you have too many factors to directly try to estimate a response surface, but you want to know if curvature is present. That is __not__ your situation. You could do a decent response surface design for 3 factors with as few as 20 runs.

I would suggest following Peter's recommendation of an I-optimal design using Custom Design and specifying your model with the RSM button. That will give you all main effects, 2-factor interaction, and quadratic terms. If you can accomodate off-face axial values, a central composite design (CCD) has significant advantages over the I-optimal Custom Design for the same number of runs. You will get much higher power for your quadratic effects in most cases. If you're set on 60+ runs, you could even replicate the full CCD a couple times. If you're unfamiliar with optimal designs, creating a CCD is probably more straight-forward as well. Just make sure you are confident the axial values have set points you can actually run.

-- Cameron Willden