The calculation of CLmax is difficult because we must deal with a flow that is viscous, compressible, and highly three-dimensional. Generally, one does not use a Navier-Stokes calculation to estimate maximum lift. This is partly because it takes a very long time to generate a grid and then solve the equations. However, it is also is difficult to estimate the effects of stall strips, fences, and vortex generators that are routinely used on wings and are essential to obtaining acceptable stalling characteristics. Thus, the more usual approach is as follows. The distribution of pressure on the wing is computed from a 3-D panel method. The pressure distributions along streamwise strips are used as input to a two dimensional boundary layer calculation in which the onset of separation is predicted. The actual maximum lift coefficient is based on the boundary layer data, sometimes supplemented with 2-D wind tunnel data.
In the absence of even 2D boundary layer computations, a variety of simpler rules are used. One of these, the pressure difference rule, has been applied frequently to aircraft with high lift systems. Experiments show that the difference between the peak Cp and the Cp at the trailing edge at Clmax varies with Mach number and Reynolds number. This has been correlated in a number of proprietary rules, but for turbulent sections at low Mach number, it varies roughly from 7 or 8 at a Reynolds number of 1 million to as much at 13 or 14 at 6 million and above. Viscous corrections are made to the results of an inviscid panel method (including a reduction in effective flap deflection due to boundary layer decambering) and then the pressure difference rule applied to each section along the span.
This method usually works well, but many approximations are made. The process of high lift prediction therefore relies strongly on wind tunnel data. But, since the flow is sensitive to changes in Reynolds number, good 3-D measurements are rare. This is often one of the areas of greatest uncertainty in aircraft design.
At the early stages of design, it is not possible to run even a panel code and 2-D boundary layer analysis. In such cases, one may compute the distribution of wing lift with a vortex lattice method or Weissinger model and compare the distribution of section Cl with the Clmax, estimated from 2-D data. One provides some margin against stalling of the outer panels to account for aileron deflections and spanwise boundary layer flow. When flaps are deflected, sections just outboard of the flap tend to stall early according to this method. In reality, the flow near the flap edge induces effective camber in the adjacent sections and so their maximum lift coefficient is increased. This effect must be included if reliable estimates of CLmax are to be obtained using this "critical section" approach.

Figure 1. Critical Section Method for CLmax Prediction: Compute CL at which most critical 2D section reaches Clmax.
One might be concerned that the use of 2-D maximum lift data is completely inappropriate for computation of wing CLmax because of 3-D viscous effects. This issue was investigated by the N.A.C.A. in Report 1339. A figure from this paper is reproduced below (Figure 2). It indicates that the "clean wing" CLmax is, in fact, rather poorly predicted by the critical section method. However, when wing fences are used to prevent spanwise boundary layer flow, the Clmax is increased dramatically and does follow the 2-D results quite well over the outer wing sections. The inboard Clmax is considerably higher than would be expected by strip theory, but inboard section Clmax values are generally reduced with the use of stall strips or other devices to make them stall before the tips. Thus, the tip Clmax and lift distribution determine what the inboard Clmax must be to obtain good stall behavior.

Figure 2. Effect of fences on the section lift coefficients of a sweptback wing. Sweep = 45° AR = 8.0, taper = .45, NACA 63(1)A-012 section. Data from NACA Rpt. 1339 Note the result that with fences, outer panel section Cl's are nearly their 2-D values.