CEE470_HW3_Solution

pdf

School

University of Illinois, Urbana Champaign *

*We aren’t endorsed by this school

Course

470

Subject

Mechanical Engineering

Date

Apr 3, 2024

Type

pdf

Pages

Uploaded by CorporalSpiderMaster145

CEE 470 HW 3 clc; clear; Problem 1 syms c s E A L Consider the following problem: As we know, the elementwise stiffness matrix is shown below where is the angle between the positive x-axis and the element. K = E*A/L*[c*c s*c -c*c -s*c; s*c s*s -s*c -s*s; -c*c -s*c c*c s*c; -s*c -s*s s*c s*s] K = 1

Now, with this mind, let us define our material constants and elementwise values (length, angle, cos(angle), sin(angle)) tol = 1e-6; %Tolerance, used for making sure that values are approximately zero %Problem Description elem = 3; %Number of Elements ldof = 4; %number of local degrees of freedom gdof = 8; %number of global degrees of freedom fdof = [1,2]; %free degrees of freedom cdof = [3,4,5,6,7,8]; %constrained degrees of freedom %Material Constants A = 0.02; %m^2 E = 100e9; %Pa = kg/m/s^2 K = subs(K); %Element 1 Values L1 = 3; %Length of Element 1 [m] a1 = pi(); %Angle between positive x-axis and Element 1 [rad] c1 = cos(a1); %cosine of angle 1 s1 = sin(a1); %sine of angle 1 %Element 2 Values L2 = 1.5; %Length of Element 2 [m] a2 = pi()/2; %Angle between positive x-axis and Element 2 [rad] c2 = cos(a2); %cosine of angle 2 s2 = sin(a2); %sine of angle 2 %Element 3 Values L3 = 5; %Length of Element 3 [m] a3 = pi()+pi()/3; %Angle between positive x-axis and Element 3 [rad] c3 = cos(a3); %cosine of angle 3 s3 = sin(a3); %sine of angle 3 Assemble the Global Stiffness Matrix: To assemble the global stiffness matrix, we'll need to build each element's local stiffness matrix and then assemble. Let's first build each individual elementwise stiffness matrix: Element 1 Local Stiffness Matrix: K1 = subs(K,{L,c,s},[L1,c1,s1]) K1 = Element 2 Local Stiffness Matrix: 2

K2 = subs(K,{L,c,s},[L2,c2,s2]) K2 = Element 3 Local Stiffness Matrix: K3 = subs(K,{L,c,s},[L3,c3,s3]) K3 = And now, to assemble, we need to relate our local element nodes to the global degrees of freedom. We can do this using an Equivalent Degree of Freedom Table (EFT). EFT = zeros(elem,ldof); To be clear, let's build our EFT element by element. The row value will represent the element while the column value will represent the local degree of freedom. The input value will correspond to the global degree of freedom associated with the local degree of freedom. %Element 1 EFT(1,1) = 1; EFT(1,2) = 2; EFT(1,3) = 3; EFT(1,4) = 4; %Element 2 EFT(2,1) = 1; EFT(2,2) = 2; EFT(2,3) = 5; EFT(2,4) = 6; %Element 3 EFT(3,1) = 1; EFT(3,2) = 2; EFT(3,3) = 7; EFT(3,4) = 8; Now, we will use an algorithm to assemble our global stiffness matrix Kg = zeros(gdof); for i = 1:elem for j = 1:ldof for k = 1:ldof 3

Your preview ends here

Eager to read complete document? Join bartleby learn and gain access to the full version

Access to all documents
Unlimited textbook solutions
24/7 expert homework help

