1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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:
      72
  2. jcsd
  3. Nov 26, 2012 #2
    so i want to change this from point load to distributed load
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Matlab and finite element
  1. Finite Element (Replies: 3)

  2. Finite element meshing (Replies: 0)

Loading...