A discrete harmonic surface is a trivalent graph which satisfies the balancing condition in the 3-dimensional Euclidean space and achieves energy minimizing under local deformations. Given a topological trivalent graph, a holomorphic function, and an associated discrete holomorphic quadratic form, a version of the Weierstrass representation formula for discrete harmonic surfaces in the 3-dimensional Euclidean space is proposed. By using the formula, a smooth converging sequence of discrete harmonic surfaces is constructed, and its limit is a classical minimal surface defined with the same holomorphic data. As an application, we have a discrete approximation of the Enneper surface.