A numerical scheme is presented that allows us to compute numerically the coupling between large-scale motions and unresolved turbulent scales for thermonuclear burning fronts in carbon-oxygen white dwarfs. In contrast to most previous calculations, the parameters of the turbulence model are determined by physical considerations and from results of laboratory combustion experiments, with the exception of a closure relation that remains largely undetermined but seems to be of great importance. As a first test case the subgrid model is coupled to a code solving the Euler equation directly by means of the piecewise parabolic method (PPM), and the propagation of a nuclear burning front in a Chandrasekhar mass white dwarf is simulated in two dimensions. We find that energy from the well-resolved Rayleigh-Taylor scale is indeed going into turbulent motions, thereby accelerating the flame speed significantly beyond its laminar value. However, with a conservative choice of the parameters of the subgrid model, the turbulent velocity never exceeds a few percent of the sound speed. Consequently, the model has very little in common with observed Type Ia supernovae. Possible ways out of this dilemma are briefly discussed.