#This program calculates the stresses and displacements for a given
#biharmonic stress function phi in Cartesian coordinates.
#
#phi:=x^3*y;
#
#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);
#
sxx:=diff(phi,y,y);
syy:=diff(phi,x,x);
sxy:=-diff(phi,x,y);
#
exx:=((kappa+1)*sxx-(3-kappa)*syy)/4;
eyy:=((kappa+1)*syy-(3-kappa)*sxx)/4;
exy:=sxy;
#
tmux1:=int(exx,x)+f(y);
tmuy1:=int(eyy,y)+g(x);
#
#
h:=simplify(expand(exy-(1/2)*(diff(tmuy1,x)+diff(tmux1,y))));
#
#The following two lines are a trick to separate the functions of x
#from the functions of y.
#
h1:=int(diff(h,x),x);
h2:=simplify(h-h1);
#
solution:=dsolve({h1=_C1,h2=-_C1},{f,g});
#
#The following solution will contain three arbitrary constants _C1, _C2, _C3
#representing rigid body displacement.
#
#
ux:=simplify(expand(subs(solution,tmux1/(2*mu))));
uy:=simplify(subs(solution,tmuy1/(2*mu)));