A numerical analysis of the initially clamped bolt joint subject to the working pressure is presented in the paper. Special, hexahedral 21- and 28-node isoparametric finite elements have been employed to model the contact zone. In this model, one takes into account loading due to the working pressure in the gap between the gasket and the flange arising as an effect of the progressing joint opening, what has not been considered in recent papers. Nonlinear stiffness characteristics of the bolt and the flange with the gasket are developed. Working pressure corresponding to the critical bolt force resulting in the joint leakage (complete opening between the gasket and the flange) is determined. FE computational results are compared with the available experimental results. The numerical results are presented using the authors' own graphical postprocessor.