kill(all);

M:80;

assume(ri>0,re>0);

dim:[ri=1,re=2,nu=3/10,mu=80000];

espansione in serie di fourier dei due carichi concentrati

sinc(x):=limit(sin(t)/t,t,%pi*x);

rappresentazione del delta dirac con o senza lanczos sigma factor

lanczos:true;

if lanczos then 
    coeff(i):= 
if oddp(i) then 
    0 
else 
    if i = 0 then
        - 1/2/%pi
    else
        - 1/%pi * sinc(i/M)
else
    coeff(i):= 
if oddp(i) then 
    0 
else 
    if i = 0 then
        - 1/2/%pi
    else
        - 1/%pi * 1
    ;

Kolosov's constant
kappa : (3-4*nu)      plane strain
kappa : (3-nu)/(1+nu) plane stress

kappa : (3-nu)/(1+nu);

srr(phi) := 1/r * diff(phi,r,1)+1/r^2*diff(phi,t,2)$
stt(phi) :=  diff(phi,r,2)$
srt(phi) := -diff(1/r*diff(phi,t,1),r,1)$

nabla2(phi):=diff(phi,r,2)+1/r*diff(phi,r,1)+1/r^2*diff(phi,t,2)$

nabla4(phi):=(nabla2(nabla2(phi)))$

soluzione di un anello caricato da azioni radiali al bordo esterno

srr_ri:0$
srr_re:0$
srt_ri:0$
srt_re:0$
stt_ri:0$
stt_re:0$
 ur_ri:0$
 ur_re:0$

quotaparte costante

phiterms : [r^2,log(r)];

urterms : [(kappa-1)*r,-1/r]/2/mu;

unks   : makelist (A[i],         i,1,length(phiterms));
myphi  : sum(unks[i]*phiterms[i],i,1,length(phiterms));

controllo

fullratsimp ( nabla4(myphi) );

define(myur (r) , sum(unks[i]* urterms[i],i,1,length(phiterms)));
define(mysrr(r) , srr(myphi));
define(mysrt(r) , srt(myphi));
define(mystt(r) , stt(myphi));

sols:linsolve(
[
    mysrr(ri)=0,
    mysrt(ri)=0,
    mysrr(re)=coeff(0),
    mysrt(re)=0
],unks
);

[srr_ri,srt_ri,stt_ri,ur_ri] : [srr_ri,srt_ri,stt_ri,ur_ri] + fullratsimp(ev([mysrr(ri),mysrt(ri),mystt(ri),myur (ri)],sols,infeval))$
[srr_re,srt_re,stt_re,ur_re] : [srr_re,srt_re,stt_re,ur_re] + fullratsimp(ev([mysrr(re),mysrt(re),mystt(re),myur (re)],sols,infeval))$

aggiungo a termine cumulativo

termini successivi

for i : 2 thru M step 2 do(

    phiterms : [            r^(i+2),            r^(-i+2),   r^(i  ),  r^(-i  )]*cos(i*t),
    urterms  : [(kappa-i-1)*r^(i+1),(kappa+i-1)*r^(-i+1),-i*r^(i-1),i*r^(-i-1)]*cos(i*t)/2/mu, 
    unks   : makelist (A[i],         i,1,length(phiterms)),
    myphi  : sum(unks[i]*phiterms[i],i,1,length(phiterms)),
    define(myur (r) , sum(unks[i]* urterms[i],i,1,length(phiterms))),
    define(mysrr(r) , srr(myphi)),
    define(mysrt(r) , srt(myphi)),
    define(mystt(r) , stt(myphi)),
    sols:linsolve(
        [
            mysrr(ri)=0,
            mysrt(ri)=0,
            mysrr(re)=coeff(i)*cos(i*t),
            mysrt(re)=0
        ],unks
    ),
    [srr_ri,srt_ri,stt_ri,ur_ri] : [srr_ri,srt_ri,stt_ri,ur_ri] + fullratsimp(ev([mysrr(ri),mysrt(ri),mystt(ri),myur (ri)],sols,infeval)),
    [srr_re,srt_re,stt_re,ur_re] : [srr_re,srt_re,stt_re,ur_re] + fullratsimp(ev([mysrr(re),mysrt(re),mystt(re),myur (re)],sols,infeval)),
    disp([i,ev(stt_ri,t=%pi/2,dim,infeval,numer),ev(ur_re,t=0,dim,infeval,numer)])
);

grandezze nulle

srr_ri;
srt_ri;
srt_re;

grandezze non nulle

plotme : ev([srr_re,stt_ri,stt_re],dim)$

plot2d(plotme,[t,-%pi,%pi],[legend,"srr_re","stt_ri","stt_re"]);

plotme : ev([ur_ri,ur_re],dim)$

wxplot2d(plotme,[t,-%pi,%pi],[legend,"ur_ri","ur_re"]);

integrate(srr_re*cos(t)*re,t,-%pi/2,%pi/2)$
%,numer;

integrate(srr_re*cos(t)*re,t,%pi/2,3*%pi/2)$
%,numer;