حل یک مساله تیر (Beam) با روشهای FEM

در این پسن قصد داریم یک مساله را با چند روش اجزا محدود (FEM) حل نماییم. در این مساله یک تیر که تحت نیروی گسترده و ممان قرار دارد، در دو نقطه مقید شده است. مساله با روشهای WRM، گلرکین (Galerkin)، روش Collocation point و روش حداقل مربعات جزئی (Least Square) حل خواهد شد.

صورت مساله

تیر زیر را با مشخصات (L, E, I, A) در نظر بگیرید که دو ممان و نیروی گسترده به آن اعمال شده است( (f(x و M0 ). هدف از حل مساله، بدست آوردن تغییر شکل (deflection)، تنش (stress) و کرنش (strain) در تیر در سه حالت زیر می باشد.
الف) حل مساله با روش معادلات دیفرانسیلی (جواب دقیق)
ب) حل با روش Weighted Residual Methods یا همان WRM
ج) رسم نمودارهای جابجایی، تنش و کرنش

f(x) = f0 [1-(x/l)0.5] otgg

beam

حل یک مساله تیر (Beam) با روشهای FEM

برای نگارش حل دقیق با استفاده از دیفرانسیل مرتبه چهارم می توانیم تابع گشتاور و شیب و خیز را در طول تیر بدست بیاوریم:

beam equation

با توجه به روابط بالا و شرایط مرزی که در تیر موجود است میتوانیم دو ثابت C3 و C4 را اینگونه بدست آوریم.

با توجه به معادله شیب و خیز در نقطه 0 = نتیجه می شود که C3 و C4 صفر بوده پس معادله ما به دو ضریب مجهول C1 و C2 تقلیل پیدا خواهد کرد. برای حل این دو مجهول دیگر با استفاده از دو شرط مرزی دیگر که در نقطه x=l تعریف می شوند میتوانیم با حل دو معادله دومجهول ضرائب C1 و C2 را بدست بیاوریم.

beam

پس از بدست آوردن این دو ضریب بر حسب پارامترهای دیگر می توانیم در نهایت یک تابع به عنوان تابع خیز در طول تیر داشته باشیم.که همان حل دقیق مسئله ماست. برای حل از نرم افزار متلب کمک گرفته و خواهیم داشت:

clear,clc
% Define some Char
syms L E I F0 yex M Ml x C1 C2 A
% Define Load Function
F=F0*(1-(x/L)^2);
% Derive Moment Function
M=int(int(F))+C1*x+C2;
% Equation from BC for Moment at x = L
equ1=subs(M,x,L)-Ml;
% Deflection function
yex=int(int(M));
% Equation from BC for Deflection at x = L
equ2=subs(yex,x,L);
% Derive C1 & C2 from equ1-equ2
[c1,c2]=solve(equ1,equ2,'C1','C2');
% Deflection Function without unknown C1 & C2
D=(subs(yex,{C1 C2},{c1 c2})/(E*I));

ضرایب C1و C2 معادله خیز به شکل زیر است:

c1 =
(180*Ml - 61*F0*L^2)/(120*L)
c2 =
(11*F0*L^2)/120 - Ml/2
D =
((F0*x^4)/24 - (x^2*(Ml/2 - (11*F0*L^2)/120))/2 - (F0*x^6)/(360*L^2) +
(x^3*(180*Ml - 61*F0*L^2))/(720*L))/(E*I)

حال برای رسم یک نمونه نمودار تغییر مکان نیازمند تعریف یکسری مقادیر اولیه برای پارامترهای مسئله هستیم که به صورت فرضی این مقادیر را اخد می نماییم:

 

% initialize some Vlaue to plot
a=0.1;I1=(a^4)/12;E1=200e+9;
F01=2000;Ml1=130;L1=1;

برای رسم نمودار بایستی مقادیری که مشخص کرده ایم را در معادلات جایگذاری نموده و نقاط رسم را تعیین کنیم.

نمودارهای تنش و کرنش هم با توجه به معادله گشتاور که بدست آورده ایم و رابطه σ=mc/I میتوان مستقیما نتیجه شود. و کرنش هم با توجه به σ=Ee بدست خواهد آمد در نتیجه داریم:

