Abstract: Wrinkling deformation is a common instability mode for flexible membrane structures. The numerical simulation on this problem is challenging. Based on the continuum and tension field theory (TFT), a complementarity co-rotational finite element method (FEM) for the wrinkling analysis of pneumatic membrane structures is proposed. By using the co-rotational approach, the finite deformation is decomposed into a rigid body motion in the global coordinate system and small strain deformation in the local coordinate system of the element. The tangent stiffness matrix of a spatial 3-node triangular membrane element is derived. It includes three parts: material stiffness, rotational stiffness and balanced projection stiffness matrices, and covers the influence of a follower load on the elemental stiffness. In the elemental local coordinate system, a wrinkling model is constructed based on the constitutive relation of bi-modulus material, which can judge the status of one element, i.e., ‘taut’, ‘wrinkled’ or ‘slack’. Furthermore, the oscillation of internal force during the iterative solution is eliminated by establishing an equivalent linear complementarity problem. The stability of the algorithm is improved. Numerical examples show that the proposed method can accurately predict the displacement, stress and wrinkling region of pneumatic membrane structures. Compared with the existing methods such as ‘quasi-dynamic’ and ‘penalty’ ways, the proposed method does not require additional solving techniques to ensure convergence. It is convenient for engineering applications.