Evaluation of the Stiffness Matrix and Stress Matrix by Gaussian Quadrature

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]

[Reza Rezvani]

Hi

Also this is my question, Could you please help me if you found the solution

Email: [e-mail]

Best;

[Brian Robinson]

Hi Reza,

This is the best I could do. You will need to modify the coordinates and properties to suit.

% Material properties
E = 50e9; v = 0.2;
% Use plane stress
D = (E/(1-v^2))*[1 v 0; v 1 0; 0 0 (1-v)/2]; 

% Coords
x1 = 1; x2 = 5; x3 = 5; x4 = 1;
X = [x1; x2; x3; x4];
y1 = 2; y2 = 2; y3 = 5; y4 = 4; 
Y = [y1; y2; y3; y4];

% Combinations of -0.5773 and 0.5773
Gauss_points = [-0.5773,-0.5773; -0.5773, 0.5773; 0.5773, -0.5773; 0.5773 0.5773];

k = zeros(8);
for i = 1:4
    s = Gauss_points(i,1);
    t = Gauss_points(i,2);
    %disp(s);
    %disp(s);
    % Jacobian
    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;

    %disp(J);
    % 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));

    %Coeffs = [a b c d];
    %disp(Coeffs);
    
    N1s = (t-1)/4;
    N2s = (1-t)/4 ;
    N3s = (1+t)/4 ;
    N4s = -(1+t)/4;
    N1t = (s-1)/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];
    ki = B.'*D*B*J;
    k = k + ki;
end
%MaxDiag = max(diag(k));
%disp(MaxDiag)
%disp(k);

[Reza Rezvani]

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]

Brian

Hi

Do you know how can do it for Q8 and Q9?

I wanna to develope my code for Q8 and more nodes.

Answer 1

[Majeed]

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;