# Matlab and finite element

1. Nov 26, 2012

1. The problem statement, all variables and given/known data
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[/COLOR];(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;

#### Attached Files:

• ###### s.PNG
File size:
17.8 KB
Views:
69
2. Nov 26, 2012