D=subs(D,{I E F0 Ml L},{I1 E1 F01 Ml1 L1});
% Derive Stress and Strain Function
sigma=-(subs(M,{C1 C2},{c1 c2})*a/(2*E*I^2));
sigma=subs(sigma,{I E F0 Ml L},{I1 E1 F01 Ml1 L1});
X=0:0.01:L1;
d=subs(D,x,X);
%-------------------------------
figure
plot(X,d,'r');xlabel('length (m)');ylabel('Deflection (m)');
title('Exact Deflection of the Beam')
sigma=subs(sigma,x,X);
%-------------------------------
figure
plot(X,sigma);xlabel('length (m)');ylabel('Stress');
title('Exact Maximum Stress')
strain=sigma/E1;
%-------------------------------
figure
plot(X,strain);xlabel('length (m)');ylabel('Strain');
title('Exact Maximum Strain')

نتایج حاصل از حل دقیق مسئله در نمودار های روبرو نمایش داده شده است.

beam

مقدار M یا گشتاور در انتهای تیر را بگونه ای گرفته ایم که ممان خمشی حاصل تغییر قابل ملاحظه ای در مقایسه با F کششی گسترده که بر روی تیر اعمال میشود، داشته باشد.

همچنین نمودار های تنش و کرنش را آنچنان که گفته شد، رسم کرده و از آنجا که تنش و کرنش با ضریب E به یکدیگر مربوطند در نتیجه نمودار های حاصل به یک شکل است منتهی اعداد بدست آمده متفاوت می باشد.

beam

تقریب با استفاده از روش Weighted Residuals Methods
در این قسمت ما نیازمند حدس یک تابع تقریب با ضرایب نامعلوم هستیم که بایستی در گام اول شرایط مرزی ما را برآورده کند. در کار حاضر ما اقدام به تخمین تابع خیز توسط یک چند جمله ای مرتبه 5 می نماییم. سپس مرتبه تابع تخمین را یک مرتبه افرایش می دهیم و آنگاه نتایج را با حالت قبل بررسی می کنیم. تابع تقریب yapp به صورت درجه 5 فرض می کنیم فرض می کنیم. اما این حدس به گونه است که شرایط مرزی را برآورده کند.

beam

تابع تقریب درجه 6 هم به صورت زیر است:

beam

