The energy functional, the governing partial differential equation(s) (PDE), and the boundary conditions need to be consistent with each other in a modeling system. In electrolyte solution study, people usually use a free energy form of an infinite domain system (with vanishing potential boundary condition) and the derived PDE(s) for analysis and computing. However, in many real systems and/or numerical computing, the objective domain is finite, and people still use the similar energy form, PDE(s) but with different boundary conditions, which may cause inconsistency. In this work, (1) we present a mean field free energy functional for electrolyte solution within a finite domain with either physical or numerically required artificial boundary. Apart from the conventional energy components (electrostatic potential energy, ideal gas entropy term and chemical potential term), new boundary interaction terms are added for both Neumann and Dirichlet boundary conditions. These new terms count for physical interactions with the boundary (for real boundary) or the environment influence on the computational domain system (for non-physical but numerically designed boundary). (2) The traditional physical-based Poisson-Boltzmann (PB) equation and Poisson-Nernst-Planck (PNP) equations are proved to be consistent with the new free energy form, and different boundary conditions can be applied. (3) In particular, for inhomogeneous electrolyte with ionic concentration-dependent dielectric permittivity, we derive the generalized Boltzmann distribution (thereby the generalized PB equation) for equilibrium case, and the generalized PNP equations for non-equilibrium case, under different boundary conditions. Numerical tests are performed to demonstrate the different consequences resulted from different energy forms and their derived PDE(s).