Bolted flange joints are widely used in various engineering structures to combine two or more substructures together. However, the understanding of mechanical mechanism of bolted flange joints is still insufficient due to the lack of efficient mathematical models. This paper proposes an effective mathematical model for bolted flange joints, which is then employed to conduct vibration analysis on bolted flange joint multiple-plate structures. During the modeling process, the flange is modeled by the Kirchhoff plate theory, and the bolted joints are modeled by the two-dimensional face constraint in bolted joint affected region (BJAR). The artificial spring technique is adopted to simulate the face constraint. The friction contact model and equivalent linearization technique are adopted to analyze the stick and slip states of bolted joints. By adopting the Chebyshev orthogonal polynomials as admissible functions, governing equations are derived. Then, the Rayleigh-Ritz method is applied to study the free vibration characteristics, and the Newmark-beta method is employed to study the vibration response under base excitation. The experimental studies are performed on the bolted flange joint two-plate structure (BFJTPS) to illustrate the effectiveness of the proposed theoretical model. The results indicate that the proposed mathematical model can well predict the vibration characteristics of multiple-plate structures under both bolt looseness and no-looseness conditions. The soft nonlinearity phenomenon under bolt looseness condition can be successfully revealed by the model.