if i == 1 jj = EFT(i,j); kk = EFT(i,k); Kg(jj,kk) = Kg(jj,kk) + K1(j,k); %Pull from Element 1 elseif i == 2 jj = EFT(i,j); kk = EFT(i,k); Kg(jj,kk) = Kg(jj,kk) + K2(j,k); %Pull from Element 2 elseif i == 3 jj = EFT(i,j); kk = EFT(i,k); Kg(jj,kk) = Kg(jj,kk)+ K3(j,k); %Pull from Element 3 else disp( "WARNING: Element is outside of our Stiffness Matrix!" ) end end end end Kg Kg = 8×8 10 9 × 0.7667 0.1732 -0.6667 0.0000 -0.0000 -0.0000 -0.1000 -0.1732 0.1732 1.6333 0.0000 -0.0000 -0.0000 -1.3333 -0.1732 -0.3000 -0.6667 0.0000 0.6667 -0.0000 0 0 0 0 0.0000 -0.0000 -0.0000 0.0000 0 0 0 0 -0.0000 -0.0000 0 0 0.0000 0.0000 0 0 -0.0000 -1.3333 0 0 0.0000 1.3333 0 0 -0.1000 -0.1732 0 0 0 0 0.1000 0.1732 -0.1732 -0.3000 0 0 0 0 0.1732 0.3000 It is always a good check on our stiff matrix to make sure our stiffness matrix is symmetric. if (abs(det(Kg-Kg'))<tol) disp( "Stiffness Matrix is Symmetric!" ) else disp( "WARNING: Stiffness Matrix is NOT Symmetric" ) end Stiffness Matrix is Symmetric! Assemble the Global Load Vector We find this to be much easier than assembling the global stiffness matrix. Make sure that your units are consistent! Fg = zeros(gdof,1); Fg(2) = -50e3; %N Fg Fg = 8×1 0 -50000 0 0 0 0 0 4

0 That was easy! Solve the System Since we KNOW that our Stiffness Matrix is singular, we need to constrain our degree's of freedom. In this way, we will reduce our stiffness matrix and only our free dof's will be calculated. Kcc = Kg(cdof,cdof); Kcf = Kg(cdof,fdof); Kfc = Kg(fdof,cdof); Kff = Kg(fdof,fdof) %Free dof Stiffness Matrix Kff = 2×2 10 9 × 0.7667 0.1732 0.1732 1.6333 Fc = Fg(cdof); %Constrained dof Force Vector Ff = Fg(fdof) %Free dof Force Vector Ff = 2×1 0 -50000 u = zeros(gdof,1); uc = u(cdof); uf = Kff\(Ff-Kfc*uc); %Free dof's u(fdof) = uf u = 8×1 10 -4 × 0.0709 -0.3136 0 0 0 0 0 0 We can interpret this solution as node 1 moving 0.00709 mm to the right and 0.03136 mm down due to the 50 kN force in the negative y-direction on Node 1. Find the Element Forces Though it seems as if we are moving out of order, it is typically easier to find our element forces first, and then calculate the element stress. We will leverage the constrained dof's to find the resultant forces from the pins. Fc = Kcf*uf + Kcc*uc; %NOTE: Kcc*uc contributes nothing as uc = zeros Fg(cdof) = Fc Fg = 8×1 10 4 × 0 5

-5.0000 -0.4724 0.0000 0.0000 4.1818 0.4724 0.8182 Okay, these values seem to make sense but how can we check? Sum of Forces MUST be zero! %Sum of Forces in X-direction sumX = 0; for i=1:length(Fg)/2 sumX = sumX + Fg(2*i-1); end sumX; %Sum of Forces in Y-direction sumY = 0; for i=1:length(Fg)/2 sumY = sumY + Fg(2*i); end sumY; %Verify Equilibirum if (abs(sumX) < tol && abs(sumY) < tol) disp( "Structure is in Equilibirum!" ) else disp( "WARNING: Structure is NOT in Equilibrium" ) end Structure is in Equilibirum! Now that the forces in the pins have been found, we can determine the forces in each bar! This will be very easy as they can only carry axial forces! RECALL: Tension will be positive while compression will be negative. F_e1 = -Fg(3) %N F_e1 = 4.7240e+03 F_e2 = Fg(6) %N F_e2 = 4.1818e+04 F_e3 = -sqrt(Fg(7)*Fg(7)+Fg(8)*Fg(8)) %N F_e3 = -9.4478e+03 Find the Element Stresses In order to determine element stresses, we only need to divide our forces by the area. EASY! Sig_e1 = F_e1/A Sig_e1 = 2.3620e+05 Sig_e2 = F_e2/A 6

Your preview ends here

Eager to read complete document? Join bartleby learn and gain access to the full version

Access to all documents
Unlimited textbook solutions
24/7 expert homework help

Sig_e2 = 2.0909e+06 Sig_e3 = F_e3/A Sig_e3 = -4.7239e+05 Find Strain and Confirm Stresses As we all know, stress can also be found by using the global displacements to calculate the axial strain and then using the modulus of elasticity to calculate stress. T_1 = [-c1 -s1 0 0; 0 0 c1 s1]; u_1 = [u(1); u(2); u(3); u(4);]; eps_1 = T_1*u_1/L1 eps_1 = 2×1 10 -5 × 0.2362 0 Sig_1 = eps_1(1)*E %Matches! Sig_1 = 2.3620e+05 T_2 = [-c2 -s2 0 0; 0 0 c2 s2]; u_2 = [u(1); u(2); u(5); u(6);]; eps_2 = T_2*u_2/L2 eps_2 = 2×1 10 -4 × 0.2091 0 Sig_2 = eps_2(1)*E %Matches! Sig_2 = 2.0910e+06 T_3 = [-c3 -s3 0 0; 0 0 c3 s3]; u_3 = [u(1); u(2); u(7); u(8);]; eps_3 = T_3*u_3/L3 eps_3 = 2×1 10 -5 × -0.4724 0 Sig_3 = eps_3(1)*E %Matches! Sig_3 = -4.7239e+05 Check Overall Equilibrium of the Structure This has already been done after calculating the force vector solution, however, we will repeat it for clarity. 7

%Sum of Forces in X-direction sumX = 0; for i=1:length(Fg)/2 sumX = sumX + Fg(2*i-1); end sumX; %Sum of Forces in Y-direction sumY = 0; for i=1:length(Fg)/2 sumY = sumY + Fg(2*i); end sumY; %Verify Equilibirum if (abs(sumX) < tol && abs(sumY) < tol) disp( "Structure is in Equilibirum!" ) else disp( "WARNING: Structure is NOT in Equilibrium" ) end Structure is in Equilibirum! Problem 2 This problem is VERY similar to Problem 1 with the only difference being the y-settlement at Node 4. This will be accomodated in the force vector. Fg = zeros(gdof,1); Fg(2) = -50e3; Ff = Fg(fdof); Fc = Fg(cdof); u = zeros(gdof,1); u(8) = -0.4e-3; %m uc = u(cdof); Now that we've built the new displacement vector (and redeclared the force vector), let's solve the free displacements as before uf = Kff\(Ff-Kfc*uc) %Free dof's uf = 2×1 10 -4 × -0.6850 -0.9682 Look! our x-displacement is now negative at node 1, as we'd expect for this type of displacement! Now, determine the new set of forces based on the new displacements. Fc = Kcf*uf + Kcc*uc; %NOTE: Kcc*uc contributes nothing as uc = zeros Fg(cdof) = Fc Fg = 8×1 10 5 × 0 -0.5000 0.4566 8

-0.0000 0.0000 1.2909 -0.4566 -0.7909 Find the Element Stress in Element 3 To find the stress, let's first find the force, and then use the stress = F/A relationship to determine stress. F_e3 = sqrt(Fg(7)*Fg(7)+Fg(8)*Fg(8)) %N F_e3 = 9.1326e+04 Sig_e3 = F_e3/A Sig_e3 = 4.5663e+06 Find the Element Force in Element 2 F_e2 = Fg(6) F_e2 = 1.2909e+05 Find the Nodal Forces at Node 4 Node4_x = Fg(7) Node4_x = -4.5664e+04 Node4_y = Fg(8) Node4_y = -7.9090e+04 9

Your preview ends here

Eager to read complete document? Join bartleby learn and gain access to the full version