(*This file gives the solution for the accelerating rotating rectangular bar, as given in Section 7.4.2*)

(*We write here `C' for dot(Ω).*)
σxx1=ρ*C*x*y/2
σxy1=-ρ*C*(x^2-y^2)/4
σyy1=-ρ*C*x*y/2

ϕ=A[1]*x^5*y +A[2]*x^3*y^3+A[3]*x*y^5+A[4]*x^3*y+A[5]*x*y^3+A[6]*x*y

σxx=D[ϕ,y,y]+σxx1
σxy=-D[ϕ,x,y]+σxy1
σyy=D[ϕ,x,x]+σyy1

t1=σyy/.y->b
t2=σxy/.y->b

t3=σxx/.x->a
t4=σxy/.x->a

(*The stress function is order 6 so the stresses are order 4 in x and y. However, only even powers occur in σxx and only odd powers in σxy. .

We calculate the coefficients of each power of x in these expressions as*)

bc[1]=Coefficient[t1,x^3]
bc[2]=Coefficient[t1,x^1]


bc[3]=Coefficient[t2,x^4]
bc[4]=Coefficient[t2,x^2]
bc[5]=t2

(*The biharmonic equation is 4th order so applying it to a 6the order polynomial generates a quadratic function, but in view of the symmetry the only term is a multiplier of (x*y) .*)

biharm=D[ϕ,x,x,x,x]+D[ϕ,y,y,y,y]+2*D[ϕ,x,x,y,y]
b1=Coefficient[biharm,x*y]

(*We also calculate the three force resultants on x=a as*)

Fx=Integrate[t3,{y,-b,b}]
Fy=Simplify[Integrate[t4,{y,-b,b}]]
M=Integrate[t3*y,{y,-b,b}]

(* Fx is zero by symmetry and Fy must take the value required by Newton's law to achieve the acceleration C. The only weak condition to be imposed is therefore M=0*)

Solve[{bc[1]==0,bc[2]==0,bc[3]==0,bc[4]==0,bc[5]==0,b1==0,M==0},{A[1],A[2],A[3],A[4],A[5],A[6]}]

(*We then substitute the solution into the original stress function and stresses and also into the applied force Fy*)

ϕ2=ϕ/.%
σxx2=Simplify[σxx/.%27]
σxy2=Simplify[σxy/.%27]
σyy2=σyy/.%27
Fy2= Simplify[Fy/.%27]