Design optimization methods are generally categorized into shape optimization and topology optimization. Shape optimization typically involves a small number of design variables, making it effective for identifying global optima or a wide range of Pareto-optimal solutions. However, its limited number of parameters restricts the representation of complex geometries, resulting in reduced design freedom. In contrast, topology optimization offers greater design freedom than shape optimization by optimizing the distribution of solid elements. The main drawback of topology optimization is that such methods often struggle to identify global optima due to the complexity of the design space. To overcome these trade-offs, this study proposes a novel aerodynamic design optimization method that combines high design freedom with the capability of global optimization by leveraging quantum annealing. Quantum annealing can obtain optimal solutions for binary variables, either qi = 0 or 1, or si = -1 or +1, by minimizing cost functions known as the Quadratic Unconstrained Binary Optimization (QUBO) and Ising models. While quantum annealing has recently been applied to topology optimization of truss and continuum structures, the present approach differs by determining the QUBO and Ising model parameters through machine learning. Computational fluid dynamics (CFD) simulations are conducted to generate training samples, and an immersed boundary method is employed to represent solid and fluid regions using binary variables. In a demonstration test, a single-element airfoil optimization problem is solved to maximize the lift-to-drag ratio. The proposed method successfully identifies an airfoil shape with superior aerodynamic performance compared to the best sample in the training dataset.