Let Omega be either R-n or an unbounded strongly Lipschitz domain of R-n, and let Phi be a continuous, strictly increasing, subadditive, and positive function on (0, infinity) of upper type 1 and of strictly critical lower type index p(Phi) is an element of (n/(n + 1), 1]. Let L be a divergence form elliptic operator on L-2(Omega) with the Neumann boundary condition, and assume that the heat semigroup generated by L has the Gaussian property (G(infinity)). In this paper, the authors introduce the Orlicz-Hardy space H-Phi,H-L (Omega) via the nontangential maximal function associated with {e(-t root L)}(t >= 0) and establish its equivalent characterization in terms of the Lusin area function associated with {e(-t root L)}(t >= 0). The authors also introduce the "geometrical" Orlicz-Hardy space H-Phi,H-z(Omega) via the classical Orlicz-Hardy space H-Phi (R-n) and prove that the spaces H-Phi,H-L (Omega) cm and H-Phi,H-z (Omega) coincide with equivalent norms, from which characterizations of H-Phi,H-L (Omega), including the vertical and the nontangential maximal function characterizations associated with {e(-tL)}(t >= 0) and the Lusin area function characterization associated with {e(-tL)}(t >= 0), are deduced. All the above results generalize the well-known results of P. Auscher and E. Russ by taking Phi(t) t for all t is an element of (0, infinity).