The present study develops a new three-dimensional Timoshenko beam finite element whose length can be varied during transient dynamic analysis. The variable-length element enables the dynamic deployment analysis of flexible appendages with nonnegligible bending stiffness. In addition, the developed scheme employs an implicit time integration whereby energy and momentum in the system are properly conserved, and no artificial numerical dissipation is introduced. The developed beam element is then used in an finite element model of a solar sailcraft, and its deployment dynamics are analyzed allowing for the nonzero bending stiffness of the bundled membranes, as well as the effect of some realistic design imperfections. Nomenclature A = beam cross-sectional area A r = 3 3 matrix defined in Eq. (13) d CM = vector pointing from the geometrical center to the center of mass of a satellite bus E = Young's modulus of beam E r i , e r i = nodal base vector in undeformed and deformed configurations at node r (i 1; 2; 3) f p xE = nodal elastic force for translational degrees of freedom at node p f p xI = nodal inertial force for translational degrees of freedom at node p f r E = nodal elastic force for rotational degrees of freedom at node r f r I = nodal inertial force for rotational degrees of freedom at node r G = shear modulus of beam G i , g i = covariant base vector in undeformed and deformed configurations (i 1; 2; 3) h = beam height I = 3 3 identity matrix I ij = moment-of-inertia components of rigid body I = beam second moment of area i i = local reference unit vector (i 1; 2; 3) J rs = elemental inertia matrix K i , D i = beam stiffness parameters defined in Eq. (35) K pq I , K rs I = stiffness matrix corresponding to inertial force vector at translational and rotational nodes L 0 , L f = initial and final length of cantilevered beam ' = beam length in undeformed configuration M, M pq = elemental mass matrix and its component m = mass of rigid body N= number of beam elements in convergence study N p , N r = shape functions n = unit eigenaxis vector R r I , R r = nodal base matrix in undeformed and deformed configurations at node r T = kinetic energy t = time U in = strain energy w = beam width X, x = position vector in undeformed and deformed configurations X p , x p = nodal position vector in undeformed and deformed configurations at node p _ x p = nodal velocity vector at node p Y i = three-dimensional local coordinate embedded to beam (i 1; 2; 3) r = nodal modified Rodrigues vector at node r ij = components of Green-Lagrange strain t = time-step size x p = incremental nodal position vector at node p r = incremental nodal modified Rodrigues vector at node r " i , i = beam strain parameters defined in Eq. (35) " ij = components of infinitesimal strain = in-plane phase angle of tip mass = beam density ij = components of stress = out-of-plane phase angle of tip mass = pseudo velocity vector for rotational degree of freedom = eigenaxis rotation angle = element domain ! r = nodal angular velocity vector at node r Subsc...