I encountered a problem in which:
I expect it to return -2*t*L^3/15.
The unfinished code is as follows, with the input of
syms L t positive real;shearcentre([0 0;0 L;L*sqrt(2) L/2],sym([3 1 3/4*t;1 2 t;2 3 3/4*t]))
on the command line.
Many thanks.
function [ec,final,qvec]=shearcentre(coor,barprop)% finding shear centre
arguments coor(:,2,1)sym % [x y]
barprop(:,3,1)sym % [n1 n2 t]
endclose all;[final,coor,barprop]=beamprop(coor,barprop); % barprop=[n1 n2 t length]
try double(coor); double(final); double(barprop);catch syms t L; assume([t L],{'real','positive'}); coor=simplify(coor); final=simplify(final); barprop=simplify(barprop);endif isAlways(final(4)~=final(3)) theta=atan(2*final(5)/(final(4)-final(3)))/2+pi;else theta=sym(pi);end coor=([cos(theta) sin(theta);-sin(theta) cos(theta)]*coor')';IX=final(3)*cos(theta)^2+final(4)*sin(theta)^2-2*final(5)*sin(theta)*cos(theta);IY=final(3)*sin(theta)^2+final(4)*cos(theta)^2+2*final(5)*sin(theta)*cos(theta);occur=sym(zeros(size(coor,1),1));for ct1=1:size(coor,1) for ct2=1:size(barprop,1) if barprop(ct2,1)==ct1 || barprop(ct2,2)==ct1 occur(ct1)=occur(ct1)+1; end endendif sum(isAlways(occur==0))>0 error('Error! Invalid inputs!');endqvec=sym(zeros(size(barprop,1),1));F=[qvec qvec qvec qvec];q=sym(zeros(size(coor,1),4));syms lv positive real;while size(find(occur==1,1),1)==1 start=find(occur==1,1); [mem,~]=find(barprop(:,1:2)==start); bar=barprop(mem,:); F(start(1),4)=-sign(sum([coor(bar(2),1)-coor(bar(1),1);coor(bar(2),2)-coor(bar(1),2)] ... .*[-coor(bar(1),2);coor(bar(1),1)])); if isAlways(start==bar(2))==1 bar(1:2)=fliplr(bar(1:2)); end % Vy
l1=coor(bar(1),2); l2=(coor(bar(1),2)+coor(bar(2),2))/2; lbar(lv)=((l2-l1)/bar(4))*lv+l1; A(lv)=bar(3)*lv; F(start,1)=int(q(start,1)+lbar*A,[0 bar(4)])/IX; q(bar(2),1)=q(start,1)+q(bar(2),1)+lbar(bar(4))*A(bar(4))/IX; qvec(start,1)=q(start,1)+lbar*A/IX; figure('NumberTitle','off'); fplot(qvec(start,1),[0 double(bar(4))]); title(sprintf('Vy, from node %g to node %g',start,bar(2))); grid minor; % Vx
l1=coor(bar(1),1); l2=(coor(bar(1),1)+coor(bar(2),1))/2; lbar(lv)=((l2-l1)/bar(4))*lv+l1; F(start,2)=int(q(start,3)+lbar*A,[0 bar(4)])/IY; q(bar(2),3)=q(start,3)+q(bar(2),3)+lbar(bar(4))*A(bar(4))/IY; qvec(start,2)=q(start,3)+lbar*A/IY; figure('NumberTitle','off'); fplot(qvec(start,2),[0 double(bar(4))]); title(sprintf('Vx, from node %g to node %g',start,bar(2))); grid minor; % Lever Arm (F(:,3))
m=(coor(bar(2),2)-coor(bar(1),2))/(coor(bar(2),1)-coor(bar(1),1)); if isAlways(m==0)==1 F(start,3)=abs(coor(bar(1),2)); elseif isAlways(isinf(m))==1 || isAlways(isinf(-m))==1 F(start,3)=abs(coor(bar(1),1)); else F(start,3)=abs((coor(bar(1),2)-m*coor(bar(1),2))/(m+1/m)*sqrt(1+m^-2)); end % finishing step
barprop=barprop([1:(mem-1) (mem+1):end],:); occur(start)=occur(start)-1; occur(bar(2))=occur(bar(2))-1;endif sum(occur)~=0 % box section
switch size(find(occur==3,1),1) case 0 syms qv positive real; qt=sym(zeros(size(barprop,1),5)); node=find(occur~=0); qt(:,1)=node'; start=min(node); qt(start,2)=-qv; qt(start,3)=-qv; [mem,~]=find(barprop(:,1:2)==start,1); barpropt=barprop; qti=sym([0;0]); while sum(occur)>0 % assume a cut in the min node
bar=barpropt(mem,:); if isAlways(start==bar(2))==1 bar(1:2)=fliplr(bar(1:2)); end % Vy l1=coor(bar(1),2); l2=(coor(bar(1),2)+coor(bar(2),2))/2; lbar(lv)=((l2-l1)/bar(4))*lv+l1; A(lv)=bar(3)*lv; qt(bar(2),2)=qt(start,2)+qt(bar(2),2)*(size(barpropt,1)>1)+lbar(bar(4))*A(bar(4)); qti(1)=qti(1)+int((qt(start,2)+lbar*A)/bar(3),lv,[0 bar(4)]); % Vx l1=coor(bar(1),1); l2=(coor(bar(1),1)+coor(bar(2),1))/2; lbar(lv)=((l2-l1)/bar(4))*lv+l1; qt(bar(2),3)=qt(start,3)+qt(bar(2),3)*(size(barpropt,1)>1)+lbar(bar(4))*A(bar(4)); qti(2)=qti(2)+int((qt(start,3)+lbar*A)/bar(3),lv,[0 bar(4)]); % finishing step occur(start)=occur(start)-1; occur(bar(2))=occur(bar(2))-1; barpropt(mem,:)=sym(zeros(1,4)); start=bar(2); if size(barpropt,1)>=1 [mem,~]=find(barpropt(:,1:2)==start,1); end end qti(1)=solve(qti(1)==0,qv); % problem encountered here.
qti(2)=solve(qti(2)==0,qv); case 2 syms q [2,1] positive real; % not finished
otherwise error('Error! Too complicated!'); endendex=sum(F(:,1).*F(:,3).*F(:,4))/IX;ey=sum(F(:,2).*F(:,3).*F(:,4))/IY;ec=vpa(([cos(theta) -sin(theta);sin(theta) cos(theta)]*[ex;ey]));endfunction [final,coor,barprop]=beamprop(coor,barprop)arguments coor(:,2,1)sym % [x y] barprop(:,3,1)sym % [n1 n2 t]endprop=sym(zeros(size(barprop,1),8));for ct=1:size(barprop,1) prop(ct,1)=sqrt((coor(barprop(ct,2),1)-coor(barprop(ct,1),1))^2+(coor(barprop(ct,2),2)-coor(barprop(ct,1),2))^2); % length
prop(ct,2)=prop(ct,1)*barprop(ct,3); % Area
prop(ct,3)=(coor(barprop(ct,1),1)+coor(barprop(ct,2),1))/2; % centre x-coor
prop(ct,4)=(coor(barprop(ct,1),2)+coor(barprop(ct,2),2))/2; % centre y-coor
if coor(barprop(ct,2),1)==coor(barprop(ct,1),1) prop(ct,5)=pi/2; % angle of bar
else prop(ct,5)=atan((coor(barprop(ct,2),2)-coor(barprop(ct,1),2))/(coor(barprop(ct,2),1)-coor(barprop(ct,1),1))); % angle of bar endendbarprop=[barprop prop(:,1)];final=sym(zeros(5,1));final(1)=sum(prop(:,2).*prop(:,3))/sum(prop(:,2));final(2)=sum(prop(:,2).*prop(:,4))/sum(prop(:,2));prop(:,3)=prop(:,3)-final(1);prop(:,4)=prop(:,4)-final(2);coor(:,1)=coor(:,1)-final(1);coor(:,2)=coor(:,2)-final(2);for ct=1:size(barprop,1) X=prop(ct,3)*cos(prop(ct,5))+prop(ct,4)*sin(prop(ct,5)); Y=-prop(ct,3)*sin(prop(ct,5))+prop(ct,4)*cos(prop(ct,5)); IX=prop(ct,1)*barprop(ct,3)*(barprop(ct,3)^2/12+Y^2); IY=prop(ct,1)*barprop(ct,3)*(prop(ct,1)^2/12+X^2); IXY=prop(ct,1)*barprop(ct,3)*X*Y; prop(ct,6)=IX*cos(-prop(ct,5))^2+IY*sin(-prop(ct,5))^2-2*IXY*sin(-prop(ct,5))*cos(-prop(ct,5)); % Ix
prop(ct,7)=IX*sin(-prop(ct,5))^2+IY*cos(-prop(ct,5))^2+2*IXY*sin(-prop(ct,5))*cos(-prop(ct,5)); % Iy
prop(ct,8)=(IX-IY)*sin(-prop(ct,5))*cos(-prop(ct,5))+IXY*(cos(-prop(ct,5))^2-sin(-prop(ct,5))^2); % Ixy
endfinal(3)=simplify(sum(prop(:,6)));final(4)=simplify(sum(prop(:,7)));final(5)=simplify(sum(prop(:,8)));% for ct=1:size(barprop,1)
% barprop(ct,6)=(coor(barprop(ct,1),1)+coor(barprop(ct,2),1))/2;
% barprop(ct,7)=(coor(barprop(ct,1),2)+coor(barprop(ct,2),2))/2;
% end
end
Best Answer