Çubuk elemanlar için rijitlik matrisinin çıkarılması oldukça basittir. Aşağıda verdiğim örnekte ise dikdörtgen ince plakların rijitlik matrisinin çıkarılması, yük vektörü ve deplasmanların hazırlanmasına ait kod verilmiştir. 

​

https://www.mathworks.com/matlabcentral/fileexchange/72624-rectangular-thin-plate

​

https://www.youtube.com/watch?v=bFn_Z4c33BQ

​​

%%---Uniformly Loaded Kirchhoff Thin Plate Analysis Using Matrix Displacement Method---%%

%%---Cagatay Tahaoglu---%%

%%---Structural Engineer---%%

%%---ac.tahaoglu@gmail.com---%%

%clear memory%

clear all;

%MATERIAL AND GEOMETRY%

E = input('modulus of elasticity, kN/m^2 : ' );

nu = input('poisson ratio : ');

g = input('unit weight, kN/m^3 : ');

h = input('plate thickness, m : ');

Lx = input('slab length -x direction, m : ');

Ly = input('slab length -y direction, m : ');

m = input('maximum mesh size, m : ');

%LOADS%

DL = g * h;

q = input('load above the plate, kN/m^2 : ');

%SUPPORTS%

top = input('top support (0,1,2)(free, simple, fix)');

bottom = input('bottom support (0,1,2)(free, simple, fix)');

left = input('left support (0,1,2)(free, simple, fix)');

right = input('right support (0,1,2)(free, simple, fix)');

%Flexural Rigidity%

D = E * (h^3) / 12 / (1- (nu*nu));

%Element Stifness Matrix%

a = Lx/(ceil(Lx/m));

b = Ly/(ceil(Ly/m));

p = a/b;

pu = b/a;

F = (4*((p^2)+(pu^2)))+((14-(4*nu))/5);

G = ((2*(p^2))+((1+(4*nu))/5))*b;

H = -((2*pu*pu)+((1+(4*nu))/5))*a;

I = (2*((p*p)-(2*pu*pu)))-((14-(4*nu))/5);

J = ((p^2)-((1+(4*nu))/5))*b;

K = -((2*pu*pu)+((1-nu)/5))*a;

L = (2*((pu*pu)-(2*p*p)))-((14-(4*nu))/5);

M = ((2*p*p)+((1-nu)/5))*b;

N = -((pu*pu)-((1+(4*nu))/5))*a;

O = (-2*((p*p)+(pu*pu)))+((14-(4*nu))/5);

P = ((p*p)-((1-nu)/5))*b;

Q = -((pu*pu)-((1-nu)/5))*a;

R = (4/3)*((p*p)+((1-nu)/5))*b*b;

S = (2/3)*((p*p)-(2*(1-nu)/5))*b*b;

T = (2/3)*((p*p)-((1-nu)/10))*b*b;

U = (1/3)*((p*p)+((1-nu)/5))*b*b;

V = (4/3)*((pu*pu)+((1-nu)/5))*a*a;

W = (2/3)*((pu*pu)-((1-nu)/10))*a*a;

X = (2/3)*((pu*pu)-((2*(1-nu))/5))*a*a;

Y = (1/3)*((pu*pu)+((1-nu)/5))*a*a;

Z = -nu*a*b;

Ke = D/a/b;

kii = [F G H ; G R Z ; H Z V ];

kij = [L M N ; -M T 0 ; N 0 X ];

kik = [O P Q ; -P U 0 ; -Q 0 Y ];

kil = [I J K ; J S 0 ; -K 0 W ];

kji = [L -M N ; M T 0 ; N 0 X ];

kjj = [F -G H ;-G R -Z ; H -Z V ];

kjk = [I -J K ; -J S 0 ; -K 0 W ];

kjl = [O -P Q ; P U 0 ; -Q 0 Y ];

kki = [O -P -Q ; P U 0 ; Q 0 Y ];

kkj = [I -J -K ; -J S 0 ; K 0 W ];

kkk = [F -G -H ; -G R Z ; -H Z V ];

kkl = [L -M -N ; M T 0 ; -N 0 X ];

kli = [I J -K ; J S 0 ; K 0 W ];

klj = [O P -Q ; -P U 0 ; Q 0 Y ];