حال با توجه به دو شرط مرزی یکی (تغییر مکان صفر و شیب صفر(در ابتدای تیر و یکی) تغییر مکان صفر( در انتهای تیر،برخی از ضرایب چند جمله ای بالا بدین گونه محاسبه می شوند.

beam

شرط مرزی دیگر گشتاور Ml در انتهای تیر است که به صورت نقطه ای وارد میشود و میتوان با معادله که بدست می آید یکی دیگر از این ضرایب را معلوم کرد.

beam

از این معادله میتوان ضریب a2 را بر حسب ضرایب دیگر محاسبه کرد و آنرا از تابع تقریب حذف نمود.در نهایت سه ضریب a3 ,a4,a5 در دو تابع تقریب ما باقی خواهند ماند که با توجه به روش های متفاوت میتوانیم آنها را بدست آوریم. حال تابع باقیمانده را برای هریک از توابع تخمین به صورت جداگانه محاسبه کرده و داریم:

beam

تمامی مراحل بالا را به صورت زیر کد نویسی کرده و داریم:

syms a2 a3 a3 a4 a5
% yapp5 = (L-x)*(a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4)
% yapp6 = (L-x)*(a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4 + a5*x^5)
% From BC's a0 & a1 = 0
yapp5=(L-x)*(a2*x^2+a3*x^3+a4*x^4);
yapp6=(L-x)*(a2*x^2+a3*x^3+a4*x^4+a5*x^5);
a225=solve(subs(E*I*diff(yapp5,2),x,L)-Ml,a2);
a226=solve(subs(E*I*diff(yapp6,2),x,L)-Ml,a2);
yapp5=subs(yapp5,a2,a225);
yapp6=subs(yapp6,a2,a226);
R5=(E*I*diff(yapp5,4)-F);
R6=(E*I*diff(yapp6,4)-F);

نتایج محاسبات بالا بدین ترتیب است:

a225 = -(Ml + E*I*(8*a4*L^3 + 6*a3*L^2))/(4*E*I*L)
a226 = -(Ml + E*I*(10*a5*L^4 + 8*a4*L^3 + 6*a3*L^2))/(4*E*I*L)
yapp5 = (L - x)*(a3*x^3 + a4*x^4 - (x^2*(Ml + E*I*(8*a4*L^3 +
6*a3*L^2)))/(4*E*I*L))
yapp6 = (L - x)*(a3*x^3 + a4*x^4 + a5*x^5 - (x^2*(Ml + E*I*(10*a5*L^4 +
8*a4*L^3 + 6*a3*L^2)))/(4*E*I*L))
R5 = F0*(x^2/L^2 - 1) - E*I*(24*a3 + 96*a4*x - 24*a4*(L - x))
R6 = F0*(x^2/L^2 - 1) - E*I*(24*a3 + 96*a4*x - (24*a4 + 120*a5*x)*(L -
x) + 240*a5*x^2)

 

Collocation point Method

در این روش کافی است بدانیم تابع تخمین ما دارای چند ضریب مجهول است تا بوسیله معادلاتی که بین بازه مورد نظر می نویسیم بتوانیم آن ضرایب را معین نماییم. در حالت تابع تخمین درجه 5، ما دو مجهول داریم a3,a4 و در حالت درجه 6 سه مجهولa3,a4,a5 را داریم. به راحتی می توانیم با نوشتن تابع باقیمانده در نقاطی بین بازه ای معادلاتی پیدا کنیم که ضرایب مجهول ما را معلوم نمایند.معادلات از این قرارند:

beam

% point Collocation Method order 5
[a3p,a4p]=solve(subs(R5,x,(L/5)),subs(R5,x,(3*L/5)),a3,a4);
ypc5=subs(yapp5,{a3 a4},{a3p a4p});
Ypc5=subs(ypc5,{x I E F0 Ml L},{X I1 E1 F01 Ml1 L1});
plot(X,Ypc5,'k+');title('Collocation point Vs. Exact')
% point Collocation Method order 6
[a3p,a4p,a5p]=solve(subs(R6,x,(L/5)),subs(R6,x,(2*L/5)),subs(R6,x,(3*L/5
)),a3,a4,a5);
ypc6=subs(yapp6,{a3 a4 a5},{a3p a4p a5p});
Ypc6=subs(ypc6,{x I E F0 Ml L},{X I1 E1 F01 Ml1 L1});
plot(X,Ypc6,'k+');title('Collocation point Vs. Exact')

نتایج تابع تخمین مرتبه 5:

a3p =
-F0/(25*E*I)
a4p =
F0/(150*E*I*L)
ypc5 =
-(L - x)*((F0*x^3)/(25*E*I) + (x^2*(Ml - (14*F0*L^2)/75))/(4*E*I*L) -
(F0*x^4)/(150*E*I*L))

نتایج تابع تخمین مرتبه 6:

a3p =
-(7*F0)/(180*E*I)
a4p =
F0/(360*E*I*L)
a5p =
F0/(360*E*I*L^2)
ypc6 =
-(L - x)*((7*F0*x^3)/(180*E*I) + (x^2*(Ml - (11*F0*L^2)/60))/(4*E*I*L)
- (F0*x^4)/(360*E*I*L) - (F0*x^5)/(360*E*I*L^2))

نتیجه:
در نمودار زیر دو تابع تخمین به همراه تابع دقیق خیز رسم شده اند که حاکی از دقت بالای مرتبه 6 دارد.

beam

در این روش هم با استفاده از انتگرال گیری از تابع باقیمانده در بازه های بین بازه ای، میتوانیم ضرایب را برای معادله جدید بدست بیاوریم. فرضا برای تابع تخمین درجه 6 داریم:

beam

با استفاده از سه معادله بالا میتوانیم ضرایب a3, a4, a5 را بدست آوریم:

% Subdomain Method order 5
[a3s,a4s]=solve(int(R5,0,(L/4)),int(R5,(L/2),(9*L/10)),a3,a4);
ysub5=subs(yapp5,{a3 a4},{a3s a4s});
Ysub5=subs(ysub5,{x I E F0 Ml L},{X I1 E1 F01 Ml1 L1});
plot(X,Ysub5,'m+');title('Subdomain Method Vs. Exact')
% subdomain Method order 6
[a3s,a4s,a5s]=solve(int(R6,0,(L/5)),int(R6,(3*L/10),(L/2)),int(R6,(L/2),
(7*L/10)),a3,a4,a5);
ysub6=subs(yapp6,{a3 a4 a5},{a3s a4s a5s});
Ysub6=subs(ysub6,{x I E F0 Ml L},{X I1 E1 F01 Ml1 L1});
plot(X,Ysub6,'m+');title('Subdomain Method Vs. Exact')

نتایج محاسبات این روش به شکل زیر است:
درجه 5 :

a3s =
-(3161*F0)/(82800*E*I)
a4s =
(193*F0)/(27600*E*I*L)
ysub5 =
-(L - x)*((3161*F0*x^3)/(82800*E*I) + (x^2*(Ml -
(2389*F0*L^2)/13800))/(4*E*I*L) - (193*F0*x^4)/(27600*E*I*L))

درجه 6 :

a3s =
-(7*F0)/(180*E*I)
a4s =
F0/(360*E*I*L)
a5s =
F0/(360*E*I*L^2)
ysub6 =
-(L - x)*((7*F0*x^3)/(180*E*I) + (x^2*(Ml - (11*F0*L^2)/60))/(4*E*I*L)
- (F0*x^4)/(360*E*I*L) - (F0*x^5)/(360*E*I*L^2))

 

در نهایت توسط رسم این توابع متوجه تفاوت مرتبه 5 و دقت مرتبه 6 خواهیم شد.

beam

حل با استفاده از روش Galerkin’s Method

این روش هم همانند روش قبلی است منتها انتگرال های میان بازه ای به شکل زیر تغییر میکنند که Wi=Qi  می باشد. Qi عباراتی هستند که در کنار ضرایب مجهول ظاهر می شوند. که در این مسئله برابرند با:

beam

برای درجه 5 و درجه 6 :

(-((3*L*x^2)/2 - x^3)*(L - x)) * a3
((x^4 - 2*L^2*x^2)*(L - x)) * a4
% -------------------------
(-((3*L*x^2)/2 - x^3)*(L - x)) * a3
((x^4 - 2*L^2*x^2)*(L - x)) * a4
((x^5 - (5*L^3*x^2)/2)*(L - x)) * a5

حال با توجه به انتگرال بین بازه ای بالا معادلات مختلفی را برای حل در بازه های مختلف بین بازه اصلی تعریف میکنیم و در نهایت داریم:

% Galerkin's Method order 5
[a3g,a4g]=solve(int((-((3*L*x^2)/2 - x^3)*(L -
x))*R5,0,(L/5)),int(((x^4 - 2*L^2*x^2)*(L -
x))*R5,(3*L/10),(L/2)),a3,a4);
ygal5=subs(yapp5,{a3 a4},{a3g a4g});
Ygal5=subs(ygal5,{x I E F0 Ml L},{X I1 E1 F01 Ml1 L1});
plot(X,Ygal5,'g*');

% Galerkin's Method order 6
[a3g,a4g,a5g]=solve(int((-((3*L*x^2)/2 - x^3)*(L -
x))*R6,0,(L/5)),int(((x^4 - 2*L^2*x^2)*(L -
x))*R6,(3*L/10),(L/2)),int(((x^5 - (5*L^3*x^2)/2)*(L -
x))*R6,(3*L/5),(4*L/5)),a3,a4,a5);
ygal6=subs(yapp6,{a3 a4 a5},{a3g a4g a5g});
Ygal6=subs(ygal6,{x I E F0 Ml L},{X I1 E1 F01 Ml1 L1});
plot(X,Ygal6,'g*');

 

a3g = -(60967493593*F0)/(1545120532200*E*I)
a4g = (1035229541*F0)/(220731504600*E*I*L)
ygal5 = -(L - x)*((60967493593*F0*x^3)/(1545120532200*E*I) - (1035229541*F0*x^4)/(220731504600*E*I*L) +
(x^2*(Ml - (153916053631*F0*L^2)/772560266100))/(4*E*I*L))
%-------------
a3g = -(7*F0)/(180*E*I)
a4g = F0/(360*E*I*L)
a5g = F0/(360*E*I*L^2)
ygal6 = -(L - x)*((7*F0*x^3)/(180*E*I) + (x^2*(Ml - (11*F0*L^2)/60))/(4*E*I*L) - (F0*x^4)/(360*E*I*L) -
(F0*x^5)/(360*E*I*L^2))

در نهایت پس از حل معادلات و معلوم کردن ضرایب داریم:

beam

حل با استفاده از روش Least Square

در این روش هم همانند Galerkin بایستی تابع وزن در باقیمانده درون انتگرال ضرب شود منتها این تابع مشتق از خود باقیمانده نسبت به ضرایب مجهول است. به عبارتی دیگر معادلات ما توسط این انتگرال ساخته می شوند:

beam

در نتیجه خواهیم داشت:

% Least Square Method order 5
[a3l,a4l]=solve(int(diff(R5,a3)*R5,0,(L/5)),int(diff(R5,a4)*R5,(3*L/10),
(L/2)),a3,a4);
yls5=subs(yapp5,{a3 a4},{a3l a4l});
Yls5=subs(yls5,{x I E F0 Ml L},{X I1 E1 F01 Ml1 L1});
plot(X,Yls5,'k*');
% Least Square Method order 6
[a3l,a4l,a5l]=solve(int(diff(R6,a3)*R6,0,(L/5)),int(diff(R6,a4)*R6,(3*L/
10),(L/2)),int(diff(R6,a5)*R6,(3*L/5),(9*L/10)),a3,a4,a5);
yls6=subs(yapp6,{a3 a4 a5},{a3l a4l a5l});
Yls6=subs(yls6,{x I E F0 Ml L},{X I1 E1 F01 Ml1 L1});
plot(X,Yls6,'k*');

 

a3l = -(533*F0)/(13680*E*I)
a4l = (49*F0)/(11400*E*I*L)
yls5 = -(L - x)*((533*F0*x^3)/(13680*E*I) + (x^2*(Ml -
(2273*F0*L^2)/11400))/(4*E*I*L) - (49*F0*x^4)/(11400*E*I*L))
% ===============
a3l =-(7*F0)/(180*E*I)
a4l =F0/(360*E*I*L)
a5l = F0/(360*E*I*L^2)
yls6 =-(L - x)*((7*F0*x^3)/(180*E*I) + (x^2*(Ml -
(11*F0*L^2)/60))/(4*E*I*L) - (F0*x^4)/(360*E*I*L) -
(F0*x^5)/(360*E*I*L^2))

و در نهایت نتایج را رسم می نماییم:

beam

 

در نهایت در یک نگاه کلی برای تمامی تقریب های مرتبه 5 داریم:

beam

مطالب مرتبط

دیدگاهتان را بنویسید

بخش های مورد نیاز علامت گذاری شده اند

نشانی ایمیل منتشر نخواهد شد

نویسنده : آدرس سایت : ایمیل :


0

شبکه های اجتماعی

دانشنامه تخصصی مهندسی ایران را در شبکه های اجتماعی دنبال کنید

0 0

عضویت در خبرنامه

برای دریافت آخرین اخبار در زمینه مهندسی شامل نرم افزارها، استانداردها و آموزش ها به سامانه اطلاع رسانی ما بپیوندید.

بدون پرداخت هزینه، تا هر وقت بخواهید.

تست

همکاران ما

گروه مپنا
گروه مپنا
دانشگاه تهران
دانشگاه تهران
سایپا
سایپا
ایران خودرو
ایران خودرو
شرکت ملی نفت ایران
شرکت ملی نفت ایران
ذوب‌آهن اصفهان
ذوب‌آهن اصفهان
فولاد خوزستان
فولاد خوزستان
مشاوره

نیاز به مشاوره دارید؟

 
                    همکاران ما پاسخگو شما خواهند بود.