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

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

صورت مساله

تیر زیر را با مشخصات (L, E, I, A) در نظر بگیرید که دو ممان و نیروی گسترده به آن اعمال شده است:

f (x , M0 )

هدف از حل مساله، بدست آوردن تغییر شکل (deflection)، تنش (stress) و کرنش (strain) در تیر در سه حالت زیر با روش FEM می باشد.
الف) حل مساله با روش معادلات دیفرانسیلی (جواب دقیق)
ب) حل با روش Weighted Residual Methods یا همان WRM
ج) رسم نمودارهای جابجایی، تنش و کرنش

f(x) = f₀ [1-(x/l)^0.5] 

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

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

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

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

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

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));

بیشتر بخوانید: دانلود پروژه FEM در MATLAB

ضرایب 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’)

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

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

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

تقریب با استفاده از روش Weighted Residuals Methods

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

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

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

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

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

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

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 را داریم. به راحتی می توانیم با نوشتن تابع باقیمانده در نقاطی بین بازه ای معادلاتی پیدا کنیم که ضرایب مجهول ما را معلوم نمایند.معادلات از این قرارند:

% 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 دارد.

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

با استفاده از سه معادله بالا میتوانیم ضرایب 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 خواهیم شد.

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

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

برای درجه 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))

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

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

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

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

% 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))

رسم نمودارها

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

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