This answer is nothing more than a variation of Sklivvz's answer. I simply wish to discuss some quantitative ideas following Sklivvz's answer and discuss what I understand (from an aerospace engineering friend) to be a common conceptual mistake - that the application of "mere surface effects" and "application of Bernoulli's principle" is wrong. These "mere surface effects and Bernoulli's principle" follow from Sklivvz's idea as I hope to make clear. Everything in aeroplane physics begins and ends with "aeroplanes thrust air downwards, so the air thrusts aeroplanes up". This answer is written to be understandable to someone like me who knows nothing about fluid dynamics - aside from:
The mathematically elegant and thoroughly enjoyable 2D problems tackled with complex variable theory (see Finding Stagnation Points from the complex potential );
I know that there is a Clay Mathematics prize up for grabs for anyone who can prove the existence of, or give counterexample against, the existence of smooth, globally well-defined solutions to the Navier-Stokes equations;
Aerospace engineering colleagues and friends tell me experimental proof is still Queen in this field: most real fluid dynamics involving aeroplane flight leans on heavily phenomenological models tuned by experiment.
I shall answer by taking these points up in turn.
#Experiment is Queen
From a particular experimental point of view, there is no mystery as to why aeroplanes fly. Rather, the better question, in my opinion, is: "how do they control the inevitable huge lift forces on them to make the latter stably lift in a constant, vertical direction?"
This experimental view is as follows: think of the Beaufort Scale and other scales used by meteorologists to get across the practical meaning of their wind and other warnings: for example, the Fujita Scale for tornadoes and Tropical Cyclone Category Systems, which describe in practical terms the effects of storms of various intensities.
Now I understand flying regulations forbid commercial jetliners from flying at slower than $300 \, \mathrm{km\,h^{-1}}$ before their final approach to the runway. Think of about $300 \, \mathrm{km\,h^{-1}}$ airspeed in terms of the scales I have just spoken of: this is an F4 tornado, category 5 cyclone and is well off the 12-class Beaufort scale. Buildings and structures of any shape the size and weight of fully laden aeroplanes are torn up and borne into the sky or utterly torn down and destroyed. There is NO shortage of lift from a $300 \, \mathrm{km\,h^{-1}}$ relative airspeed to hold almost anything the size and weight of a fully laden commercial jetliner up: at these airspeeds, almost anything of this size and weight and lighter flies. At least it does so fleetingly: if it isn't designed like an aeroplane, as it moves, its attitude changes, and so does the direction of ram pressure: it is then likely to be flipped over and dashed catastrophically onto the ground. Put simply, almost everything flies at this airspeed, but only very special things do so stably.
#Simple Mathematical Models
We can do a back-of-the-envelope estimation of ram pressure in this case: see my drawing below of a simple aerofoil with a significant angle of attack being held stationary in a wind tunnel. I'm going to put some numbers to Sklivvz's description:
Let's suppose the airflow is deflected through some angle $\theta$ radians to model an aeroplane's attitude (not altitude!) on its last approach to landing or as it takes off, flying at $300 \, \mathrm{km\,h^{-1}}$ airspeed or roughly $80 \, \mathrm{m\,s^{-1}}$. I have drawn it with a steep angle of attack. Air near sea-level atmospheric pressure has a density of about $1.25 \, \mathrm{kg\,m^{-3}}$ (molar volume of $0.0224 \, \mathrm{m^{-3}})$. The change in the momentum diagram is shown, whence the change in vertical and horizontal momentum components are (assuming the speed of flow stays roughly constant):
$$\Delta p_v = p_b \sin\theta;\quad\quad\Delta p_h = p_b \,(1-\cos\theta).$$
At the same time, the deflecting wing presents an effective blocking area to the fluid of $\alpha\,A\,\sin\theta$ where $A$ is the wing's actual area and $\alpha$ a scale factor to account for the fact that in the steady state not only fluid right next to the wing is disturbed so that the wing's effective area will be bigger than its actual area. Therefore, the mass of air deflected each second is $\rho\,\alpha\,A\,v\,\sin\theta$ and the lift $L$ and drag $D$ (which force the engines must afford on takeoff) must be:
$$L = \rho\,\alpha\,A\,v^2\,(\sin\theta)^2;\quad\quad D = \rho\,\alpha\,A\,v^2\,(1-\cos\theta)\, \sin\theta.$$
If we plug in an angle of attack of $30$ degrees, assume $\alpha = 1$ and use $A = 1000 \, \mathrm{m^3}$ (roughly the figure for an Airbus A380 wing area), we get a lifting force $L$ for $\rho = 1.25 \, \mathrm{kg\,m^{-3}}$ and $v = 80 \, \mathrm{m\,s^{-1}}$ of $200$ tonne weight. This is rather less than the takeoff weight of a fully laden A380 Airbus (which is $592$ tonnes, according to the A380 Wikipedia page), but it is an astonishingly high weight just the same and within the right order of magnitude. As I said, experiment is Queen here. We see that the wing's effective vertical cross-section is bigger than the actual wing's by a factor of 2 to 3. This is not surprising at steady state, well below the speed of sound flow: the fluid bunches up, and the disturbance is much bigger than just around the wing's neighbourhood. So, plugging in an $\alpha = 3$ (given the experimental fact that the A380 can lift off at 592 tonnes gross laden weight), we get a drag $D$ of
54 tonne weight (538kN) - about half of the Airbus's full thrust of 1.2MN, so this ties in well with the Airbus's actual specifications, given there must be a comfortable margin to lift the aeroplane out of difficulty when needed.
In these F4 / C5 grade winds (and up to three times faster in normal flight), we see, therefore, there simply isn't any shortage of lift. The aeronautical engineering problem is more about keeping this plentiful lift stably directed upwards and allowing the aeroplane to hold a steady attitude and keep any torques arising from lift nonuniformity from flipping the plane over.
As the aeroplane picks up speed, the ram pressure calculated above is proportional to the square of the airspeed (see my answer to Drag force at high speeds ), so that at full speed the effect more than accounts for the drop in air density and the shallower angle of attack - we cannot make this downwards ram pressure without overcoming the much greater horizontal hindwards component - drag - so that is important to fly with low angle of attack for good fuel efficiency.
#Refining the Mathematical Model
It is important to heed that the above description in terms of the momentum difference between incoming air and the downwash begotten by the wing is exactly the same physics as the "more popular" descriptions given in terms of the Bernoulli equation and the integration of pressure around the wing. This is easy to see: the Navier-Stokes equation(See the Wikipedia page for the derivation of the Navier-Stokes equation ), is a very simple application of nothing more than Newton's second and third laws to infinitesimal volumes of fluid, notwithstanding the lack of knowledge about its fundamental mathematical properties (as bespoken by the Clay Mathematics Millenium Prize's unclaimed status: I love the Navier-Stokes equation- such a simple, readily grasped idea so baldly just an embodiment of Newton's laws, yet throwing up profound mysteries that show us, scientists, how little we yet know about the World). The steady-state Navier Stokes equation for a perfect, incompressible fluid is (here $\vec{v}$ is the steady state velocity field and $p$ the scalar pressure field):
$$(\vec{v}\cdot \nabla) \vec{v} = \nabla \left(\frac{|\vec{v}|^2}{2}\right) + \nabla\wedge(\nabla\wedge\vec{v}) = -\nabla p$$
which gives $\nabla\left(p + \frac{|\vec{v}|^2}{2}\right) = 0$ or $p + \frac{|\vec{v}|^2}{2} = \text{const}$ for an irrotational flow ($\nabla\wedge\vec{v} = \vec{0}$) when integrated along the integral curve of $\vec{v}$, i.e. a streamline. Or alternatively, we can argue in a more first principles way in this simple case: the force on an infinitesimal volume is $-\nabla p$, and the acceleration of a particle on the streamline is, by application of the Serret-Frenet formulas (here $s$ is the arc length along the streamline through the particle and $\kappa$ the path's curvature):
$$\mathrm{d}_t (v \hat{\mathbf{t}}) = \mathrm{d}_s v \times \mathrm{d}_t s\, \hat{\mathbf{t}} + v\,\mathrm{d}_s(\hat{\mathbf{t}})\,\mathrm{d}_t s=v\,\mathrm{d}_s v, \hat{\mathbf{t}} - \kappa\,v^2\,\hat{\mathbf{n}}=\mathrm{d}_s \left(\frac{v^2}{2}\right)\, \hat{\mathbf{t}} - \kappa\,v^2\,\hat{\mathbf{n}}$$
whence, on applying $\vec{F} = m \vec{a} \Rightarrow -\nabla p \,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z = \rho\,\vec{a}\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z$, we get:
$$-\nabla p = \rho \left(\mathrm{d}_s \left(\frac{v^2}{2}\right)\, \hat{\mathbf{t}} - \kappa\,v^2\,\hat{\mathbf{n}}\right)$$
which again yields $p + \frac{|\vec{v}|^2}{2} = \text{const}$ when integrated along a streamline (here we can see the sideways (normal to streamline) centripetal force $-v^2\,\hat{\mathbf{n}} / R$ given by the wonted $v^2/R$ formula). So, we can (and will, below), for example, apply the Theorem of Blasius to calculate lift, and be assured it is no more than a quantification of Sklivvz's idea that "aeroplanes thrust air downwards, so the air thrusts aeroplanes up". The pressure difference between the upper and lower surface of a wing exists because the wing is pushing the air down, not a separate phenomenon. Often, one hears that the Bernoulli principle applied to wings is wrong: this is not true. There is a fallacy (to be discussed below) as shown by experiment (and, hand-waving, by theory) in the wonted demonstration of lift using Bernoulli's principle, but the idea is basically sound, as it must be from its derivation from the Navier-Stokes equation and Newton's laws shown above.
#A Joukowsky Aerofoil Calculation and Errors in Wonted Application of Bernoulli's Principle to Wings
We look at a 2D calculation of lift by Bernoulli's principle or, equivalently, by application of the Theorem of Blasius. The common misconception here is that airflows split at the wing's leading edge, and two neighbouring particles will reach the wing's lagging edge at the same time so that the upper particles must fare the curved surface at higher speeds and, therefore, the pressure on the upper wing surface is less. Actually, the upper path particles are sped up much more than this explanation implies and reach the wing's lagging edge well before their lower-path-faring neighbours. See this wonderful video from the University of Cambridge, particularly at about 50 seconds in. This fact shows that the circulation $\oint_\Gamma \vec{v}\cdot\mathrm{d}\vec{r}$ around the wing's surface $\Gamma$ is nonzero, a fact which we intuitively expect from simple theory (as shown below) and which is amply confirmed in experiment: see the video, or go to the end of a runway of a large airport on a damp day so that you can let big commercial jetliners fly over you at about 50m height (take your earmuffs). On a damp day, you will see vortices breaking off the swing' outer edges, will see them swirling in the damp air for many seconds in the aeroplane's wake and, if you take your hearing protection off after the aeroplane has passed, you will hear the vortices crackling in the air, sounding a little like waves washing on the beach. This is much more fun than it sounds when your children are badgering you to do such a thing, and, from the sights and sounds, I learnt heaps more from doing it than I thought I would. Even though the following calculation has an air of theoretical soundness and "first principles" to it, it is important to understand that it too is an experimental model: the circulation is forced into our description, motivated by the confirmation of the former's existence by experiment. The Kutta-Joukowski Condition (see Wikipedia page for Kutta Condition) as well as the Wikipedia page for the Kutta-Joukowski Theorem is little more than an ad-hoc experimentally motivated fix: it is simply this. When we model the flow with a Joukowski aerofoil (described below), there is a sharp, lagging edge on the wing. This begets a singularity with unphysical, infinite velocities. However, by postulating and choosing the right circulation in the flow, we can put a stagnation point at the lagging edge, thus cancelling the singularity, regularising our solution and also forcing the experimentally observed condition that there is only ever one stagnation point at the wing's leading edge, never elsewhere.
Another way to look at this experimentally motivated condition is well explained in this answer to the Physics SE question Does a wing in a potential flow have lift?. An irrotational, inviscid, incompressible flow cannot alone lift a wing. We add circulation to "fudge" a compensation for this theoretical lack: viscosity is "nature's way of enforcing the Kutta-Joukowski condition".
So we begin with the complex variable method (see the Wikipedia page for "Potential flow" in the section "Analysis for two-dimensional flow") to study a potential flow, i.e. irrotational ($\nabla \wedge = \vec{0}$) velocity field $\vec{v}$ with a potential $\psi$ such that $\vec{v} = -\nabla \psi$ that is also incompressible (continuity equation $\nabla\cdot \vec{v} = \nabla^2 \psi = 0$). See also the Physics SE questions Finding Stagnation Points from the complex potential ).
The main method here is to use the Joukowski transform:
$$\omega(z,\,s_z,\,s_\omega) = \frac{ s_\omega }{2}\left(\frac{z}{ s_z } + \frac{ s_z }{z}\right)$$
to map the potential flow corresponding to a spinning, offset cylinder (see the NASA page "Lift of a Rotating Cylinder") into the flow around the image of this cylinder under the Joukowsky transform. The truly weird Flettner Aeroplane actually used spinning cylinders rather than wings to fly successfully. The Joukowsky transform maps the circle $|z| = s_z$ onto the real axis between the points $\omega = \pm s_\omega$ in the $\omega$-plane; this section of the real axis between $\omega = \pm s_\omega $ is then the branch cut for the inverse Joukowski transform. The Joukowsky transform is a two to one mapping, and the branches of the inverse Joukowski transform map the whole $\omega$-Riemann sphere (if we define the stereographic projection so that $|z| = s_\omega $ is the $\omega$-Riemann sphere’s equator) separately to the inside and outside of the circle $|z| = s_z$ in the $z$-plane (which outside and inside can be thought of the Northern and Southern hemispheres of the $z$-Riemann sphere, if the stereographic projection is chosen so that the circle $|z| = s_z$ is the $z$-Riemann sphere’s equator). The $\omega$-Riemann surface is made by slitting two copies of the Riemann sphere along the branch cut and stitching the edges together to get a double cover of genus nought for the $\omega$-Riemann sphere. For this problem, I define the branch cut as slightly different from the real axis section between the $\pm s_\omega$. I define it as the path:
$$\operatorname{Im}(\omega) = h \cos\left(\frac{\pi}{2} \operatorname{Re}(\omega)\right)$$
between the two branch points with an adjustable height parameter $h$, for reasons that will become clear.
The radius $r$ of the spinning cylinder radius is chosen so that the cylinder surface passes through the point $z=+s_z$, which is the image of one of the branch points in the $\omega$ plane. This achieves the sharp edge that becomes the lagging edge of our aerofoil.
The complex potential for the spinning cylinder is:
$$\Omega(z) = v \,e^{-i\alpha}\,\left(z- \delta\right) + \frac{r^2 \,v\, e^{+i\alpha }}{z- \delta } + i\,a\,\log\left(z - \delta \right)$$
where $\alpha$ is the angle of attack, $\delta = \delta_r + i\,\delta_i$ is the offset, and $r$ is the radius of the cylinder steeped in a uniform flow which converges to $v$ metres per second along the positive real axis, as $z\to\infty$. The logarithm and dipole terms put a branch point and pole at the cylinder’s centre, so the flow is perfectly valid outside and on the cylinder. $a$ is the circulation. If we let $\phi$ stand for the angular coordinate labelling the cylinder’s edge, there are two stagnation points on the cylinder with angular coordinates $\phi_\pm$ where $\mathrm{d}_z \Omega(z) = 0$, i.e. when:
$$e^{i\,(\phi_\pm - \alpha)} = -i\frac{a}{2\,v\,r}\pm\sqrt{1-\left(\frac{a}{2\,v\,r }\right)^2} = \exp\left(-\arcsin\frac{a}{2\,v\,r }\right).$$
Now, we map this flow to the $\omega$ plane and apply the Theorem of Blasius to the image of the offset circle so as to work out the lift on this image. The image can be plotted with the Mathematica command:
$$\small{\mathrm{P[\delta_r\_, \delta_i\_] := \\
ParametricPlot[\{Re[\omega[\delta_r + i \delta_i +
\sqrt{(1 - \delta_r)^2 + \delta_i^2} Exp[i \theta]],
Im[\omega[\delta_r + i \delta_i +
\sqrt{(1 - \delta_r)^2 + \delta_i^2} Exp[i \theta]]\}, \{\theta, 0, 2 \pi\}]}}$$
and the result is drawn below in the $\omega$-plane for $s_z = s_\omega = 1$, $\delta_r = -0.1$, $\delta_i = 0.3$ (i.e. the spinning circle offset so that its centre is at $-0.1+i\,0.2$ and with a radius $r = \sqrt{(1 - \delta_r)^2 + \delta_i^2}$ so that its image passes through the branch point $\omega = +s_\omega = 1$ in the $\omega$-plane:
Now we come to the crucial Kutta-Joukowski postulate, an experimental "fudge". The sharp edge on the aerofoil above would normally map the flow in the $z$-plane so that there was an unphysical infinite velocity at this sharp point. In practice, it is seen in wind tunnel tests that the streamlines stay tangent to the upper surface and that there is one stagnation point at the wing's leading edge (intuitively, the air "crashes" here) and no other stagnation points on either the top or bottom of the wing. Sometimes, there is a small region of turbulence around the wing's lagging edge (as in the University of Cambridge video) (i.e. the incompressible potential flow model fails here), or the flow peels smoothly off the lagging edge. The way we achieve effects similar to experiment and "renormalise" our solution is to add the right amount of circulation $a$ to the flow so that one of the stagnation points on the spinning cylinder is mapped to the sharp edge (the branch point at $\omega = +s_\omega$) in the $\omega$-plane: the stagnation thus cancels the otherwise unphysical infinite velocities there and "regularises" our solution. With the radius of the cylinder chosen as $r = \sqrt{(1 - \delta_r)^2 + \delta_i^2}$, it can readily be shown from the equation above for the stagnation point positions that the circulation needed is:
$$a = 2 v\,\delta_i \cos\alpha + 2\,v\,(1-\delta_r) \sin\alpha.$$
This, then, is the wholly experimentally motivated Kutta-Joukowski condition. It is motivated by the knowledge that circulation is observed around wings, that there is experimentally only one stagnation point on the leading edge of the wing and the fact that the right amount of circulation can reproduce these experimentally seen results.
When this is done, the Blasius theorem lift calculation done around the transformed Joukowski aerofoil in the $\omega$-plane is:
$$\begin{array}{lcl}D_\ell - i\,L_\ell &=& \frac{i\,\rho}{2}\oint_{\Gamma_\omega} (\mathrm{d}_\omega \Omega)^2 \,\mathrm{d} \omega\\
&=& \frac{i\,\rho}{2}\oint_{\Gamma_z} (\mathrm{d}_z \Omega)^2 \frac{1}{\mathrm{d}_z \omega}\,\mathrm{d} z\\
&=& -\pi\,\rho \Sigma[\,\mathrm{residues\,of\,}\,(\mathrm{d}_z \Omega)^2 \frac{1}{\mathrm{d}_z \omega}\,\mathrm{at\,poles\,within\,}\Gamma]\\
&=& -4\,\pi\,i\,\rho\,a\,v\,e^{-i\,\alpha}\end{array}$$
where $\Gamma_\omega$ is the Joukowski aerofoil and $\Gamma_z$ the transformed aerofoil (i.e. the spinning cylinder). So there is no lift without circulation. It's worth stating again:
An irrotational, inviscid, incompressible flow cannot alone lift a wing. We add circulation to "fudge" a compensation for this theoretical lack: viscosity is "nature's way of enforcing the Kutta-Joukowski condition".
Now we substitute the Kutta-Joukowski condition to get:
$$D_\ell + i\,L_\ell = 8\,\pi\,i\,\rho\,v^2\,\left(\delta_i\,\cos\alpha + (1-\delta_r)\,\sin\alpha\right) \frac{s_z^2}{s_\omega} e^{+i\alpha}.$$
We now need to scale the velocities so that the relative airspeeds are equal in the $\omega$- and $z$-planes.
The above is the force per unit length (in a direction normal to the page) on the wing, and its direction is the direction in the $\omega$-plane. We have:
$$\lim\limits_{\omega\to\infty} \left(\mathrm{d}_\omega \Omega(\omega(z))\right) = \lim\limits_{z\to\infty} \left(\mathrm{d}_z\Omega(\omega(z))\right) \lim\limits_{\omega\to\infty} \left(\mathrm{d}_\omega z\right) = 2 \,e^{-i\alpha} v \frac{s_z}{s_\omega}$$
so we need $s_\omega = 2$ and $s_z = 1$, then $\delta$ will be a dimensionless parameter defining the offset of the $z$-plane cylinder as a fraction of its radius. But now the $\omega$-plane planform width of the wing is 4 units. Moreover, the above calculation yields the force per unit length (normal to the 2D flow). So, we divide the result for $s_\omega = 2$ and $s_z = 1$ by 4 and then scale up by the total wing area to get the total force on the wing. Furthermore, we need to rotate the flow in the sketch below so that the incoming flow is horizontal (i.e. in the direction of the aeroplane's relative air velocity) in the $\omega$-total force on the wing above becomes:
$$D + i\,L = \pi\,i\,\rho\,v^2\,A\,\left(\delta_i\,\cos\alpha + (1-\delta_r)\,\sin\alpha\right)$$
We witness the d'Alembert paradox: the perfect flow cannot model the drag. Now, let's put some numbers in. If we put $\delta = 0$, then the wing is simply the straight branch cut between $\omega = \pm 1$, so we have a version of the calculation I began with but now refined to take account of the full flow pattern. With $\alpha = 0.3$ (a little less than 20 degrees), $\rho = 1.25 \, \mathrm{kg\,m^{-3}}$, $v=80 \, \mathrm{m\,s^{-1}}$ and $A = 850 \, \mathrm{m^2}$, we get $L=643 \, \mathrm{tonne}$, pretty near to the Airbus's fully laden takeoff weight. If we chose the parameters $\delta_i = 0.2$, $\delta_r =-0.1$ to give a wing shape that does not seem too fanciful for a jetliner wing with the lagging edge flaps wound fully out for takeoff and landing (see the plot below) we get about $1200$ tonnes to lift for our $300 \, \mathrm{km\,h^{-1}}$ airspeed. Clearly, this is optimistic, and the over-reckoning arises from the assumption of equal effectiveness of the whole wingspan, whereas the tips will clearly not be well modelled by the 2D flow. Not all of the wings will work as modelled, thus the $A$ in this formula is somewhat less than the planform area. What the flow model does show (see below), however, is that the effective vertical cross-section presented to the incoming air is much greater than the tilted area $A \,\sin\theta$ assumed in the very simple model at the beginning of my answer. At steady state, a considerable cross-section of air both above and below the vertical cross-section is bent downwards and contributes to the effect "aeroplanes thrust air downwards, so the air thrusts aeroplanes up" described in Sklivvz's answer.
Now, to plot the complete transformed flow in the $\omega$-plane, we must use the inverse Joukowski transform. To do this successfully, one must use the right branches of the inverse transform in the right coordinate patches. For Mathematica, which puts the branch cut for the square root function along the negative real axis (the namespace std::sqrt in Microsoft Visual C++ puts it along the positive real axis), we define the following chart functions, which are particular branches of the inverse transform:
$$\zeta_1(\omega) = \frac{s_z}{s_\omega}\left(\omega- i \sqrt{\omega-s_\omega}\,\sqrt{-\left(\omega+s_\omega\right)}\right)$$
$$\zeta_2(\omega) = \frac{s_z}{s_\omega}\left(\omega+ i \sqrt{\omega-s_\omega}\,\sqrt{-\left(\omega+s_\omega\right)}\right)$$
$$\zeta_3(\omega) = \frac{s_z}{s_\omega}\left(\omega- \sqrt{\omega^2-s_\omega^2}\right)$$
$$\zeta_4(\omega) = \frac{s_z}{s_\omega}\left(\omega+ \sqrt{\omega^2-s_\omega^2}\right)$$
and then the following Mathematica commands will plot the full flow:
$$\small{\mathrm{\Omega[z\_,\,\delta\_,\,v\_,\,r\_,\,a\_,\,\alpha\_,\,s\_]:= v\,e^{-i\,\alpha}\left(\frac{z}{s}-\delta\right) + \frac{r^2\,v\,e^{i\,\alpha}}{\frac{z}{s}-\delta} + i\,a\,Log\left[\frac{z}{s}-\delta\right]}}$$
$$\small{\mathrm{G[z\_,\,\delta_r\_,\,\delta_i\_,\,\alpha\_]:=\Omega\left[z,\,\delta_r+i\,\delta_i,\,1,\,\sqrt{(1-\delta_r)^2 + \delta_i^2},2\,\delta_i Cos[\alpha] + 2\,(1-\delta_r)\,Sin[\alpha],\,\alpha,\,1\right]}}$$
$$\small{\mathrm{S[\delta_r\_, \delta_i\_, \alpha\_, h\_, c\_] := \\
Show[ContourPlot[
Im[If[(Abs[x] < 1 ) \wedge (y > 0) \wedge (y < h\, Cos[\pi x/2]),
G[\zeta_1[x + i y], \delta_r, \delta_i, \alpha]],
If[x < 0, G[\zeta_3[x + i y], \delta_r, \delta_i, \alpha]],
G[\zeta_4[x + i y], \delta_r, \delta_i, \alpha]]]]], \{x, -2, 2\}, \{y, -2,
2\}, Contours \to c, MaxRecursion\to 2, PlotPoints \to 300, AspectRatio \to 1],
P[\delta_r, \delta_i, \{Black, Thick\}]]}}$$
where $\mathrm{P}[]$ is the parametric plot command above used to plot the aerofoil. The above use of the branch functions works for $\delta_r < 0$: other branches are needed for correct results when $\delta_r > 0$. The parameter $h$ bends the branch cut so that it bows upwards and stays inside the aerofoil, thus allowing the branches of the inverse Joukowsky transform to plot the mapped cylinder flow properly. Drawn below is the outcome from the command $\mathrm{S[-0.1, 0.2, 0.2, 0.2, 100]}$, i.e., the flow around the wing for an angle of attack of $0.2$ radians, the circle offset parameters of $-0.1 + 0.2\,i$, a bow in the branch cut so that $h=0.2$. Witness the branch cut inside the aerofoil below and also how far from the wing's surface its effect stretches. The effective vertical component of the wing's area that is presented to the flow is clearly much greater than the actual vertical component of the wing's area, so the factor of $2$ to $3$ scaling in the A380 Airbus lift as reckoned by the simple fluid deflexion calculation seems highly plausible and unsurprising.
Lastly, to come the full circle, here is an animation to be found on the web page "Irrotational plane flows of an inviscid fluid" at the University of Genoa's environmental engineering department; see https://www.diam.unige.it/~irro/. The animation shows the progress of fluid particles for the Joukowski aerofoil flow, illustrates the assertion that the flow above the wing traverses the wing much more quickly than the flow underneath and lastly, shows very well the main thesis that "aeroplanes thrust air downwards".
![Joukowsky Aerofoil Animation][4]
