Numerical instabilities in bone remodeling simulations: the advantages of a node-based finite element approach

Long bone structure occurs in two distinct forms. The bone mass near the joint is primarily found in a distributed, porous trabecular structure, while in the diaphyses a tubular cortical structure is formed. It seems likely that these two observed morphologies come about, at least in part, as a mechanical adaptation to the different mechanical demands in the two regions. Mathematical formulations of this dependency have been proposed, thus facilitating numerical simulations of bone adaptation. Recently two types of discontinuities have been observed in these simulations. The first type (near-field) appears in areas near distributed load application and is characterized by a 'checkerboard' pattern of density wherein adjacent remodeled elements alternate between low and high density. The second type of discontinuity (far-field) appears remote from the load application and is characterized by strut or column-like regions of elements which become fully compact bone while adjacent regions are fully resorbed. In fact, the far-field discontinuity is an accurate representation of bone physiology and morphology since it is consistent with the appearance of cortical bone in the diaphysis. On the other hand, the near-field discontinuity, appears in a region where continuous distributions of intermediate apparent densities (trabecular bone) are expected. This finding may cause some to question whether a single continuum formulation of bone remodeling can predict both discontinuous far-field behavior and continuous near-field behavior. We describe a node-based implementation of current continuum bone remodeling theories which eliminates the spurious near-field discontinuities and preserves the anatomically correct far-field discontinuities, thus indicating that a single biological process may be at work in forming and maintaining both far-field and near-field morphologies.