## 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?

### Shubham Gorde (view profile)

on 8 May 2018
Latest activity Commented on by Shubham Gorde

on 17 Jun 2018

### Walter Roberson (view profile)

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)

### Walter Roberson (view profile)

on 8 May 2018

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)

Shubham Gorde

### Shubham Gorde (view profile)

on 16 Jun 2018
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?
Walter Roberson

### Walter Roberson (view profile)

on 17 Jun 2018
I think I am a bit too tired to look at this today, sorry.
Shubham Gorde

### Shubham Gorde (view profile)

on 17 Jun 2018
Will you able to check this today?