kill(all);

vista la specifica natura delle condizioni al contorno,
considero solo i termini elementari della soluz. di 
Michell t.c.
- il campo di spostamenti sia continuo a t=0,t=2%pi
- siano costanti in t
- prevedano stato tensionale limitato (= non infinito) 
  a remoto;

philist : [ 
    r^2 , 
    log(r) ,
    t
];

definisco le componenti di spostamento associate
alla soluzione parametrica che vado a cercare

twomu_ur_list : [ 
    (kappa-1)*r , 
    -1/r ,
    0
];
twomu_ut_list : [ 
    0 , 
    0 ,
    -1/r
];

creo una combinazione lineare dei termini elementari
con coeff. incogniti

m : length(philist);

unks : makelist( A[i], i, 1 , m );

myphi : sum( unks[i]*philist[i] , i , 1 , m  );

mysrr :   1/r  *diff( myphi , r , 1) 
        + 1/r^2*diff( myphi , t , 2 );
mystt : diff( myphi , r , 2 );
mysrt :   1/r^2*diff( myphi, t,1) 
        - 1/r  *diff( myphi,r,1,t,1);

mysrr : expand(mysrr);

mystt : expand(mystt);

mysrt : expand(mysrt);

my_twomu_ur : sum( unks[i] * twomu_ur_list[i] ,i,1,m );

my_twomu_ut : sum( unks[i] * twomu_ut_list[i] ,i,1,m );

definisco le equazioni associate alle condizioni al
contorno, in forma di residuo da annullare

a=b -> a-b [=0]

ipotesi di completa aderenza, interferenza imposta

res_cc_bi :[
    ev(my_twomu_ur , r=1 ) - eta*2*mu,
    ev(my_twomu_ut , r=1 ) - 0
];

ipotesi di assenza di attrito, interferenza imposta 
res_cc_bi :[
    ev(my_twomu_ur , r=1 ) - eta*2*mu,
    ev(mysrt , r=1 ) - 0
];

ipotesi di bordo pressurizzato al foro
res_cc_bi :[
    ev(mysrr , r=1 ) + p,
    ev(mysrt , r=1 ) - 0
];

ipotesi di aderenza, ma a pressione (e non spostamento) radiale imposti
res_cc_bi :[
    ev(mysrr , r=1 ) + p,
    ev(my_twomu_ut , r=1 ) - 0
];

scelgo una delle precedenti

res_cc_bi :[
    ev(my_twomu_ur , r=1 ) - eta*2*mu,
    ev(my_twomu_ut , r=1 ) - 0
];

condizioni a remoto

res_cc_ff : [
    limit(  mysrr , r , infinity) - 0,
    limit(  mystt , r , infinity) - 0,
    limit(  mysrt , r , infinity) - 0
];

creo un'unica lista di residui

res_cc : append( res_cc_bi , res_cc_ff );

risolvo

linsolve(res_cc, unks), globalsolve=true;

controllo che la soluzione annullante i coeff. di fourier
del residuo selezionati anullino pure il residuo stesso
per ogni t.

fullratsimp( ev(res_cc) );

mysrr : ev(mysrr);

mystt : ev(mystt);

mysrt : ev(mysrt);

ricavo le tensioni al bordo interno

[mysrr_bi,mystt_bi,mysrt_bi]:ev( [ mysrr , mystt , mysrt ] , r=1 );

definisco un valore di riferimento per la costante di Kolosov

defkappa:kappa=(3-0.3)/(1+0.3);

interferenza radiale necessaria per ottenere una data pressione di contatto
al bordo interno

linsolve(mysrr_bi=-p,eta);