There are two approaches being discussed here.
(1) Rigid foundation - The foundation is assmed rigid and soil is elastic but does not offer any tensile resistance. When a bending is applied together with an axial load the tensile zone is assumed to have no resistance contribution. The standard method of calculation is the normal stress equation.
stress=axial load/ares+-moment/I-value*neutral distance to stress point.
(2) Elastic foundation - The foundation has a finite stiffness and can bent with the soil. The latter is modeled as springs according to the subgrade reaction in a computer simulation.
I believe akastud could be interested in the first approach although most of the posters offer solutions for the latter.
My post is intended for the first approach only.
One can calculate any geomertical shape of a foundation subjected to a combination of axial load and biaxial bending. In foundation design the biaxial bending is usually computed as a resultant applied at an angle to the othorgonal axes when performed by hand calculation.
The technique is known as the Saville variation and based on the full stress equation which has the product moment of areas Ixy. In stress analysis we normally have a section symmetrical about at least one axis and so the Ixy vanishes but the actual equation is more complicated. This equation is described in the earlier versions of JE Bowles book on "Foundation analysis and design" with an iterative method given. For completion I list it out as
stress=axial load/area +or- (Mxx-Myy*Ixy/Ixx)*x/(Ixx-Ixy*Ixy/Ixx) +or- (Myy-Mxx*Ixy/Iyy)*y/(Iyy-Ixy*Ixy/Iyy)
note : The Ixx, Iyy and Ixy refer to the centroidal axes of the stressed section, i.e. with tensile area already deducted. +or- depending on whether the applied moment is in the same direction of chosen sign convention.
With computer, the full stress equation can be solved automatically by coordinate geometery. Since the stress equation is in the form of stress=A*x+B*y+C, where A, B and C are functions of section properties of the stressed area of the bearing surface (area, Ixx, Iyy, Ixy), moments (Mxx and Myy) and axial load, by equating stress=0 the stress equation becomes the a straight line equation for the neutral axis of the form y=mx+c, where m is the slope and c is the intercept with the y axis.
The algorithm is started simply by any guess of the neutral axis. The section properties are calculated and straight line equation is solved to obtain an improved neutral axis position (slope and intercept). The computation is stopped once the improvement is small enough to satisfy a specified tolerance. The procedure is inherently stable and always converges.
The difficult part of the algorithm is the computation of the section properties and the standard technique used is to arrange the irregular section as a polygon with a series of triangles. The triangular elements'section propertis are summed together globally and if needed transferred back to centroidal axes, based on which the stress equation holds true.
The section properties of a triangular element, with respect to its own centroid and having nodes 1, 2 & 3 are
area=0.5[x1(y2-y3)+x2(y3-y1)+x3(y1-y2)
xbar (mean x) =1/3(x1+x2+x3)
ybar (mean y) =1/3(y1+y2+y3)
Ixx=area/12[y1*y1+y2*y2+y3*y3]
Iyy=area/12[x1*x1+x2*x2+x3*x3]
Ixy=area/12[x1*y1+x2*y2+x3*y3]
When transferring from centroidal to global axes
Ixx(global)=Ixx(centroid)+area*ybar*ybar etc
By noding the perimeter consecutively in anti-clockwise direction and voids in clockwise direction, with imaginery cut introduced as necessary, the above described procedure works satisfactorily for any polygon and suitable for any geometrical shape foundation.
The above procedure forms the basis of elastic design of a nonhomogeneous material (i.e. compression elastic but zero tension). It is widely used in the old days (about 20 years ago!) for the elastic method design for reinforced concrete where the steel reinforcement is transformed into equivalent areas multipied by the modulus ratio.