تحلیل خرپا دوبعدی به روش اجزاء محدود در متلب
در این بخش به تحلیل خرپای دوبعدی شکل 1 با استفاده از یک برنامه کامپیوتری که در محیط نرم افزار MATLAB نوشته و اجرا میشود، پرداخته خواهد شد. مقصود از نوشتن برنامه کامپیوتری، به کارگیری آن برای محاسبه مقادیر گره ای جابجایی و نیروهای تکیه گاهی سپس محاسبه نیرو، تنش و کرنش در تمامی اعضای خرپا دوبعدی میباشد. در این برنامه بخش های مختلفی وجود دارد. این بخشها شامل بدنه اصلی برنامه، فایل های ورودی، توابع و خروجی ها می باشند که در ادامه به آنها پرداخته خواهد شد. لازم به ذکر است پروژه متلب ذکر شده در این محتوا توسط تیم حامی پیپر انجام شده است.
1-3- بدنه اصلی برنامه:
بدنه اصلی این برنامه، «Mainbody» نام دارد و تابع فرمان برنامه میباشد. وظیفه مدیریت کلیه وضایف برنامه، برعهده این قسمت از برنامه می باشد. کلیه توابع موجود در این برنامه توسط این بخش فراخوانی میشوند و ورودیها و خروجیهای این توابع در این قسمت ثبت و فراخوانی میشوند. در این برنامه موقعیت گرهها، جدول اطلاعات المانی و خواص مکانیکی اعضای خرپا و شرایط مرزی مسئله به عنوان ورودی به صورت فایل .txt وجود دارند و از توابع مختلفی برای محاسبه طول المانها، محاسبه زاویه المانها، محاسبه ماتریس سفتی کلی، اعمال شرایط مرزی و محاسبه نیروهای تکیه گاهی، نیروهای المانی، تنش و کرنش المانها استفاده شده است. روند کلی برنامه در روندنمای شکل 9 قابل مشاهده است.
2-3- فایلهای ورودی:
از 4 فایل به عنوان ورودی به فرمت.txt استفاده شده است. این فایلها به عنوان ورودی مورد استفاده قرار میگیرند و دادههای مورد نیاز اولیه نظیر شماره گرهها، شماره المانها، مشخصات المانها و شرایط مرزی مسئله برای خرپای شکل 1 از این فایلها استخراج میشوند.
1-2-3- فایل Nodes Positions.txt: در این فایل شماره گرهها و موقعیت آنها در سه ستون (شکل10) آورده شده است. ستون اول شماره گره، ستون دوم موقعیت در راستای محور X و ستون سوم موقعیت در راستای محور Y می باشند. دادههای این فایل در متغیر Nodes در بدنه اصلی برنامه ذخیره میشوند.
2-2-3- فایل Element Property And connectivity Data.txt: در این فایل شماره المانها به همراه شماره گرههایی که المان بین آنها قرار دارد در کنار مساحت سطح مقطع و مدول یانگ اعضای خرپا آورده شده است (شکل 11). ستون اول شماره المان، ستون دوم و سوم شماره گرههای المان، ستون چهارم مدول یانگ و ستون پنجم مساحت مقطع المان میباشد. دادههای این فایل در متغیر Elem در بدنه اصلی برنامه ذخیره میشوند.
3-2-3- فایل BC_Force.txt: در این فایل شرایط مرزی از نوع نیرو مشخص میشود و شامل سه ستون است (شکل 12). ستون اول مربوط به شماره گره، ستون دوم مقدار نیرو در راستای X و ستون سوم مقدار نیرو در راستای Y میباشد. دادههای این فایل در متغیر BCF در بدنه اصلی برنامه ذخیره میشوند.
شکل 9 روند کلی برنامهMainbody
شکل 10 فایل ورودی Nodes Positions.txt
شکل 11 فایل ورودی Element Property And connectivity Data.txt
شکل 12 فایل ورودی BC_Force.txt
3-2-4- فایل BC_Displacement.txt: در این فایل شرایط مرزی از نوع جابجایی مشخص میشود و شامل سه ستون است (شکل 13). ستون اول مربوط به شماره گره، ستون دوم مقدار جابجایی در راستای X و ستون سوم مقدار جابجایی در راستای Y میباشد. برای گرههایی که به آنها نیرو وارد میشود عدد 99999 به عنوان یک نماد در ستونهای دوم و سوم قرار داده شده است. دادههای این فایل در متغیر BCD در بدنه برنامه ذخیره میشوند.
شکل 13 فایل ورودی BC_Displacement.txt
3-3- توابع استفاده شده در برنامه:
در این برنامه برای محاسبه پارامترهای مختلف از برخی توابع استفاده شده است. این کار هم باعث درک راحتتر روند برنامه و هم زیباتر شدن متن برنامه میشود.
1-3-3- تابع محاسبه زاویه المانها:
در این برنامه با استفاده از تابع Alfa، زاویه المانها محاسبه میشود. ورودیهای این تابع اطلاعات ورودی گرهها (Nodes)، تعداد المانها (e) و اطلاعات ورودی المانها (Elem) میباشند. متن این تابع در ادامه مشاهده میشود. خروجی این تابع در متغیر ElemAlfa ذخیره میشود.
شاید این مطالب نیز برای شما جذاب باشد، پیشنهاد میکنیم به این صفحات نیز سر بزنید:
functionElemAlfa= Alfa(e, Elem, Nodes)
ElemAlfa=zeros(e,1);%The angle of elements
for k=1:e
i= Elem(k,2); %node i of element e
j= Elem(k,3); %node j of element e
ElemAlfa(k,1)=atand((Nodes(j,3)-Nodes(i,3))/(Nodes(j,2)-Nodes(i,2)));
% arctan((yj- yi)/(xj- xi))
end
2-3-3- تابع محاسبه طول المانها:
در این برنامه با استفاده از تابع ElemLength، طول المانها محاسبه میشود. ورودیهای این تابع تعداد المانها (e)، اطلاعات ورودی گرهها (Nodes) و اطلاعات ورودی المانها (Elem) میباشند. متن این تابع در ادامه مشاهده میشود. خروجی این تابع در متغیر L ذخیره می شود.
function L=ElemLength(e,Nodes,Elem)
L=zeros(e,1);
for m=1:e
i= Elem(m,2); %node i of element e
j= Elem(m,3); %node j of element e
L(m,1)=sqrt((Nodes(j,3)-Nodes(i,3))^2+(Nodes(j,2)-Nodes(i,2))^2);
end
3-3-3- تابع استخراج ماتریس سفتی کل:
در این برنامه با استفاده از تابع ElemStiffAssamblage، ابتدا از رابطه 8 ماتریس سختی المانی در مختصات محلی (KLoc) استخراج میشود. سپس با استفاده از رابطه 18 ماتریس سختی المانی در مختصات کلی(KElem) محاسبه و در نهایت در ماتریس سختی کل (KTotal) مونتاژ میشود. این فرآیند برای تمام المانها تکرار شده و ماتریس سختی کلی تشکیل میشود. ورودیهای این تابع تعداد المانها (e)، تعداد گرهها (n)، اطلاعات ورودی المانها (Elem)، زاویه المانها (ElemAlfa) و طول المانها (L) میباشند. متن این تابع در ادامه مشاهده میشود. خروجی این تابع در متغیر KTotal ذخیره میشود.
functionKTotal= ElemStiffAssamblage(n,e,Elem,ElemAlfa,L)
KTotal= zeros(2*n,2*n)
for m=1:e;
kLoc=((Elem(m,4)*Elem(m,5))/L(m))*[1 0 -1 0; 0 0 0 0; -1 0 1 0; 0 0 0 0];
s= sind(ElemAlfa(m));
c= cosd(ElemAlfa(m));
T= [c s 0 0; -s c 0 0; 0 0 c s; 0 0 -s c]; %rotation tensor
KElem= inv(T)*kloc*T;
i= Elem(m,2);
j= Elem(m,3);
KTotal(2*i-1:2*i, 2*i-1:2*i)= KTotal(2*i-1:2*i, 2*i-1:2*i)+ kelm(1:2,1:2);
KTotal(2*i-1:2*i, 2*j-1:2*j)= KTotal(2*i-1:2*i, 2*j-1:2*j)+ kelm(1:2,3:4);
KTotal(2*j-1:2*j, 2*i-1:2*i)= KTotal(2*j-1:2*j, 2*i-1:2*i)+ kelm(3:4,1:2);
KTotal(2*j-1:2*j, 2*j-1:2*j)= KTotal(2*j-1:2*j, 2*j-1:2*j)+ kelm(3:4,3:4);
% stiffness matrix assamblage
end
end
3-3-4- تابع اعمال شرایط مرزی:
برای اعمال شرایط مرزی در این برنامه از دو تابع ForceBoundary و KBoundary استفاده شده است. این دو تابع شرایط مرزی مسئله را طبق توضیحات ارائه شده در قسمت 2-5 اعمال میکنند. ForceBoundary شرایط مرزی مسئله را بر روی بردار نیرو و KBoundary شرایط مرزی مسئله را بر روی ماتریس سختی اعمال میکند.
ورودیهای تابع ForceBoundary، تعداد گرهها (n)، ماتریس سختی کل(KTotal)، اطلاعات شرایط مرزیجابجایی و نیرو (BCD و BCF)، تعداد گرههایی که شرایط مرزی جابجایی دارند (NDBC) و تعداد گرههایی که شرط مرزی نیرو دارند (NFBC) میباشند. متن تابع ForceBoundary در ادامه آورده شده است.
function F= ForceBoundary(n,BCD,BCF,NDBC,NFBC,KTotal)
F=zeros(2*n,1);
fori=1:NDBC;
if BCD(i,2)~= 99999;
for k= 1:2*n;
F(k)= F(k)- KTotal(k,2*BCD(i,1)-1)*BCD(i,2);
end
end
if BCD(i,3)~= 99999;
for k= 1:2*n;
F(k)= F(k)- KTotal(k,2*BCD(i,1))*BCD(i,3);
end
end
end
if NFBC~=0;
fori=1:NFBC;
F(2*BCF(i,1)-1)= F(2*BCF(i,1)-1)+ BCF(i,2);
F(2*BCF(i,1))= F(2*BCF(i,1))+ BCF(i,3);
end
end
fori=1:NDBC;
if BCD(i,2)~= 99999;
F(2*BCD(i,1)-1)= BCD(i,2);
end
if BCD(i,3)~= 99999;
F(2*BCD(i,1))= BCD(i,3);
end
end
ورودیهای تابع KBoundary، اطلاعات شرایط مرزی جابجایی (BCD)، تعداد گرههایی که شرایط مرزی جابجایی دارند (NDBC)) و ماتریس سختی کل (KTotal) میباشند. متن تابع KBoundary در ادامه مشاهده میشود.
function K= KBoundary(KTotal,BCD,NDBC)
K_Temp= KTotal;
fori= 1:NDBC;
if BCD(i,2)~= 99999;
K_Temp(2*BCD(i,1)-1,:)= 0;
K_Temp(:,2*BCD(i,1)-1)= 0;
K_Temp(2*BCD(i,1)-1,2*BCD(i,1)-1)= 1;
end
if BCD(i,3)~= 99999;
K_Temp(2*BCD(i,1),:)= 0;
K_Temp(:,2*BCD(i,1))= 0;
K_Temp(2*BCD(i,1),2*BCD(i,1))= 1;
end
end
K= K_Temp;
4-3- مرحله پس پردازش و محاسبه سایر خروجیهای مورد نیاز:
در این مرحله از برنامه با استفاده از مقادیر گرهای جابجاییها، سایر خروجیهای مورد نیاز مسئله محاسبه میشوند.
1-4-3- تابع محاسبه تنش:
با استفاده از تابع ElemStress، تنشهای المانهای خرپا محاسبه میشوند.
function Stress= ElemStress(e,a,ElemAlfa,L,Elem)
Stress= zeros(e,1);
for k= 1:e;
s= sind(ElemAlfa(k));
c= cosd(ElemAlfa(k));
T= [c s 0 0; -s c 0 0; 0 0 c s; 0 0 -s c];
i= Elem(k,2);
j= Elem(k,3);
xlocal= T* [a(2*i-1:2*i); a(2*j-1:2*j)];
% displacement in local coordinate.
Stress(k)= (Elem(k,4)/L(k))* (xlocal(3)- xlocal(1));
end
;
2-4-3- تابع محاسبه کرنش:
با استفاده از تابع ElemStrain، کرنشها در المانهای خرپا محاسبه میشوند.
function Strain=ElemStrain(e,a,ElemAlfa,L,Elem
Strain=zeros(e,1);
for m= 1:e;
s= sind(ElemAlfa(m));
c= cosd(ElemAlfa(m));
T= [c s 0 0; -s c 0 0; 0 0 c s; 0 0 -s c];
i= Elem(m,2);
j= Elem(m,3);
xlocal= T* [a(2*i-1:2*i); a(2*j-1:2*j)];
Strain(m)= (1/ L(m))* (xlocal(3)- xlocal(1));
End
3-4-3- تابع محاسبه نیروهای داخلی:
با استفاده از تابع ElemForce، نیروهای داخلی المانهای خرپا محاسبه می شوند.
function Force= ElemForce(e,a,alf,l,elm)
Force=zeros(e,1);
for m= 1:e;
s= sind(alf(m));
c= cosd(alf(m));
T= [c s 0 0; -s c 0 0; 0 0 c s; 0 0 -s c];
i= elm(m,2);
j= elm(m,3);
xlocal= T* [a(2*i-1:2*i); a(2*j-1:2*j)];
Force(m)= ((elm(m,4)* elm(m,5))/l(m))* (xlocal(3)- xlocal(1));
end
4-4-3- تابع محاسبه نیروهای تکیه گاهی:
با استفاده از تابع ReactionForce، نیروهای تکیه گاهی محاسبه میشود.
functionReacForce= ReactionForce(n,NDBC,BCD,KTotal,a)
ReacForce= zeros(2*n,1);
fori= 1:NDBC;
if BCD(i,2)~= 99999;
ReacForce(2*BCD(i,1)-1)= (KTotal(2*BCD(i,1)-1,:)*a);
end
if BCD(i,3)~= 99999;
ReacForce(2*BCD(i,1))= (KTotal(2*BCD(i,1),:)*a);
end
end
5-3- نتایج تحلیل خرپای دوبعدی با استفاده از برنامه
دادده های ورودی مطابق بخش 3-2، برای خرپای شکل 14 وارد و برنامه Mainbody اجرا شده است. نتایج اجرای برنامه در ادامه گزارش شده است.
شکل 14 خرپای تحلیل شده به وسیله برنامه کامپیوتری
مقادیر گره ای جابجایی ها در دستگاه مختصات کلی در جدول 3 مشاهده می شوند.
در گره های 1 و 10 تکیه گاه ساده وجود دارد که هر دو درجه آزادی گره را مقید می کند بنابراین بدیهی است مقادیر جابجایی در گره 1 و 10 در جدول 3 صفر باشد. همانطور که در شکل 14 مشخص است، خرپا هم از نظر هندسی و هم از نظر بارگذاری متقارن است که این امر در جدول 3 نیز مشاهده می شود.
نیروهای تکیه گاهی در گره های 1 و 10 نیز محاسبه شده و نتیجه در جدول 5 گزارش شده است.
جدول 5 نیروهای تکیه گاهی حاصل از برنامه
مجموع نیروهای تکیه گاهی در راستای Y برابر با KN10 است که برابر است با نیروی خارجی وارد بر خرپا که نشان دهنده ارضاء شدن رابطه تعادل میباشد.
سایر خروجی خواسته شده مانند تنشها، کرنش ها و نیروهای داخلی المانها محاسبه شده و در جدول 6 گزارش شده اند.
جدول 6 مقادیر تنش، کرنش و نیرو در المان ها حاصل از برنامه
همانطور که در جدول 6 ملاحظه می شود، تقارن مسئله در مقادیر تنش ها، کرنش ها و نیروهای داخلی المان ها نیز قابل مشاهده است.