There are really two distinct questions here - I will wave my hands to answer one, and give a "better" answer for the other.
First - "real" drag is very, very messy. Flow of liquid around objects (even mathematically simple spheres) changes depending on the velocity of the flow, density and viscosity of the liquid, and the "characteristic length" of the object. This is captured in the dimensionless Reynolds number
$$\rm{Re} = \frac{\rho v L}{\mu}$$
Where $\rho$ is density, $v$ is velocity, $L$ is the characteristic length (diameter, for a sphere) and $\mu$ the viscosity.
When the Reynolds number is very small (microscopic objects, or larger objects at very low velocities) the flow is laminar and can be analyzed mathematically; but as soon as the flow becomes turbulent, the mathematics becomes approximate. A full description of this fills many volumes - and is completely outside the scope of this answer.
For laminar flow, the drag is proportional to velocity - because the drag is viscous, and viscous shear scales with velocity gradient. For a given shape, the distance over which the velocity gradient approaches zero is roughly constant, from which it follows drag is linear with velocity.
When you go to a higher velocity regime, flow becomes turbulent. This is the case for most macroscopic objects (tennis balls etc). When flow is turbulent, you leave "a mass of disturbed air" in the wake of the object. If you consider the object (say a ball) punching a hole in the fluid in front of it, you can see that the mass of air you interact with per second is equal to the projected area of your object $A$ multiplied by the velocity $v$ (to the the length of a tube traversed in one second) and multiplied by the density $\rho$. To bring that mass of air up to a velocity $v$ would require a force equal to $mv$ per unit time; this tells us that drag will be of the form $\rho A v^2$. In fact, the equation is
$$F_d = \frac12 \rho A v^2 C_D$$
Two factors crept in: a factor $\frac12$, and $C_D$. Both these factors capture the fact that not all the fluid in front of the object needs to be accelerated to the velocity of the object: some of it is pushed sideways. The shaped factor $C_D$ captures the geometric parameter - going from a value of 0.5 for a sphere, it is larger for "less aerodynamic" shapes and lower for things that are deliberately shaped to minimize drag. See for example the wiki article with a list of drag coefficients. Note that $C_D$ is really a function of Reynolds number - this further complicates matters.
Usually the drag coefficient is determined experimentally - if you need results that are better than within 50% or so of the true value. But if you need an approximate value, you can look at the table in the linked article, and find "a shape that matches" best. The biggest factor is the projected area - a smaller area will almost always result in lower drag.
That gets us to part 2: calculating the velocity of an object subject to drag. As you surmised, this could be a very hard problem; but for certain situations you can simplify the equations.
One such situation is free fall. It turns out that free fall with quadratic drag is a problem for which there is an exact mathematical solution (assuming $C_D$ is independent of velocity - otherwise you don't have "true" quadratic drag). This solution is given in numerous places - for example here. I give a bit more detail in this earlier answer
So in the case of free fall with quadratic drag, there is a closed form solution if you know the $C_D$ of your object. Unfortunately, the moment the force and velocity are not aligned, the math becomes more complex, and you need to use numerical methods to solve (this is the case, for example, for a projectile fired at an angle, and subject to gravity). A simple consequence of quadratic drag is the fact that a cross wind will increase your drag (see this answer).
There are good (stable) numerical techniques for tackling this kind of problem; alternatively you can get a "quick and dirty" feeling for it by looking at the code I wrote for this earlier answer. As you can see, by alternating the updates of position and velocity you can improve the stability of the solution without needing to take tiny time steps. I'm providing that link for information only - not because it's a good primer in numerical integration of differential equations. But if you google "Runge Kutta" you will find all the material you need...
