The integration of multiple energy systems (MES) provides opportunities to explore the flexibility to accommodate more renewable energy. The concept and methodology of energy hub (EH) enable the standardized modeling of district MES. However, current EH formulations introduce nonlinearities into the modeling and present challenges to analysis and optimization. This paper proposes an automatic and linearized modeling method to formulate energy conversion in EHs, which simplifies the optimization of EH operations. On this basis, the flexibility of an EH is analyzed and quantitatively evaluated based on the ranks of the coupling matrices of the EH and its feasible operational region. Finally, an application of the linearized model on the interval optimization model is illustrated to show how it can suppress uncertainties and fluctuations in distributed renewable energy. A case study is used to demonstrate the effectiveness of the proposed model and the rationality of the flexibility analysis by comparing two EHs with different flexibilities.