Muscle wrapping influences its moment arm and alters the dynamic response of the musculoskeletal system. Conventional musculoskeletal multibody models often simplify muscle wrapping as the shortest path problem, while these models cannot consider the three-dimensional muscular geometry. In comparison, the finite element description of skeletal muscles using three-dimensional solid element discretization leads to expensive computational costs. This paper established a novel flexible multibody modeling method of muscle wrapping based on absolute nodal coordinate formulation (ANCF), which can simultaneously deal with its fusiform geometry, distributed inertia, and active contractions. The skeletal muscle is discretized by full-parameterized ANCF beam elements with variable cross section. The three-dimensional shape of the muscle belly can be expressed by an explicit function, and its time-dependent change can be captured via the surface-to-surface contact detection algorithm. The muscle stress contains the effect of Hill-type active-passive contractions, muscular fascia hyperelasticity, and volume constancy based on continuum mechanics description. As a validation, a forward dynamics simulation of the sagittal elbow rotation driven by surface electromyography (sEMG) data was performed and compared with literature data and experimental results. The established numerical method provides an efficient analysis tool to describe muscle-bone force transmissions, presenting an approach for understanding the hidden biomechanics during human locomotion and human-machine interactions.