||In this work, we present a theoretical and computational framework for the 3D fully non-linear dynamics of lipid bilayers. Lipid bilayers are unique soft materials, operating in general in the low Reynolds limit. While their shape is predominantly dominated by curvature elasticity as in a solid shell, their in-plane behavior is that of a largely inextensible viscous fluid. This mechanical duality provides structural stability and adaptability, allowing membranes to build relatively stable structures that can nevertheless undergo dynamic shape transformations. The equilibrium behavior of lipid bilayers can be understood to a large extent with the classical bending model of Helfrich. For that reason, modeling and simulation of lipid bilayers have mainly focused on exercising Helfrich model to investigate equilibrium configurations of lipid bilayers. Beyond the Helfrich model, more general models are required to describe the dynamic transformations that bilayers undergo, which should capture the interfacial dissipative mechanisms that dominate at sub-cellular scales. Seifert and Langer developed a continuum model explicitly accounting for the bilayer architecture and capturing the major energetic driving forces and dissipative drag forces involved in the dynamics of lipid membranes. However, this model was proposed for small deviations of the geometry from a planar state and derived through a procedure based on the interfacial stress tensor, which is difficult to generalize. Here, we systematically derive a nonlinear SL theory for the dynamics of bilayer membranes in 3D. The cornerstone in this theory is Onsager’s variational principle, which provides a unified framework for the dissipative dynamics of soft-matter systems. In addition to the multi-physics and geometric aspects of the theory, the three-dimensional simulation of lipid bilayers requires unconventional numerical methods since the resulting equations (1) involve higher order derivatives of the parametrization, (2) lead to a mixed system of elliptic and hyperbolic partial differential equations and (3) are stiff and difficult to integrate in time. Here, we propose a discretization based on subdivision surfaces. Regarding time integration, we show that the functional being minimized in Onsager’s principle can be time-discretized, yielding variational integrators that inherit qualitative features of the time- continuous system, in particular the existence of a Lyapunov function and therefore a notion of nonlinear stability.