kill(all);

derivo le componenti di tensione a remoto
per la lastra in trazione

define(
    srr_from_xy(sxx,syy,sxy,t),
    sxx*c^2 + syy*s^2+2*sxy*s*c
), [c = cos(t) , s = sin(t)];

define(
    srt_from_xy(sxx,syy,sxy,t),
    sxy*(c^2-s^2)+(syy-sxx)*s*c
), [c = cos(t) , s = sin(t)];

define(
    stt_from_xy(sxx,syy,sxy,t),
    syy*c^2 + sxx*s^2-2*sxy*s*c
), [c = cos(t) , s = sin(t)];

calcolo i valori specifici di tali componenti
dello stato tensionale a remoto

lastra a trazione:
srr_ff : trigreduce( srr_from_xy(1,0,0,t) );
stt_ff : trigreduce( stt_from_xy(1,0,0,t) );
srt_ff : trigreduce( srt_from_xy(1,0,0,t) );

lastra a taglio remoto
srr_ff : trigreduce( srr_from_xy(1,-1,0,t) );
stt_ff : trigreduce( stt_from_xy(1,-1,0,t) );
srt_ff : trigreduce( srt_from_xy(1,-1,0,t) );

srr_ff : trigreduce( srr_from_xy(1,0,0,t) );

stt_ff : trigreduce( stt_from_xy(1,0,0,t) );

srt_ff : trigreduce( srt_from_xy(1,0,0,t) );

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, o modulati in seconda armonica
  in t (ossia prop. a cos(2t) o sin(2t))
- siano simmetrici rispetto a t=0, asse di simmetria
  per dominio e cond. al contorno.
- prevedano stato tensionale limitato (= non infinito) 
  a remoto;

philist : [ 
    r^2 , 
    log(r) ,
    t, 
    r^(-n+2)*cos(n*t), 
    r^n*cos(n*t) , 
    r^(-n)*cos(n*t)
],n=2;

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);

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

twomu_ur_list : [ 
    (kappa-1)*r , 
    -1/r ,
    0, 
    (kappa+n-1)*r^(-n+1)*cos(n*t), 
            -n *r^( n-1)*cos(n*t) , 
             n *r^(-n-1)*cos(n*t)
],n=2;

twomu_ut_list : [ 
    0 , 
    0 ,
    -1/r, 
   -(kappa-n+1)*r^(-n+1)*sin(n*t), 
             n *r^( n-1)*sin(n*t) , 
             n *r^(-n-1)*sin(n*t)
],n=2;

le compongo linearmente sulla base dei coeff. incogniti

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

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

ipotesi di assenza di attrito
res_cc_bi :[
    ev(my_twomu_ur , r=1 ) - 0,
    ev(mysrt , r=1 ) - 0
];

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

scelgo una delle precedenti

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

condizioni a remoto

res_cc_ff : [
    limit(  mysrr , r , infinity) - srr_ff,
    limit(  mystt , r , infinity) - stt_ff,
    limit(  mysrt , r , infinity) - srt_ff
];

creo un'unica lista di residui

res_cc : append( res_cc_bi , res_cc_ff );

estraggo i coeff. dei termini di residuo espressi
in serie di Fourier

res_A0 : integrate( res_cc , t, 0, 2*%pi )/2/%pi;

res_A2 : integrate( res_cc *cos(2*t) ,t,0,2*%pi )/%pi;

res_B2 : integrate( res_cc *sin(2*t) ,t,0,2*%pi )/%pi;

genero la lista di equazioni da passare al linsolve

eqns : fullratsimp(
    append( res_A0,res_A2,res_B2  )
);

risolvo

linsolve(eqns, 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);

plotto l'andamento delle componenti di tensione a bordo foro

wxplot2d(

    ev( 
        [ mysrr_bi , mystt_bi , mysrt_bi ], 
        defkappa
    ),
    [t,0,2*%pi],
    [legend, "srr","stt","srt" ]
);

suppongo sovrapposta una pressione di contatto non nulla, dovuta ad interferenza;
tale pressione p è definita in senso relativo per confronto con la tensione unitaria
applicata a remoto

pressione minima per non avere tensione radiale trattiva

pmin:fullratsimp( ev(mysrr_bi,t=0) );

coefficiente di attrito minimo per garantire aderenza

fmin: fullratsimp( mysrt_bi/(-mysrr_bi+p) );

plotto quest'ultimo a pmin

ev(fmin,p=pmin);
wxplot2d( ev(%,defkappa) , [t,0,2*%pi]);

provo con pmin*(1+epsilon)

[
    ev(fmin,p=pmin*1.01),
    ev(fmin,p=pmin*1.1),
    ev(fmin,p=pmin*1.2)
]$
wxplot2d( ev(%,defkappa) , [t,0,2*%pi]);

[
    ev(fmin,p=pmin*2),
    ev(fmin,p=pmin*3),
    ev(fmin,p=pmin*4),
    ev(fmin,p=pmin*5),
    ev(fmin,p=pmin*6)
]$
wxplot2d( ev(%,defkappa) , [t,0,2*%pi],[legend,"2","3","4","5","6"]);