What are D-branes, really?
There is an abstract algebraic formulation of QFT: this says that an $n$-dimensional QFT is a consistent assignment of spaces of states and of maps between them to $n$-dimensional cobordisms.
If one allows cobordisms with boundary here, one speaks of open-closed QFT. A D-brane in this context is the type of data assigned by the QFT to these boundaries.
There is also a geometric aspect to this: many abstractly defined QFTs are imagined to be "sigma-models". They are supposed to be induced by a process called "quantization" from a functional (the "action") on a space of maps $\Sigma \to X$ from a cobordism $\Sigma$ into a smooth manifold $X$ equipped with extra geometric data (such as metric, connections, etc.)
Under this correspondence, one may ask which abstract algebraic properties of the QFT derive from which geometric aspects of these "background structures". One finds that the data that the QFT assigns to boundaries comes from geometric data on $X$ that tends to look like submanifolds with their own geometric data on them (but may be considerably more general than that!). If so, this geometric data on $X$ is called a D-brane of the sigma-model.
There are many instances of this that are understood at the rough level at which quantum field theory was understood in the 20th century. One special case that is by now under fairly complete mathematical control and which serves as a good guide to the general concept of D-branes is what is called "2d rational CFT" .
There is a complete mathematical classification of 2d rational CFTs in their abstract algebraic form: they are given by special symmetric Frobenius algebra objects internal to a modular tensor category of representation of a vertex operator algebra.
Under this classification theorem, the boundary data = D-branes in the algebraic formulation can be proven to be precisely modules over this Frobenius algebra object.
In special nice cases, one understands where these come from geometrically. The notable example is the Wess-Zumino-Witten model, where the target space is a group manifold. Here one finds that in the simplest case the geometric data corresponding to these D-branes are submanifolds given by conjugacy classes, and carrying twisted vector bundles. More generally, though, the D-branes are given by cocycles in the twisted differential K-theory of the group. So the identification "D-brane = submanifold" is too naive, in general. The correct identification is:
geometric D-brane = geometric data on the sigma-model target space that induces boundary data of the corresponding algebraically defined worldvolume QFT.
For more see
http://ncatlab.org/nlab/show/brane
I'm going to attempt a short, partial answer written for pure mathematicians.
The word "brane" in high-energy physics means "submanifold". The word is short for "membrane". More precisely, it means a submanifold of space that moves in time. A $p$-dimensional brane in space is therefore also a $p+1$-dimensional brane in spacetime.
The letter "D" is for "Dirichlet". In classical differential equations, a Dirichlet boundary condition is a condition that some of the function of the equation are held constant at the boundary of the domain. Geometrically, if $M$ is a manifold and $B \subseteq M$ is a "brane", then $B$ is used as a D-brane if a map $f:\Sigma \to M$ is (a) is mentioned in a differential equation, and (b) is required to satisfy the Dirichlet condition $f(\partial \Sigma) \subseteq D$. For example, if $\Sigma$ is a surface and $f$ is meant to satisfy the minimal surface equation, then the condition $f(\partial \Sigma) \subseteq D$ asks for a minimal surface whose boundary is a loop in $D$ (say in some homotopy class).
This picture is now subject to two reinterpretations. First, the quantum reinterpretation. Instead of demanding that $f$ satisfy its equation, we assume some Lagrangian functional $L$, and we study the formal integral of $\exp(iL(f))$ over all choices of $f$. For instance, instead of imposing the minimal surface equation, we can let $L$ be the area functional of $f(\Sigma)$, or the energy functional. This is often a non-rigorous integral (a Feynman path integral), but nonetheless in favorable cases it appears to be a would-be-rigorous integral and it can be studied. Looking at the space of all maps $f:\Sigma \to M$ for a manifold $M$, with or without a Dirichlet boundary condition, is called a "sigma model" in quantum field theory.
The second reinterpretation is that of string theory. Assume that $\Sigma$ is two-dimensional and that the sigma model is a conformal quantum field theory. Instead of just viewing this model as a quantum field theory, it is viewed as the perturbative expansion (as a power series in the genus of $\Sigma$) of some dynamical theory of $M$ itself, and maybe $M$ and some D-branes in $M$. In order for this second reinterpration to be viable, both $M$ and D-branes in $M$, and the other rules of the sigma model, have to satisfy a number of special conditions. (For instance, without extra decorations such as D-branes, $M$ must satisfy the Einstein equation.) These special conditions are sometimes part of the de facto definition of a D-brane.
I haven't read it (and I wouldn't be able to), but... have a look to this!
http://arxiv.org/pdf/1003.1178.pdf