Combined heat and electricity operation with variable mass flow rates promotes flexibility, economy, and sustainability through synergies between electric power systems (EPSs) and district heating systems (DHSs). Such combined operation presents a highly nonlinear and nonconvex optimization problem, mainly due to the bilinear terms in the heat flow model—that is, the product of the mass flow rate and the nodal temperature. Existing methods, such as nonlinear optimization, generalized Benders decomposition, and convex relaxation, still present challenges in achieving a satisfactory performance in terms of solution quality and computational efficiency. To resolve this problem, we herein first reformulate the district heating network model through an equivalent transformation and variable substitution. The reformulated model has only one set of nonconvex constraints with reduced bilinear terms, and the remaining constraints are linear. Such a reformulation not only ensures optimality, but also accelerates the solving process. To relax the remaining bilinear constraints, we then apply McCormick envelopes and obtain an objective lower bound of the reformulated model. To improve the quality of the McCormick relaxation, we employ a piecewise McCormick technique that partitions the domain of one of the variables of the bilinear terms into several disjoint regions in order to derive strengthened lower and upper bounds of the partitioned variables. We propose a heuristic tightening method to further constrict the strengthened bounds derived from the piecewise McCormick technique and recover a nearby feasible solution. Case studies show that, compared with the interior point method and the method implemented in a global bilinear solver, the proposed tightening McCormick method quickly solves the heat-electricity operation problem with an acceptable feasibility check and optimality.