klk = [L M -N ; -M T 0 ; -N 0 X ];

kll = [F G -H ; G R -Z ; -H -Z V ];

K_element = [kii kij kik kil; kji kjj kjk kjl; kki kkj kkk kkl; kli klj klk kll];

mesh = [(ceil(Ly/m)+1) (ceil(Lx/m)+1)];

joints = flip(reshape((1:(mesh(1)*mesh(2))),[mesh(1) mesh(2)]));

[y, x] = size(joints);

for i = 1:((x-1)*(y-1))

fj = flip(joints);

KL = flip(fj(i-((ceil(i/(y-1))-1)*(y-1)):1+i-((ceil(i/(y-1))-1)*(y-1)) , ceil(i/(y-1)):1+ceil(i/(y-1))));

%%Code Parameters½½

KM(i,:) = [KL(2) KL(1) KL(3) KL(4)];

end

Sii = zeros(x*y);

col = [(KM(:,1)) (KM(:,1))];

for i = 1:((x-1)*(y-1))

Sii(col(i),col(i+((x-1)*(y-1))))=1;

end

SS_ii = repmat(kii,x*y,x*y).*repelem(Sii,3,3);

Sij = zeros(x*y);

col = [(KM(:,1)) (KM(:,2))];

for i = 1:((x-1)*(y-1))

Sij(col(i),col(i+((x-1)*(y-1))))=1;

end

SS_ij = repmat(kij,x*y,x*y).*repelem(Sij,3,3);

Sik = zeros(x*y);

col = [(KM(:,1)) (KM(:,3))];

for i = 1:((x-1)*(y-1))

Sik(col(i),col(i+((x-1)*(y-1))))=1;

end

SS_ik = repmat(kik,x*y,x*y).*repelem(Sik,3,3);

Sil = zeros(x*y);

col = [(KM(:,1)) (KM(:,4))];

for i = 1:((x-1)*(y-1))

Sil(col(i),col(i+((x-1)*(y-1))))=1;

end

SS_il = repmat(kil,x*y,x*y).*repelem(Sil,3,3);

Sji = zeros(x*y);

col = [(KM(:,2)) (KM(:,1))];

for i = 1:((x-1)*(y-1))

Sji(col(i),col(i+((x-1)*(y-1))))=1;

end

SS_ji = repmat(kji,x*y,x*y).*repelem(Sji,3,3);

Sjj = zeros(x*y);

col = [(KM(:,2)) (KM(:,2))];

for i = 1:((x-1)*(y-1))

Sjj(col(i),col(i+((x-1)*(y-1))))=1;

end

SS_jj = repmat(kjj,x*y,x*y).*repelem(Sjj,3,3);

Sjk = zeros(x*y);

col = [(KM(:,2)) (KM(:,3))];

for i = 1:((x-1)*(y-1))

Sjk(col(i),col(i+((x-1)*(y-1))))=1;

end

SS_jk = repmat(kjk,x*y,x*y).*repelem(Sjk,3,3);

Sjl = zeros(x*y);

col = [(KM(:,2)) (KM(:,4))];

for i = 1:((x-1)*(y-1))

Sjl(col(i),col(i+((x-1)*(y-1))))=1;

end

SS_jl = repmat(kjl,x*y,x*y).*repelem(Sjl,3,3);

Ski = zeros(x*y);

col = [(KM(:,3)) (KM(:,1))];

for i = 1:((x-1)*(y-1))

Ski(col(i),col(i+((x-1)*(y-1))))=1;

end

SS_ki = repmat(kki,x*y,x*y).*repelem(Ski,3,3);

Skj = zeros(x*y);

col = [(KM(:,3)) (KM(:,2))];

for i = 1:((x-1)*(y-1))

Skj(col(i),col(i+((x-1)*(y-1))))=1;

end

SS_kj = repmat(kkj,x*y,x*y).*repelem(Skj,3,3);

Skk = zeros(x*y);

col = [(KM(:,3)) (KM(:,3))];

for i = 1:((x-1)*(y-1))

Skk(col(i),col(i+((x-1)*(y-1))))=1;

end

SS_kk = repmat(kkk,x*y,x*y).*repelem(Skk,3,3);

Skl = zeros(x*y);

