In our previous work, a strong-coupling approach was successfully used to simulate highly flexible wings interacting with surrounding fluid flows in a globally Eulerian framework for both fluid and solid. However there was a strong assumption to simplify the formulation: the wing structure has the same density as the surrounding fluid. Here, using a fast-tracking method to efficiently identify the solid and fluid areas and solving modified momentum equations and pressure Poisson equation with variable (but still incompressible) density, we were able to apply the same idea to solve flapping-wing system with arbitrary density ratio. In the investigation of effects by different density ratio, for the same shear modulus, higher structural density actually leads to a “softer” appearance with larger deformation of the wing, as a result of relatively stronger inertial effect. Meanwhile, the wake structure changes from inverse Kármán vortex street to Kármán vortex street and results in a transition from thrust producing to drag producing. Using Euler-Bernoulli beam theory, we analytically obtained the first three natural frequencies and their corresponding modes for the same wing used in numerical simulation. We found that the basic dynamics can be adequately described by these three modes and the leading frequency. The linear analysis also predicts the appearance of higher modes at larger density ratio in our nonlinear numerical simulations. On the other hand, the structure shows responses at higher frequency, especially for higher modes, as a clear indication of the system’s nonlinearity.