-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathStiff_HCT_Sub_T3_SP.m
58 lines (46 loc) · 5.07 KB
/
Stiff_HCT_Sub_T3_SP.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
function [ ke ] = Stiff_HCT_Sub_T3_SP( edof,H,x,y,gauss,we,ke )
L1=gauss(1); L2=gauss(2); L3=gauss(3);
x1=x(1); x2=x(2); x3=x(3); y1=y(1); y2=y(2); y3=y(3);
x0=(x1+x2+x3)/3; y0=(y1+y2+y3)/3;
% A32=2*area_Sub_triangle_3
% o: The centroid of triangle element
% o-node (3)
% o
% / \
% L2=0 / \ L1=0
% / Sub-T_1 \
% 1-node (1) o-------------o 3-node (2)
% L3=0
% x0=(x1+x2+x3)/3; y0=(y1+y2+y3)/3;
b20= y2-y0; b01=y0-y1; b12=y1-y2;
c20= x2-x0; c01=x0-x1; c12=x1-x2;
l20=sqrt(b20^2+c20^2); l01=sqrt(b01^2+c01^2); l12=sqrt(b12^2+c12^2);
% mu20=(l12^2-l01^2)/l20^2; mu01=(l20^2-l12^2)/l01^2;
mu12=(l01^2-l20^2)/l12^2;
A32=b20*(-c01)-(-c20)*b01; % A22=2*area_Sub_triangle_3
A34=2*A32;
Q=[ b20 b01 b12 ;...
-c20 -c01 -c12];
% Qui uoc , w_x=-dw/dy, w_y=dw/dx
D2L11 =[ 6*L1 - (6*L3*(b01*c20 - b20*c01))/(b01*c12 - b12*c01) + (6*L2*(b12*c20 - b20*c12))/(b01*c12 - b12*c01), (2*A32*L3*b01)/(b01*c12 - b12*c01) - (2*A32*L2*b12)/(b01*c12 - b12*c01), (2*A32*L2*c12)/(b01*c12 - b12*c01) - (2*A32*L3*c01)/(b01*c12 - b12*c01), 0, 0, 0, 0, 0, 0];
D2L22 =[ 0, 0, 0, 6*L2 + (6*L1*(b01*c12 - b12*c01))/(b12*c20 - b20*c12) - (6*L3*(b01*c20 - b20*c01))/(b12*c20 - b20*c12), (2*A32*L1*b12)/(b12*c20 - b20*c12) - (2*A32*L3*b20)/(b12*c20 - b20*c12), (2*A32*L3*c20)/(b12*c20 - b20*c12) - (2*A32*L1*c12)/(b12*c20 - b20*c12), 0, 0, 0];
D2L33 =[ 0, 0, 0, 0, 0, 0, 6*L3 - (6*L1*(b01*c12 - b12*c01))/(b01*c20 - b20*c01) - (6*L2*(b12*c20 - b20*c12))/(b01*c20 - b20*c01), (2*A32*L1*b01)/(b01*c20 - b20*c01) - (2*A32*L2*b20)/(b01*c20 - b20*c01), (2*A32*L2*c20)/(b01*c20 - b20*c01) - (2*A32*L1*c01)/(b01*c20 - b20*c01)];
D2L12 =[ (3*L3*(b01*c12 - b12*c01 + 2*b01*c20 - 2*b20*c01 + 3*b12*c20 - 3*b20*c12 - b01*c12*mu12 + b12*c01*mu12 - b12*c20*mu12 + b20*c12*mu12))/(2*(b01*c12 - b12*c01)) + (6*L1*(b12*c20 - b20*c12))/(b01*c12 - b12*c01), (A34*L3*c12)/l12^2 - (A32*L3*(2*b01 + 3*b12 - b12*mu12))/(2*(b01*c12 - b12*c01)) - (2*A32*L1*b12)/(b01*c12 - b12*c01), (A32*L3*(2*c01 + 3*c12 - c12*mu12))/(2*(b01*c12 - b12*c01)) + (A34*L3*b12)/l12^2 + (2*A32*L1*c12)/(b01*c12 - b12*c01), (3*L3*(3*b01*c12 - 3*b12*c01 + 2*b01*c20 - 2*b20*c01 + b12*c20 - b20*c12 + b01*c12*mu12 - b12*c01*mu12 + b12*c20*mu12 - b20*c12*mu12))/(2*(b12*c20 - b20*c12)) + (6*L2*(b01*c12 - b12*c01))/(b12*c20 - b20*c12), (A32*L3*(3*b12 + 2*b20 + b12*mu12))/(2*(b12*c20 - b20*c12)) + (A34*L3*c12)/l12^2 + (2*A32*L2*b12)/(b12*c20 - b20*c12), (A34*L3*b12)/l12^2 - (A32*L3*(3*c12 + 2*c20 + c12*mu12))/(2*(b12*c20 - b20*c12)) - (2*A32*L2*c12)/(b12*c20 - b20*c12), 0, 0, 0];
D2L13 =[ (3*L2*(b01*c12 - b12*c01 + 2*b01*c20 - 2*b20*c01 + 3*b12*c20 - 3*b20*c12 - b01*c12*mu12 + b12*c01*mu12 - b12*c20*mu12 + b20*c12*mu12))/(2*(b01*c12 - b12*c01)) - (6*L1*(b01*c20 - b20*c01))/(b01*c12 - b12*c01), (A34*L2*c12)/l12^2 - (A32*L2*(2*b01 + 3*b12 - b12*mu12))/(2*(b01*c12 - b12*c01)) + (2*A32*L1*b01)/(b01*c12 - b12*c01), (A32*L2*(2*c01 + 3*c12 - c12*mu12))/(2*(b01*c12 - b12*c01)) + (A34*L2*b12)/l12^2 - (2*A32*L1*c01)/(b01*c12 - b12*c01), (3*L2*(3*b01*c12 - 3*b12*c01 + 2*b01*c20 - 2*b20*c01 + b12*c20 - b20*c12 + b01*c12*mu12 - b12*c01*mu12 + b12*c20*mu12 - b20*c12*mu12))/(2*(b12*c20 - b20*c12)), (A32*L2*(3*b12 + 2*b20 + b12*mu12))/(2*(b12*c20 - b20*c12)) + (A34*L2*c12)/l12^2, (A34*L2*b12)/l12^2 - (A32*L2*(3*c12 + 2*c20 + c12*mu12))/(2*(b12*c20 - b20*c12)), -(6*L3*(b01*c12 - b12*c01))/(b01*c20 - b20*c01), (2*A32*L3*b01)/(b01*c20 - b20*c01), -(2*A32*L3*c01)/(b01*c20 - b20*c01)];
D2L23 =[ (3*L1*(b01*c12 - b12*c01 + 2*b01*c20 - 2*b20*c01 + 3*b12*c20 - 3*b20*c12 - b01*c12*mu12 + b12*c01*mu12 - b12*c20*mu12 + b20*c12*mu12))/(2*(b01*c12 - b12*c01)), (A34*L1*c12)/l12^2 - (A32*L1*(2*b01 + 3*b12 - b12*mu12))/(2*(b01*c12 - b12*c01)), (A32*L1*(2*c01 + 3*c12 - c12*mu12))/(2*(b01*c12 - b12*c01)) + (A34*L1*b12)/l12^2, (3*L1*(3*b01*c12 - 3*b12*c01 + 2*b01*c20 - 2*b20*c01 + b12*c20 - b20*c12 + b01*c12*mu12 - b12*c01*mu12 + b12*c20*mu12 - b20*c12*mu12))/(2*(b12*c20 - b20*c12)) - (6*L2*(b01*c20 - b20*c01))/(b12*c20 - b20*c12), (A32*L1*(3*b12 + 2*b20 + b12*mu12))/(2*(b12*c20 - b20*c12)) + (A34*L1*c12)/l12^2 - (2*A32*L2*b20)/(b12*c20 - b20*c12), (A34*L1*b12)/l12^2 - (A32*L1*(3*c12 + 2*c20 + c12*mu12))/(2*(b12*c20 - b20*c12)) + (2*A32*L2*c20)/(b12*c20 - b20*c12), -(6*L3*(b12*c20 - b20*c12))/(b01*c20 - b20*c01), -(2*A32*L3*b20)/(b01*c20 - b20*c01), (2*A32*L3*c20)/(b01*c20 - b20*c01)];
% Extract from formulation (4.58), FEM - Zienkiewicz vol.2
for i=1:edof
DNL=zeros(3,3);
DNL=[D2L11(i) D2L12(i) D2L13(i);...
D2L12(i) D2L22(i) D2L23(i);...
D2L13(i) D2L23(i) D2L33(i)];
% Second derivatives of shape function coresponding to x,y in (4.58), pp.132
DNxy=Q*DNL*Q';
B(1,i)=DNxy(1,1);
B(2,i)=DNxy(2,2);
B(3,i)=2*DNxy(1,2);
end
B=(1/A32^2)*B;
% Formula (27.32) in IFEM,chapter 27.
ke=ke+0.5*A32*we*B'*H*B;
end