col = [(KM(:,3)) (KM(:,4))];

for i = 1:((x-1)*(y-1))

Skl(col(i),col(i+((x-1)*(y-1))))=1;

end

SS_kl = repmat(kkl,x*y,x*y).*repelem(Skl,3,3);

Sli = zeros(x*y);

col = [(KM(:,4)) (KM(:,1))];

for i = 1:((x-1)*(y-1))

Sli(col(i),col(i+((x-1)*(y-1))))=1;

end

SS_li = repmat(kli,x*y,x*y).*repelem(Sli,3,3);

Slj = zeros(x*y);

col = [(KM(:,4)) (KM(:,2))];

for i = 1:((x-1)*(y-1))

Slj(col(i),col(i+((x-1)*(y-1))))=1;

end

SS_lj = repmat(klj,x*y,x*y).*repelem(Slj,3,3);

Slk = zeros(x*y);

col = [(KM(:,4)) (KM(:,3))];

for i = 1:((x-1)*(y-1))

Slk(col(i),col(i+((x-1)*(y-1))))=1;

end

SS_lk = repmat(klk,x*y,x*y).*repelem(Slk,3,3);

Sll = zeros(x*y);

col = [(KM(:,4)) (KM(:,4))];

for i = 1:((x-1)*(y-1))

Sll(col(i),col(i+((x-1)*(y-1))))=1;

end

SS_ll = repmat(kll,x*y,x*y).*repelem(Sll,3,3);

%SYSTEM STIFFNESS MATRIX%%

SS = (SS_ii + SS_ij + SS_ik + SS_il + SS_ji + SS_jj + SS_jk + SS_jl + SS_ki + SS_kj + SS_kk + SS_kl + SS_li + SS_lj + SS_lk + SS_ll);

%borders%

border_top = joints(1,:);

border_bottom = joints(end,:);

border_left = joints(:,1)';

border_right = joints(:,end)';

%Edge Supports%

free = zeros(3);

simple = [1000000 0 0; 0 0 0; 0 0 0];

fix = 1000000.*eye(3);

if top == 0

for k = 1:length(border_top)

SS(3*border_top(k)-2:3*border_top(k),3*border_top(k)-2:3*border_top(k))=SS(3*border_top(k)-2:3*border_top(k),3*border_top(k)-2:3*border_top(k))+free;

end

elseif top == 1

for k = 1:length(border_top)

SS(3*border_top(k)-2:3*border_top(k),3*border_top(k)-2:3*border_top(k))=SS(3*border_top(k)-2:3*border_top(k),3*border_top(k)-2:3*border_top(k))+simple;

end

elseif top == 2

for k = 1:length(border_top)

SS(3*border_top(k)-2:3*border_top(k),3*border_top(k)-2:3*border_top(k))=SS(3*border_top(k)-2:3*border_top(k),3*border_top(k)-2:3*border_top(k))+fix;

end

end

if bottom == 0

for k = 1:length(border_bottom)

SS(3*border_bottom(k)-2:3*border_bottom(k),3*border_bottom(k)-2:3*border_bottom(k))=SS(3*border_bottom(k)-2:3*border_bottom(k),3*border_bottom(k)-2:3*border_bottom(k))+free;

end

elseif bottom == 1

for k = 1:length(border_bottom)

SS(3*border_bottom(k)-2:3*border_bottom(k),3*border_bottom(k)-2:3*border_bottom(k))=SS(3*border_bottom(k)-2:3*border_bottom(k),3*border_bottom(k)-2:3*border_bottom(k))+simple;

end

elseif bottom == 2

for k = 1:length(border_bottom)

SS(3*border_bottom(k)-2:3*border_bottom(k),3*border_bottom(k)-2:3*border_bottom(k))=SS(3*border_bottom(k)-2:3*border_bottom(k),3*border_bottom(k)-2:3*border_bottom(k))+fix;

end

end

if left == 0

for k = 1:length(border_left)

SS(3*border_left(k)-2:3*border_left(k),3*border_left(k)-2:3*border_left(k))=SS(3*border_left(k)-2:3*border_left(k),3*border_left(k)-2:3*border_left(k))+free;

end

elseif left == 1

for k = 1:length(border_left)

