This paper presents the numerical analysis modeling of thrust air bearings and its application to optimum design. In the modeling, the boundary fitted coordinates and divergence formulation method are suitably combined for improving the applicability to the complicated bearing geometry, such as spiral grooved and herring bone bearings. The load capacity, friction torque, spring and damping coefficients are obtained as a function of design variables such as groove depth, groove angle, bearing width ratio and groove width ratio. The calculated results are compared with the experimental results for the verification of validity. Then, using the numerical method, the design variables which optimize the objective function defined as the ratio of friction torque to spring coefficient are determined for a wide range of rotational speeds under various constraints by applying the hybrid optimization technique combining the direct search method and SQP. Some numerical examples are given in the graphical form.