% Kreisringplatte, konstante Flaechenlast, am Innenrand eingespannt

clear all
ri = 1000 ; ra = 3000 ; t = 30 ; E = 2.1e5 ; ny = 0.3 ; p0 = 0.004 ;
K  = E*t^3/(12*(1-ny^2)) ;

A = [ 1        log(ri/ra)         ri^2          ri^2*log(ri/ra)    ; 
      0          1/ri             2*ri        ri*(2*log(ri/ra)+1)  ; 
      0       -(1-ny)/ra^2      2*(1+ny)              3+ny         ;
      0             0              0                  4/ra         ] ;

b = - [ ri^4 ; 4*ri^3 ; 4*(3+ny)*ra^2 ; 32*ra ] * p0/(64*K) ;   

x = A \ b ;                        % Loesung des Gleichungssystems
C1 = x(1) ;  C2 = x(2) ;  C3 = x(3) ;  C4 = x(4) ;

r  = ri:(ra-ri)/100:ra ;
w  = p0*r.^4/(64*K) + C1 + C2*log(r/ra) + C3*r.^2 + C4*r.^2.*log(r/ra) ;
mr = -K*((3+ny)*p0*r.^2/(16*K) - C2*(1-ny)./r.^2 + C3*2*(1+ny) + C4*(3+ny+2*(1+ny)*log(r/ra))) ;
qr = -K*(p0*r/(2*K) + C4*4./r) ;

subplot (3,1,1) ; plot (r , w)  , axis ij , grid on , title ('w(r)') , xlim([0 ra]) 
subplot (3,1,2) ; plot (r , mr) , grid on , title ('m_r(r)') , xlim([0 ra])
subplot (3,1,3) ; plot (r , qr) , grid on , title ('q_r(r)') , xlim([0 ra])

disp(['Maximale Durchbiegung: w_max  = ', num2str(max(w))]) ;
disp(['Maximales Biegemoment: mr_max = ', num2str(max(abs(mr)))]) ;
disp(['Maximale Querkraft:    qr_max = ', num2str(max(abs(qr)))]) ;