SS(3*border_left(k)-2:3*border_left(k),3*border_left(k)-2:3*border_left(k))=SS(3*border_left(k)-2:3*border_left(k),3*border_left(k)-2:3*border_left(k))+simple;

end

elseif left == 2

for k = 1:length(border_left)

SS(3*border_left(k)-2:3*border_left(k),3*border_left(k)-2:3*border_left(k))=SS(3*border_left(k)-2:3*border_left(k),3*border_left(k)-2:3*border_left(k))+fix;

end

end

if right == 0

for k = 1:length(border_right)

SS(3*border_right(k)-2:3*border_right(k),3*border_right(k)-2:3*border_right(k))=SS(3*border_right(k)-2:3*border_right(k),3*border_right(k)-2:3*border_right(k))+free;

end

elseif right == 1

for k = 1:length(border_right)

SS(3*border_right(k)-2:3*border_right(k),3*border_right(k)-2:3*border_right(k))=SS(3*border_right(k)-2:3*border_right(k),3*border_right(k)-2:3*border_right(k))+simple;

end

elseif right == 2

for k = 1:length(border_right)

SS(3*border_right(k)-2:3*border_right(k),3*border_right(k)-2:3*border_right(k))=SS(3*border_right(k)-2:3*border_right(k),3*border_right(k)-2:3*border_right(k))+fix;

end

end

SS = Ke*SS;

%%Fixed-end Moments%%

Q_DL = h*g*a*b;

Q_q = q*a*b;

Qab = Q_DL + Q_q;

Mx = Qab*a*a*b/24;

My = Qab*a*b*b/24;

Tz = Qab/4;

LoadElement = [ -Tz -Mx My -Tz Mx My -Tz Mx -My -Tz -Mx -My ]';

P1 = LoadElement(1:3);

P2 = LoadElement(4:6);

P3 = LoadElement(7:9);

P4 = LoadElement(10:12);

L1 = zeros(1,x*y)';

L1(KM(:,1)) = 1;

L2 = zeros(1,x*y)';

L2(KM(:,2)) = 1;

L3 = zeros(1,x*y)';

L3(KM(:,3)) = 1;

L4 = zeros(1,x*y)';

L4(KM(:,4)) = 1;

LOAD1 = repmat(P1,x*y,1).*repelem(L1,3,1);

LOAD2 = repmat(P2,x*y,1).*repelem(L2,3,1);

LOAD3 = repmat(P3,x*y,1).*repelem(L3,3,1);

LOAD4 = repmat(P4,x*y,1).*repelem(L4,3,1);

LOAD = LOAD1 + LOAD2 + LOAD3 + LOAD4;

%%Displacements%%

DELTA = inv(SS)*LOAD;

Dz = flip(reshape(DELTA(1:3:end),[],x));

%Internal Forces%

sig1 = zeros(16,12);

step = 0;

for ksi = [0 0 1 1]

for ita = [0 1 1 0]

step = step + 1;

sig1_1 = ((6/a/a)*(1-ksi-ksi)*(1-ita))+((6*nu/b/b)*(1-ksi)*(1-ita-ita));

sig1_2 = (2*nu/b)*(1-ksi)*(2-ita-ita-ita);

sig1_3 = (-2/a)*(2-ksi-ksi-ksi)*(1-ita);

sig1_4 = ((6/a/a)*(1-ksi-ksi)*ita)-((6*nu/b/b)*(1-ksi)*(1-ita-ita));

sig1_5 = (2*nu/b)*(1-ksi)*(1-ita-ita-ita);

sig1_6 = (-2/a)*(2-ksi-ksi-ksi)*ita;

sig1_7 = (-(6/a/a)*(1-ksi-ksi)*(ita))-((6*nu/b/b)*(ksi)*(1-ita-ita));

sig1_8 = (2*nu/b)*(1-ita-ita-ita)*ksi;

sig1_9 = (-2/a)*(1-ksi-ksi-ksi)*ita;

sig1_10= (-(6/a/a)*(1-ksi-ksi)*(1-ita))+((6*nu/b/b)*(ksi)*(1-ita-ita));

sig1_11= (2*nu/b)*(2-ita-ita-ita)*ksi;

sig1_12= (-2/a)*(1-ksi-ksi-ksi)*(1-ita);

