When neural networks are applied to solve complex engineering problems, the lack of training data can make the predictions of the surrogate inaccurate. Recently, physics-constrained neural networks were introduced to integrate physical models in the data-driven surrogate to improve the training efficiency with limited data. Nevertheless, the model-form and parameter uncertainty associated with the neural networks can still lead to unreliable predictions. In this article, a new physics-constrained Bayesian neural network (PCBNN) framework is proposed to quantify the uncertainty in physics-constrained neural networks. The bias and variance of predictions are considered simultaneously during the PCBNN training process. The variance and Kullback–Leibler divergence of neural network parameters are incorporated in the total loss function. The weights associated with the different losses are adjusted adaptively. The training of PCBNNs is also formulated as solving a minimax problem where the loss function for the worst-case scenario is minimized. The new PCBNN framework is demonstrated with engineering examples of heat transfer and phase transition based on both simulation data and experimental measurements. The results show that the accuracy and precision of predictions can be improved with the variance consideration in the PCBNN.