The well-known Lyapunov direct method provides precise and rigorous theoretical backgrounds for stability analysis of nonlinear systems. However, incorporating damping effects into accurate estimates of the domain of attraction has proved difficult, which is pointed as a minor but crucial drawback. This paper presents a new method to derive a group of Lyapunov functions to reflect the damping effects by considering the first integral of motion obtained from the system governing equations when transfer conductances are neglected. A systematical approach of converting some part of damping loss into some appropriate system energy terms is developed. Examples show that the proposed method remarkably improves the estimation of the region of attraction compared with the conventional approach.