sig1(step,:) = D*[ sig1_1 sig1_2 sig1_3 sig1_4 sig1_5 sig1_6 sig1_7 sig1_8 sig1_9 sig1_10 sig1_11 sig1_12 ];

end

end

sig1_i = sig1(1,:) ; sig1_j = sig1(6,:); sig1_k = sig1(11,:) ; sig1_l = sig1(16,:);

sig2 = zeros(16,12);

step = 0;

for ksi = [0 0 1 1]

for ita = [0 1 1 0]

step = step + 1;

sig2_1 = ((6*nu/a/a)*(1-ksi-ksi)*(1-ita))+((6/b/b)*(1-ksi)*(1-ita-ita));

sig2_2 = (2/b)*(1-ksi)*(2-ita-ita-ita);

sig2_3 = (-2*nu/a)*(2-ksi-ksi-ksi)*(1-ita);

sig2_4 = ((6*nu/a/a)*(1-ksi-ksi)*ita)-((6/b/b)*(1-ksi)*(1-ita-ita));

sig2_5 = (2/b)*(1-ksi)*(1-ita-ita-ita);

sig2_6 = (-2*nu/a)*(2-ksi-ksi-ksi)*ita;

sig2_7 = (-(6*nu/a/a)*(1-ksi-ksi)*(ita))-((6/b/b)*(ksi)*(1-ita-ita));

sig2_8 = (2/b)*(1-ita-ita-ita)*ksi;

sig2_9 = (-2*nu/a)*(1-ksi-ksi-ksi)*ita;

sig2_10= (-(6*nu/a/a)*(1-ksi-ksi)*(1-ita))+((6/b/b)*(ksi)*(1-ita-ita));

sig2_11= (2/b)*(2-ita-ita-ita)*ksi;

sig2_12= (-2*nu/a)*(1-ksi-ksi-ksi)*(1-ita);

sig2(step,:) = D*[ sig2_1 sig2_2 sig2_3 sig2_4 sig2_5 sig2_6 sig2_7 sig2_8 sig2_9 sig2_10 sig2_11 sig2_12 ];

end

end

sig2_i = sig2(1,:) ; sig2_j = sig2(6,:); sig2_k = sig2(11,:) ; sig2_l = sig2(16,:);

sig3 = zeros(16,12);

step = 0;

for ksi = [0 0 1 1]

for ita = [0 1 1 0]

step = step + 1;

sig3_1 = ((1-nu)/a/b)*((1-(6*ksi)*(1-ksi))-((6*ita)*(1-ita)));

sig3_2 = ((1-nu)/a)*(1-ita-ita-ita)*(1-ita);

sig3_3 = ((nu-1)/b)*(1-ksi-ksi-ksi)*(1-ksi);

sig3_4 = ((nu-1)/a/b)*((1-(6*ksi)*(1-ksi))-((6*ita)*(1-ita)));

sig3_5 = ((nu-1)/a)*(2-ita-ita-ita)*ita;

sig3_6 = ((1-nu)/b)*(1-ksi-ksi-ksi)*(1-ksi);

sig3_7 = ((1-nu)/a/b)*((1-(6*ksi)*(1-ksi))-((6*ita)*(1-ita)));

sig3_8 = ((1-nu)/a)*(2-ita-ita-ita)*ita;

sig3_9 = ((nu-1)/b)*(2-ksi-ksi-ksi)*ksi;

sig3_10= ((nu-1)/a/b)*((1-(6*ksi)*(1-ksi))-((6*ita)*(1-ita)));

sig3_11= ((nu-1)/a)*(1-ita-ita-ita)*(1-ita);

sig3_12= ((1-nu)/b)*(2-ksi-ksi-ksi)*ksi;

sig3(step,:) = D*[ sig3_1 sig3_2 sig3_3 sig3_4 sig3_5 sig3_6 sig3_7 sig3_8 sig3_9 sig3_10 sig3_11 sig3_12 ];

end

end

sig3_i = sig3(1,:) ; sig3_j = sig3(6,:); sig3_k = sig3(11,:) ; sig3_l = sig3(16,:);

sig_i = [sig1_i ; sig2_i; sig3_i];

sig_j = [sig1_j ; sig2_j; sig3_j];

sig_k = [sig1_k ; sig2_k; sig3_k];

