
(*                     uxy                     *)

(*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])

\[Sigma]xx = D[\[Phi], y, y]
\[Sigma]xy = -D[\[Phi], x, y]
\[Sigma]yy = D[\[Phi], x, x]

exx = ((\[Kappa] + 1)*\[Sigma]xx - (3 - \[Kappa])*\[Sigma]yy)/4
eyy = ((\[Kappa] + 1)*\[Sigma]yy - (3 - \[Kappa])*\[Sigma]xx)/4
exy = \[Sigma]xy
tmux = Integrate[exx, x] + f[y]
tmuy = Integrate[eyy, y] + g[x]

h = exy - (1/2)*(D[tmuy, x] + D[tmux, y])

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


DSolve[h2 == -C[1], f[y], y, GeneratedParameters -> A]
ux = tmux/(2*\[Mu]) /. %

DSolve[h1 == C[1], g[x], x, GeneratedParameters -> B]
uy = tmuy/(2*\[Mu]) /. %

