kill(all);

a specific sizing of the shaft

sizing : [ D=36.0, a=54.0, b=78.0 ]$

/home/ebertocchi/_shared/ffcd26/001_shaftdeflection/shaftdeflection_withpeak_image__sAB_is_actually_sAC.png

FIXME     : s_AB in figure is actually s_AC!!
FIXME also: the reaction force due to P at A is actually P*[...]

Let's type bending moment expressions

MfCD : +        F*b/(a+b)*(a+sCD) - F*sCD
                 -C/(a+b)*(a+sCD) + C
       + P*(1-xi)*b/(a+b)*(a+sCD);
        

MfBD : ( F*a/(a+b) + C/(a+b) + P*(a+xi*b)/(a+b) ) * sBD;

MfAC : ( 
   +        F*b/(a+b)
   -          C/(a+b)
   + P*(1-xi)*b/(a+b)
)*sAC;

let's check [dis-]continuity of the typed bending moment expressions;

continuity at D

shouldbezero:
ratsimp(
  ev(MfCD, sCD=   xi *b) 
- ev(MfBD, sBD=(1-xi)*b)
);

wow... now check the consistency of the bending moment jump at C with the locally applied
external torque

shouldbeC:
ratsimp(
- ev(MfAC, sAC=a) 
+ ev(MfCD, sCD=   0)
);

let's evaluate the strain energy

U :   integrate( MfAC^2/E/J/2 , sAC , 0 , a )
    + integrate( MfBD^2/E/J/2 , sBD , 0 , (1-xi)*b )
    + integrate( MfCD^2/E/J/2 , sCD , 0 , (  xi)*b ) $

tentative to simplify U somehow...

U $
ratsimp(%)$
factor(%)$
U:%;

evaluate deflections and rotations by resorting to the Castigliano theorem

f : diff(U,F);

r : diff(U,C);

v : diff(U, P);

now declare the fictitious nature of C and P

ev( [f,r,v] , C=0 , P=0)$
factor(%)$
[f,r,v] : %; 

search for a stationary point of the v(xi) deflection

diff( v , xi ) = 0 ;
expand(%);
sols : solve( % , xi );

Now I have to check whether those solutions are feasible
the D point must be within the B-C segment 
(otherwise the expression of "v" I employed is not correct)
 
0 <= xi <= 1

I try to manipulate the sols expressions in order to point out
whether the previous requirements is satisfied or not

expand(sols);

the second solution is clearly 1 +  a positive quantity, so it's >1
I also can prove that the first solution is >0  if a<b, but that proof is not automatic.
However, I can check anedoctically for specific values of a and b

ev(sols, sizing, numer);

I choose the first solution, so I obtain the xi value which 
identifies the D point at which the v deflection is stationary, and 
in the specific case maximal

xi_vmax : ev(xi, sols[1]);

then I evaluate the deflection there

vmax : factor(ratsimp(ev( v , xi = xi_vmax)));

In order to discuss the behaviour of the structure, I want to plot
the [f,r,vmax] results as a function of the wheel position along the shaft

2d plots are the only readable plots
1 variable/parameter may be monitored in continuous (the abscissa)
other parameters may be monitored in a discrete fashion, if I plot more than one curve
e.g. five different curves for 5 different values of a parameter

in order to reduce the parameters to be discussed, I group some of them in
dimensional groups of the deflection and of the rotation kind

fdim : F*l^3/E/J /* mm*/ $ 

rdim : F*l^2/E/J /*rad*/ $

in order to normalize the results by dividing them with respect to dimensional groups,

I redefine 
a =    eta *l
b = (1-eta)*l

eta = a/(a+b)

[a>0]         [a<b]
  0  <= eta <= 1/2

subst( I , II , III )
I   : is the replacement expression
II  : is the expression to be searched for
III : is the expression the substitution task is to be performed on

assume( l>0 );

I start from the [f,r,vmax] expressions, and I substitute combinations of (eta,l) in place of (a,b)

[f,r,vmax];
subst(    eta *l , a , % )$
subst( (1-eta)*l , b , % )$
factor(ratsimp(%));

/*

now I can conveniently normalize the aforementioned quantities based on the dimensional groups

*/

% / [ fdim ,rdim , fdim ]$
factor(%)$
[ nf , nr , nvmax  ] : %;

now plot the results for discussion

wxplot2d(
    [ nf , nr , nvmax], /* the list of curves (actually expressions) to be plotted */
    [ eta , 0 , 1/2 ],  /* the common abscissa, and its interval*/
    /* further instructions for plotting in the form [keyword , arg1,arg2, arg3, ... ] */
    [ legend , "norm.f" , "norm.r" , "norm.vmax" ] 
);

search for the worst case rotation

diff(nr,eta)=0;
sols:solve( % ,eta);

I doublecheck if any of the solutions is feasible

expand(sols);

the first is feasible (0<=eta<=1/2)

then evaluate the worst case wheel tilting with respect all possible positioning of the wheel
within the simply supported shaft.

ev(nr,sols[1]);
ratsimp(%);
%,numer;
%*rdim;