This paper is concerned with integrated design and operation of energy systems that are subject to significant uncertainties. The problem is cast as a two-stage stochastic programming problem, which can be transformed into a large-scale nonconvex mixed-integer nonlinear programming problem (MINLP). The MINLP exhibits a decomposable structure that can be exploited by nonconvex generalized Benders decomposition (NGBD) for efficient global optimization. This paper extends the NGBD method developed by the authors recently, such that the method can handle non-separable functions and integer operational decisions. Both the standard NGBD algorithm and an enhanced one with piecewise convex relaxations are discussed. The advantages of the proposed formulation and solution method are demonstrated through case studies of two industrial energy systems, a natural gas production network and a polygeneration plant. The first example shows that the two-stage stochastic programming formulation can result in better expected economic performance than the deterministic formulation, and that NGBD is more efficient than a state-of-the-art global optimization solver. The second example shows that the integration of piecewise convex relaxations can improve the efficiency of NGBD by at least an order of magnitude.