A three-dimensional constraint-driven numerical dynamic model of a flapping wing structure called the Dynamic Spar Numerical Model (DSNM) is introduced and implemented. The model currently includes a leading edge spar and a diagonal spar, attached to a body by revolute and spherical joints, respectively. The spars consist of a user-specified number of rigid links connected by compliant joints (CJs): spherical joints with distributed masses and three axis nonlinear torsional spring-dampers. The goal of this model is to quickly simulate mechanisms in a test platform to see how their CJ design properties and spatial distribution affect passive shape change and physical performance metrics. The results of this model can be used as a starting point for further refinement in compliant joint design for passive shape change. Previous research leading to and assumptions made for modeling CJ are presented. The constraints are established, followed by the formulation of a state model used in conjunction with a forward time integrator, and finally several example runs. Modeling the CJs as linear springs produces a nearly symmetric rotation angles through the flapping cycle, while bi-linear springs show the wing is able to flex more during upstroke than downstroke. Increasing damping ratio reduces high frequency oscillations during the flapping cycle and the number of cycles required to reach steady state. Coupling the spring stiffnesses allows an angle about one axis to induce an angle about another axis, where the magnitude is proportional to the coupling term. Modeling both the leading edge and diagonal spars show that the diagonal spar changes the kinematics of the leading edge spar verses only considering the leading edge spar, causing much larger axial rotations in the leading edge spar. The kinematics are very sensitive to CJ location, where moving the CJ toward the wing root causes a stronger response, and adding multiple CJs on the leading edge spar with a CJ on the diagonal spar allow the wing to deform with larger magnitude in all directions. Future work includes implementing a performance metric, experimental verification, applying loads to represent ambient and flight conditions, and using the model as an optimization tool for parameter and spatial optimization.