
(*                     ur\[Theta]                     *)

(*This program calculates the stresses and displacements for a given \
biharmonic stress function \[Phi] in Cartesian coordinates. To test it, 
  uncomment the next line.*)

\[Phi] = A*r^2*Log[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])

\[Sigma]rr = 
  Simplify[(1/r)*D[\[Phi], r] + (1/r^2)*D[\[Phi], \[Theta], \[Theta]]]
\[Sigma]\[Theta]\[Theta] = D[\[Phi], r, r]
\[Sigma]r\[Theta] = -Simplify[D[D[\[Phi], \[Theta]]/r, r]]

err = Expand[((\[Kappa] + 
                1)*\[Sigma]rr - (3 - \[Kappa])*\[Sigma]\[Theta]\[Theta])/4]
e\[Theta]\[Theta] = 
  Expand[((\[Kappa] + 
                1)*\[Sigma]\[Theta]\[Theta] - (3 - \[Kappa])*\[Sigma]rr)/4]
er\[Theta] = Expand[\[Sigma]r\[Theta]]
tmur = Integrate[err, r]
e1 = Expand[r*e\[Theta]\[Theta] - tmur]
tmu\[Theta] = Expand[Simplify[Integrate[e1, \[Theta]]]]
tmur1 = tmur + D[f[\[Theta]], \[Theta]]
tmu\[Theta]1 = tmu\[Theta] + g[r] - f[\[Theta]]

h = Expand[
    Simplify[r*(er\[Theta] - (1/2)*((1/r)*D[tmur1, \[Theta]] + 
                  D[tmu\[Theta]1, r] - tmu\[Theta]1/r))]]

(*The following two lines are a trick to separate the functions of r from the \
functions of \[Theta].*)
h1 = Integrate[D[h, \[Theta]], \[Theta]]
h2 = Simplify[h - h1]


DSolve[h1 == 0, f[\[Theta]], \[Theta], GeneratedParameters -> A]
f1[\[Theta]] = f[\[Theta]] /. %
tmur2 = tmur + D[f1[\[Theta]], \[Theta]]
tmu\[Theta]2 = tmu\[Theta] + g[r] - f1[\[Theta]]
DSolve[h2 == 0, g[r], r, GeneratedParameters -> B]
u\[Theta]1 = tmu\[Theta]2/(2*\[Mu]) /. %
ur1 = tmur2/(2*\[Mu])

ur = Expand[ur1]
u\[Theta] = Expand[Simplify[u\[Theta]1]]
