MATLAB Answers

0

for the force vector part ri and stiffness matrix ki,error is getting displayed while plotting the plots for displacement graph analytically and by approximation methods.What seems to be the probelm,if found ,please explain?

Asked by Shubham Gorde on 8 May 2018
Latest activity Commented on by Shubham Gorde on 17 Jun 2018
syms pi p q E A l x
n=3;
k=(n:n);
r=(1:n);
a=(1:n);
N=sym(1:n);
pi=eval(pi);
p=q*x/E*A*l;
for i=1:n  
      syms x
      N(i)=(x.^2-2*x).^i;
      %N(i)=sin((2*i-1)/2*pi*x);
      %N(i)=piecewise(x>=(i-1)/n & x<1/n+(i-1)/n,n*(x-(i-1)/n),x>=i/n & x<=(i+1)/n & x<1,-n*x+i+1,0);
      Y1=(N(i));
      Y2=diff(N(i),x);
      Y3=diff(N(i),x,2);
      x=linspace(0,1);
      plot(x,subs(Y1))
      hold on
      plot(x,subs(Y2))
      hold on
      plot(x,subs(Y3))
      hold off
      legend('N','N''','N''''');
      pause(1)
end
for i=1:n
      r(i)=-int(N(i)*p,x=linspace(0,1))*l ;   %Galerkin 
      %r(i)=-int(diff(N(i),x,2)*p,x=linspace(0,1))/l ;     %Least squares
      for j=1:n
          %k(i,j)=1/l*int(N(i)*diff(N(i),x,2),x=linspace(0,1));      %Galerkin
           %k(i,j)=1/l^3*int(diff (N(j),x,2)*diff(N(i),x,2),x=linspace(0,1));     %Least Squares
           k(i,j)=-1/l*int(diff(N(i),x)*diff(N(j),x),x=linspace(0,1));     %Displacement formulation
      end
end
display(k)
display(r)
a=linsolve(k,r);
U=q*l^2*(x/2-x^3/6)/E*A;
u=o;
for i=1:n
      u=u+a(i)*N(i);
end
u=simplify(u);
plot((U/(q*l^2/E*A)),(u/(q*l^2/E*A),x=linspace(0,1))
X1=0:0.9999;
Z3=(U/(q*l^2/E/A)-u/(q*l^2/E/A));
plot(X1,Z3)

  0 Comments

Sign in to comment.

1 Answer

Answer by Walter Roberson
on 8 May 2018
 Accepted Answer

This is as close as you can get:

syms pi p q E A l x
n=3;
k = sym(n:n);
r = sym(1:n);
a=(1:n);
N=sym(1:n);
pi=eval(pi);
p=q*x/E*A*l;
for i=1:n  
        syms x
        N(i)=(x.^2-2*x).^i;
        %N(i)=sin((2*i-1)/2*pi*x);
        %N(i)=piecewise(x>=(i-1)/n & x<1/n+(i-1)/n,n*(x-(i-1)/n),x>=i/n & x<=(i+1)/n & x<1,-n*x+i+1,0);
        Y1=(N(i));
        Y2=diff(N(i),x);
        Y3=diff(N(i),x,2);
        x=linspace(0,1);
        plot(x,subs(Y1))
        hold on
        plot(x,subs(Y2))
        hold on
        plot(x,subs(Y3))
        hold off
        legend('N','N''','N''''');
        pause(1)
end
syms x
for i=1:n
        r(i)=-int(N(i)*p,x,0, 1)*l ;   %Galerkin 
        %r(i)=-int(diff(N(i),x,2)*p,x=linspace(0,1))/l ;     %Least squares
        for j=1:n
            %k(i,j)=1/l*int(N(i)*diff(N(i),x,2),x=linspace(0,1));      %Galerkin
             %k(i,j)=1/l^3*int(diff (N(j),x,2)*diff(N(i),x,2),x=linspace(0,1));     %Least Squares
             k(i,j)=-1/l*int(diff(N(i),x)*diff(N(j),x),x, 0, 1);     %Displacement formulation
        end
end
display(k)
display(r)
a=linsolve(k,r.');
U=q*l^2*(x/2-x^3/6)/E*A;
u = sym(0);
for i=1:n
        u=u+a(i)*N(i);
end
u=simplify(u);
X = linspace(0,1);
expr1 = subs((U/(q*l^2/E*A)), x, X);
expr2 = subs((u/(q*l^2/E*A)), x, X);
plot(expr1, expr2)
X1=0:0.9999;
Z3 = subs((U/(q*l^2/E/A)-u/(q*l^2/E/A)), x, X1);
plot(X1,Z3)

The two plot() at the end will fail. The first one will fail because expr2 involves the unresolved symbolic variable l (lower-case L). The second one will fail because Z3 involves the unresolved symbolic variables A and l (lower-case L)

  44 Comments

Actually I am solving plane stress triangular element problem.I solved this in this manner but want to standardized it for more number of elements.My code is this:

syms x1 x2 x3 x4 y1 y2 y3 y4 x4 z1 z2 z3 z4
A =[z1 x1 y1; z2 x2 y2 ;z3 x3 y3];
E=10000; %Modulus of elasticity
h=1; %thickness
V=0.1; %poissons ratio
%for Element 1
Area1 = det(A);
X11=0; Y11=0; Z11=1; %these are the coordinates for element 1
X21=100; Y21=0; Z21=1;
X31=100; Y31=10; Z31=1;
this_Area1 = 0.5*subs(Area1, [x1, x2, x3, y1, y2, y3, z1, z2, z3],
 [X11, X21, X31, Y11, Y21,Y31, Z11, Z21, Z31]);
b11= X21*Y31-X31*Y21;
b21= X31*Y11-X11*Y31;
b31= X11*Y21-X21*Y11;
c11= Y21-Y31;
c21= Y31-Y11;
c31= Y11-Y21;
d11= X31-X21;
d21= X11-X31;
d31= X21-X11;
%For Element 2
A =[z1 z3 z4; x1 x3 x4 ; y1 y3 y4];
Area2 = det(A);
X12=0; Y12=0; Z12=1;     % these are the coordinates for element 2
X32=100; Y32=10; Z32=1;
X42=0; Y42=10; Z42=1;
this_Area2 =0.5* subs(Area2, [x1, x3, x4, y1, y3, y4, z1, z3, z4], 
[X12, X32, X42, Y12, Y32,Y42, Z12, Z32, Z42]);
b12= X32*Y42-X42*Y32;      % b12 means b1 value for element 2,likewise all b,c d values for element 2
b32= X42*Y12-X12*Y42;
b42= X12*Y32-X32*Y12;
c12= Y32-Y42;
c32= Y42-Y12;
c42= Y12-Y32;
d12= X42-X32;
d32= X12-X42;
d42= X32-X12;
%for element 1 and 2
      k221=E*h/4/this_Area1/(1-V.^2)*([c21*c21+d21*d21*(1-V)/2  c21*d21*V+c21*d21*(1-V)/2;
                              c21*d21*V+c21*d21*(1-V)/2  d21*d21+c21*c21*(1-V)/2]); % matrix for k22 matrix for element 1
      k231=E*h/4/this_Area1/(1-V.^2)*[c21*c31+d21*d31*(1-V)/2  c21*d31*V+c31*d21*(1-V)/2;
                              c31*d21*V+c21*d31*(1-V)/2  d21*d31+c21*c31*(1-V)/2]; % matrix for k23 matrix for element 1
    k321=E*h/4/this_Area1/(1-V.^2)*[c31*c21+d31*d21*(1-V)/2  c31*d21*V+c21*d31*(1-V)/2;
                              c21*d31*V+c31*d21*(1-V)/2  d31*d21+c31*c21*(1-V)/2]; % matrix for k32 matrix for element 1
   k331=E*h/4/this_Area1/(1-V.^2)*[c31*c31+d31*d31*(1-V)/2  c31*d31*V+c31*d31*(1-V)/2;
                              c31*d31*V+c31*d31*(1-V)/2  d31*d31+c31*c31*(1-V)/2]; % matrix for k33 matrix for element 1
k332=E*h/4/this_Area2/(1-V.^2)*[c32*c32+d32*d32*(1-V)/2  c32*d32*V+c32*d32*(1-V)/2;
                            c32*d32*V+c32*d32*(1-V)/2  d32*d32+c32*c32*(1-V)/2]; % matrix for k33 matrix for element 2
KG=[k221 k231; k321 k331+k332];
R=[100; 0 ;100; 0];
a=(KG)\R

I am getting the same answer as pdf but b,c,d,stiffness matrix are of manual input and this is limited to this particular problem. so wanted to make standard procedure for many problems.

The solution above u have given its not giving b1,b2,b3,c1,c2,c3,d1,d2,d3 values?

Sign in to comment.