#This program calculates the stresses and displacements for a given
#biharmonic stress function phi in polar coordinates.
#
#phi:=A*r^2*ln(r);
#
#This will calculate the displacements for plane stress. For plane strain,
#uncomment kappa:=3-4*nu and comment kappa:=(3-nu)/(1+nu). Alternatively,
#you can comment both lines and get the resuts in terms of kappa.
#
#kappa:=3-4*nu;
#kappa:=(3-nu)/(1+nu);
#
srr:=simplify(diff(phi,r)/r+diff(phi,theta,theta)/r^2);
stt:=simplify(diff(phi,r,r));
srt:=simplify(-diff((diff(phi,theta)/r),r));
#
err:=((kappa+1)*srr-(3-kappa)*stt)/4;
ett:=((kappa+1)*stt-(3-kappa)*srr)/4;
ert:=srt;
#
tmur1:=int(err,r);
#
e1:=r*ett-tmur1;
#
tmut1:=int(e1,theta);
#
tmur2:=tmur1+diff(f(theta),theta);
tmut2:=tmut1+g(r)-f(theta);
#
h:=simplify(expand(r*(ert-(1/2)*((1/r)*diff(tmur2,theta)+diff(tmut2,r)-tmut2/r))));
#
#The following two lines are a trick to separate the functions of r
#from the functions of theta.
#
h1:=int(diff(h,theta),theta);
h2:=h-h1;
#
solution:=dsolve({h1=0,h2=0},{f,g});
#
#The following solution will contain three arbitrary constants _C1, _C2, _C3
#representing rigid body displacement.
#
ur:=simplify(subs(solution,tmur2/(2*mu)));
#ur:=simplify(expand(subs(solution,tmur2/(2*mu))));
ut:=simplify(subs(solution,tmut2/(2*mu)));