Evaluation of the Stiffness Matrix and Stress Matrix by Gaussian Quadrature
30 views (last 30 days)
Show older comments
Hi guys,
Any help appreciated.
I am looking to build a code to get the stiffness matrix as in the example 10.4 of textbook by Logan. See the example attached.
I have made code to build one of the B matrices B(-0.5773, -0.5773) but how would I make a code to easily build the [k] matrix? I.e I would need to get the remaining B matrices for the various s and t, then build the [k] matrix as per the example.
Could anyone give me a suggestion?
This is my current code:
% Material properties
E = 30e6; v = 0.25;
D = (E/(1-v^2))*[1 v 0; v 1 0; 0 0 (1-v)/2];
% Coords
x1 = 1; x2 = 5; x3 = 5; x4 = 2;
X = [x1; x2; x3; x4];
y1 = 2; y2 = 2; y3 = 4; y4 = 4;
Y = [y1; y2; y3; y4];
% Jacobian
s = -0.5773; t = -0.5773;
J = 1/8*X.'*[0 1-t t-s s-1; t-1 0 s+1 -s-t; s-t -s-1 0 t+1; 1-s s+t -t-1 0]*Y;
% B Matrix
a = 1/4*(y1*(s-1) + y2*(-1-s) + y3*(1+s) + y4*(1-s));
b = 1/4*(y1*(t-1) + y2*(1-t) + y3*(1+t) + y4*(-1-t));
c = 1/4*(x1*(t-1) + x2*(1-t) + x3*(1+t) + x4*(-1-t));
d = 1/4*(x1*(s-1) + x2*(-1-s) + x3*(1+s) + x4*(1-s));
N1s = -(1-t)/4;
N2s = (1-t)/4 ;
N3s = (1+t)/4 ;
N4s = -(1+t)/4;
N1t = -(1-s)/4;
N2t = -(1+s)/4 ;
N3t = (1+s)/4 ;
N4t = (1-s)/4;
B1 = [a*N1s - b*N1t 0; 0 c*N1t - d*N1s; c*N1t - d*N1s a*N1s - b*N1t];
B2 = [a*N2s - b*N2t 0; 0 c*N2t - d*N2s; c*N2t - d*N2s a*N2s - b*N2t];
B3 = [a*N3s - b*N3t 0; 0 c*N3t - d*N3s; c*N3t - d*N3s a*N3s - b*N3t];
B4 = [a*N4s - b*N4t 0; 0 c*N4t - d*N4s; c*N4t - d*N4s a*N4s - b*N4t];
B = [B1 B2 B3 B4]
4 Comments
Reza Rezvani
on 9 Jun 2021
Hi Brian
Thank you for your answer.
I can help you, If you have a question about fracture, fatigue and crack propagation.
Best;
Reza
Reza Rezvani
on 12 Jun 2021
Brian
Hi
Do you know how can do it for Q8 and Q9?
I wanna to develope my code for Q8 and more nodes.
Answers (1)
Majeed
on 19 Apr 2022
Edited: Majeed
on 16 Sep 2022
Hey Brian,
I edited your code and it works well for me. Here it is and thank you.
% Material properties
E = 1; v = 0.3;
D = (E./(1-v.^2))*[1 v 0;v 1 0;0 0 ((1-v)./2)];
s1=-1/sqrt(3);
s2=1/sqrt(3);
t1=-1/sqrt(3);
t2=1/sqrt(3);
% Coords
x1 = 1; x2 = 5; x3 = 5; x4 = 2;
X = [x1; x2; x3; x4];
y1 = 2; y2 = 2; y3 = 4; y4 = 4;
Y = [y1; y2; y3; y4];
% Jacobian
syms s t;
J = 1/8*X.'*[0 1-t t-s s-1; t-1 0 s+1 -s-t; s-t -s-1 0 t+1; 1-s s+t -t-1 0]*Y;
% B Matrix
a = 1/4*(y1*(s-1) + y2*(-1-s) + y3*(1+s) + y4*(1-s));
b = 1/4*(y1*(t-1) + y2*(1-t) + y3*(1+t) + y4*(-1-t));
c = 1/4*(x1*(t-1) + x2*(1-t) + x3*(1+t) + x4*(-1-t));
d = 1/4*(x1*(s-1) + x2*(-1-s) + x3*(1+s) + x4*(1-s));
N1s = -(1-t)/4;
N2s = (1-t)/4 ;
N3s = (1+t)/4 ;
N4s = -(1+t)/4;
N1t = -(1-s)/4;
N2t = -(1+s)/4 ;
N3t = (1+s)/4 ;
N4t = (1-s)/4;
B1 = [a*N1s - b*N1t 0; 0 c*N1t - d*N1s; c*N1t - d*N1s a*N1s - b*N1t];
B2 = [a*N2s - b*N2t 0; 0 c*N2t - d*N2s; c*N2t - d*N2s a*N2s - b*N2t];
B3 = [a*N3s - b*N3t 0; 0 c*N3t - d*N3s; c*N3t - d*N3s a*N3s - b*N3t];
B4 = [a*N4s - b*N4t 0; 0 c*N4t - d*N4s; c*N4t - d*N4s a*N4s - b*N4t];
B = [B1 B2 B3 B4];
Bmat1= double(subs(B,[s,t],[s1,t1]));
Bmat2= double(subs(B,[s,t],[s1,t2]));
Bmat3=double(subs(B,[s,t],[s2,t1]));
Bmat4=double(subs(B,[s,t],[s2,t2]));
J1= double(subs(J,[s,t],[s1,t1]));
J2= double(subs(J,[s,t],[s1,t2]));
J3=double(subs(J,[s,t],[s2,t1]));
J4=double(subs(J,[s,t],[s2,t2]));
KE=Bmat1'*D*Bmat1*J1+Bmat2'*D*Bmat2*J2+Bmat3'*D*Bmat3*J3+Bmat4'*D*Bmat4*J4;
0 Comments
See Also
Categories
Find more on Regression in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!