Summary
A mixture‐theory‐based model for multi‐constituent solids is presented where each constituent is governed by its own balance laws and constitutive equations. Interactive forces between constituents that emanate from maximization of entropy production inequality provide the coupling between constituent‐specific balance laws and constitutive models. The deformation of multi‐constituent mixtures at the Neumann boundaries requires imposing inter‐constituent coupling constraints such that the constituents deform in a self‐consistent fashion. A set of boundary conditions is presented that accounts for the non‐zero applied tractions, and a variationally consistent method is developed to enforce inter‐constituent constraints at Neumann boundaries in the finite deformation context. The new method finds roots in a local multiscale decomposition of the deformation map at the Neumann boundary. Locally satisfying the Lagrange multiplier field and subsequent modeling of the fine scales via edge bubble functions result in closed‐form expressions for a generalized penalty tensor and a weighted numerical flux that are free from tunable parameters. The key novelty is that the consistently derived constituent coupling parameters evolve with material and geometric nonlinearity, thereby resulting in optimal enforcement of inter‐constituent constraints. Various benchmark problems are presented to validate the method and show its range of application. Copyright © 2016 John Wiley & Sons, Ltd.