Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Homework Help: Matlab and finite element

  1. Nov 26, 2012 #1
    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
      s.PNG
      File size:
      17.8 KB
      Views:
      81
  2. jcsd
  3. Nov 26, 2012 #2
    so i want to change this from point load to distributed load
     
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook