- #1
Sadeq
- 103
- 0
Homework Statement
The problem picture is attached(file 1),its a beam subjected to horizonatal ditributed load
2. Relevant examples
the MATLAB solution for rectangular shape with vertical load on the upper right corner is like follow, i try to modify it according to the new picture, but i struggle in the load
clear all
close all
clc
%% Constants
L = 2;
D = 0.4;
t = 0.005;
E = 200e9;
nu = 0.3;
P = 10000;
n_x = 2;
n_y = 2;
C = E/(1-nu^2)*[1 nu,0;nu,1,0;0,0,(1-nu)/2];
K_G=zeros(2*(n_x+1)*(n_y+1),2*(n_x+1)*(n_y+1)); % The global stiffness matrix
%% Gauss-Legendre sampling points and weights
alpha = [1,1];
r = [-0.57735026919,0.57735026919];
s = [-0.57735026919,0.57735026919];
%% Calculate the global stiffness matrix
for i=1:n_x
for j=1:n_y
K_local = zeros(8,8);
%Note that the nested loop here is essentially cycling through the
%summation required for the numerical integration
for p=1:2
for q=1:2
x=[i*L/n_x,(i-1)*L/n_x,(i-1)*L/n_x,i*L/n_x]; %Global x-coordinate
y=[j*D/n_y,j*D/n_y,(j-1)*D/n_y,(j-1)*D/n_y]; %Global y-coordinate
dxdr=1/4*(1+s(1,q))*x(1,1)-1/4*(1+s(1,q))*x(1,2)-1/4*(1- s(1,q))*x(1,3)+1/4*(1-s(1,q))*x(1,4);
dxds=1/4*(1+r(1,p))*x(1,1)+1/4*(1-r(1,p))*x(1,2)-1/4*(1-r(1,p))*x(1,3)-1/4*(1+r(1,p))*x(1,4);
dydr=1/4*(1+s(1,q))*y(1,1)-1/4*(1+s(1,q))*y(1,2)-1/4*(1-s(1,q))*y(1,3)+1/4*(1-s(1,q))*y(1,4);
dyds=1/4*(1+r(1,p))*y(1,1)+1/4*(1-r(1,p))*y(1,2)-1/4*(1-r(1,p))*y(1,3)-1/4*(1+r(1,p))*y(1,4);
J=[dxdr,dydr;dxds,dyds];
invJ=inv(J);
detJ=det(J);
I1=.25*[1+s(1,q),0,-1-s(1,q),0,-1+s(1,q),0,1-s(1,q),0;1+r(1,p),0,1-r(1,p),0,-1+r(1,p),0,-1-r(1,p),0];
I2=.25*[0,1+s(1,q),0,-1-s(1,q),0,-1+s(1,q),0,1-s(1,q);0,1+r(1,p),0,1-r(1,p),0,-1+r(1,p),0,-1-r(1,p)];
B=[1,0;0,0;0,1]*invJ*I1+[0,0;0,1;1,0]*invJ*I2;
%Note here that we're adding to the previous value calculated for K_local (Gauss-Legendre scheme)
K_local=K_local+alpha(1,p)*alpha(1,q)*B.'*C*B*t*detJ;
end
end
nb=1+2*(i-1)+2*(j-1)*(n_x+1); %The bottom left global DoF number for element i,j
nt=1+2*(i-1)+2*j*(n_x+1); %The top left global DoF number for element i,j
%Assign the elements of the local stiffness matrix to the relevant location in the global stiffness matrix
K_G([[nb:1:nb+3],[nt:1:nt+3]],[[nb:1:nb+3],[nt:1:nt+3]])=K_G([[nb:1:nb+3],[nt:1:nt+3]],[[nb:1:nb+3],[nt:1:nt+3]])+K_local([5,6,7,8,3,4,1,2],[5,6,7,8,3,4,1,2]);
end
end
%% Construct the force vector
F=zeros(2*(n_x+1)*(n_y+1),1); %Construct the global force vector (zeros only)
F(2*(n_x+1)*(n_y+1),1)=-P;(how can i put the distributed load here) %Add the applied load to the global force vector
dof=[1:2*(n_x+1)*(n_y+1)]; %Construct a vector of DoF numbers
%Note
%C = setdiff(A,B) for vectors A and B, returns the values in A that are not in B with no repetitions.
dof=setdiff(dof,[1:2*(n_x+1):2*(n_x+1)*(n_y+1)]); %Extract/elimiate the horizontal DoF at the support
dof=setdiff(dof,[2:2*(n_x+1):2*(n_x+1)*(n_y+1)]); %Eliminate the vertical DoF at the support, leaving only the 'active' DoF
F2=F(dof,1); %Determine the reduced force matrix
K_G2=K_G(dof,dof); %Determine the reduced stiffness matrix
d2=inv(K_G2)*F2; %Calculate the vector of displacements
d=zeros(2*(n_x+1)*(n_y+1),1); %Produce a vector of zeros corresponding to the number of DoF
d(dof,1)=d2; %Add the calculated displacements to the previous vector
uy1=d;
