As the amplitude of the unsteady flow oscillation is large or large changes occur in the mean background flow such as limit cycle oscillation, the traditional proper orthogonal decomposition reduced order model based on linearized time or frequency domain small disturbance solvers can not capture the main nonlinear features. A new nonlinear reduced order model based on the dynamically nonlinear flow equation was investigated. The nonlinear second order snapshot equation in the time domain for proper orthogonal decomposition basis construction was obtained from the Taylor series expansion of the flow solver. The NLR 7301 airfoil configuration and Goland+ wing/store aeroelastic model were used to validate the capability and efficiency of the new nonlinear reduced order model. The simulation results indicate that the proposed new reduced order model can capture the limit cycle oscillation of aeroelastic system very well, while the traditional proper orthogonal decomposition reduced order model will lose effectiveness. reduced-order model, limit cycle oscillation, proper orthogonal decomposition, aeroelasticity PACS: 47.52.+j, 47.40.Hg With the development of computational aeroelasticity, the aeroelastic response can be accurately predicted by the high-fidelity physics-based mathematical model such as computational fluid dynamics (CFD) and computational structure dynamics (CSD) couple solver. However, the computation cost is too large for these high-fidelity methods to be applied to the multidisciplinary conception design with many iterations. Unlike high-fidelity couple solver, the reduced order model (ROM) aims to construct a simple mathematical representation model, which can capture the dominating behavior of aeroelastic system and can be conveniently used in conception design, control and data-driven systems [1].Many approaches for constructing linear flow based aeroelastic ROMs have been developed and were shown to produce numerical results which compare well with highfidelity nonlinear solver. Among these approaches, the reduced order model based on proper orthogonal decomposition (POD/ROM)has become the most popular method. For example, the POD/ROM was successfully applied to the CFD-based aeroelastic analysis of airfoil [2], wing [3] and full aircraft [4], especially in flutter prediction [3][4][5]. Most of the currently proposed ROMs are typically constructed based on the linearized time or frequency domain small disturbance solvers, which are dynamically linearized about the nonlinear steady mean background flow.Most aeroelastic phenomena such as flutter and gust response can be dealt with these ROMs based on dynamically linearization methods. But unfortunately, some important strong nonlinear dynamics with large structure deformation can not be simulated by the small disturbance solvers, for example, limit cycle oscillation (LCO) [6,7]. There are seldom successful reports about the traditional time-domain POD/ROM to predict LCO generated by nonlinear aerodynamics. Although the first order frequency-do...