sig_l = [sig1_l ; sig2_l; sig3_l];

for i = 1:((x-1)*(y-1))

displacements(:,i) = DELTA(reshape([3*KM(i,:)-2; 3*KM(i,:)-1; 3*KM(i,:)],[],1));

end

mxi = sig_i*displacements;

mxi = flip(reshape(mxi(1,:),(y-1),(x-1)));

Mxi = zeros(y,x);

Mxi(2:y,1:x-1) = mxi;

mxj = sig_j*displacements;

mxj = flip(reshape(mxj(1,:),(y-1),(x-1)));

Mxj = zeros(y,x);

Mxj(1:y-1,1:x-1) = mxj;

mxk = sig_k*displacements;

mxk = flip(reshape(mxk(1,:),(y-1),(x-1)));

Mxk = zeros(y,x);

Mxk(1:y-1,2:x) = mxk;

mxl = sig_l*displacements;

mxl = flip(reshape(mxl(1,:),(y-1),(x-1)));

Mxl = zeros(y,x);

Mxl(2:y,2:x) = mxl;

Mxx = (Mxi + Mxj + Mxk + Mxl)/-4;

Mxx(1,:) = (Mxj(1,:) + Mxk(1,:))/-2;

Mxx(end,:) = (Mxi(end,:) + Mxl(end,:))/-2;

Mxx(:,1) = (Mxi(:,1) + Mxj(:,1))/-2;

Mxx(:,end) = (Mxk(:,end) + Mxl(:,end))/-2;

Mxx(end,1) = -Mxi(end,1);

Mxx(1,1) = -Mxj(1,1);

Mxx(1,end) = -Mxk(1,end);

Mxx(end,end)= -Mxl(end,end);

myi = sig_i*displacements;

myi = flip(reshape(myi(2,:),(y-1),(x-1)));

Myi = zeros(y,x);

Myi(2:y,1:x-1) = myi;

myj = sig_j*displacements;

myj = flip(reshape(myj(2,:),(y-1),(x-1)));

Myj = zeros(y,x);

Myj(1:y-1,1:x-1) = myj;

myk = sig_k*displacements;

myk = flip(reshape(myk(2,:),(y-1),(x-1)));

Myk = zeros(y,x);

Myk(1:y-1,2:x) = myk;

myl = sig_l*displacements;

myl = flip(reshape(myl(2,:),(y-1),(x-1)));

Myl = zeros(y,x);

Myl(2:y,2:x) = myl;

Myy = (Myi + Myj + Myk + Myl)/-4;

Myy(1,:) = (Myj(1,:) + Myk(1,:))/-2;

Myy(end,:) = (Myi(end,:) + Myl(end,:))/-2;

Myy(:,1) = (Myi(:,1) + Myj(:,1))/-2;

Myy(:,end) = (Myk(:,end) + Myl(:,end))/-2;

Myy(end,1) = -Myi(end,1);

Myy(1,1) = -Myj(1,1);

Myy(1,end) = -Myk(1,end);

Myy(end,end)= -Myl(end,end);

f1= figure('NumberTitle', 'off', 'Name', 'Vertical Displacements (m)');

[XX,YY] = meshgrid(0:a:Lx, 0:b:Ly);

surf(XX,YY,Dz,'FaceColor', 'interp'); xlabel('x'); ylabel('y'); set(gca,'Ydir','reverse'); colorbar; colormap(flipud(jet)); view(2);

f2= figure('NumberTitle', 'off', 'Name', 'Mxx Moment Diagram (kNm/m)');

surf(XX,YY,Mxx,'FaceColor', 'interp'); xlabel('x'); ylabel('y'); set(gca,'Ydir','reverse'); colorbar; colormap(flipud(jet)); view(2);

f3= figure('NumberTitle', 'off', 'Name', 'Myy Moment Diagram (kNm/m)');

surf(XX,YY,Myy,'FaceColor', 'interp'); xlabel('x'); ylabel('y'); set(gca,'Ydir','reverse'); colorbar; colormap(flipud(jet)); view(2);

disp('Vertical Displacements - m')

disp(Dz)

disp('Mxx Moment Values - kNm/m')

disp(Mxx)

disp('Myy Moment Values - kNm/m')

disp(Myy);