ҚР Ұ 218-149-2018 ЖОЛ ТӨСЕМЕСІНІҢ КӨП ҚАБАТТЫ ҚҰРЫЛЫМЫНЫҢ КЕРНЕУЛІ-ДЕФОРМАЦИЯЛАНҒАН КҮЙІН АНЫҚТАУ БАҒДАРЛАМАСЫН ТӘЖІРИБЕЛІК ҚОЛДАНУ БОЙЫНША ӘДІСТЕМЕЛІК ҰСЫНЫМДАР

Қазақстан Республикасы Инвестициялар және даму министрлігі Автомобиль жолдары комитеті Төрағасының 2018 жылғы 21 желтоқсандағы № 122 бұйрығымен бекітілген.

МАЗМҰНЫ

Алғысөз

1	"Қазақстан жол ғылыми-зерттеу институты" акционерлік қоғамы ("ҚазжолҒЗИ" АҚ) ДАЙЫНДАП ЕНГІЗДІ
2	Қазақстан Республикасы Инвестициялар және даму министрлігі Автомобиль жолдары комитеті Төрағасының 2018 жылғы 21 желтоқсандағы № 122 бұйрығымен БЕКІТІЛІП, ҚОЛДАНЫСҚА ЕНГІЗІЛДІ
3	"ҚазАвтоЖол" ҰК" Акционерлік қоғамымен 14.11. 2018 жылғы № 03/14-2-2623-И хатымен КЕЛІСІЛДІ
4	БІРІНШІ ТЕКСЕРУ МЕРЗІМІ	2023 жыл
5	ТЕКСЕРУ КЕЗЕҢДІЛІГІ	5 жыл
6	АЛҒАШ РЕТ ЕНГІЗІЛДІ

Кіріспе

Бұл әдістемелік ұсынымдар Б01.02 "Жол төсемесінің көп қабатты құрылымының кернеулі-деформацияланған күйін анықтауға арналған теориялық шешімдер мен бағдарламаларды әзірлеу" тақырыбының шеңберінде қабылданған жұмыстар жоспарына сәйкес әзірленді.

Құжат автомобиль салмағымен статикалық жүктелген көп қабатты жол құрылымының кернеулі-деформацияланған күйін анықтаудың үлгілік міндет қойылымдарынан тұрады және негізгі формулалар шекті элементтер әдісімен есептеу алгоритмінен алынған.Әдісте сегіз түйінді дәлдігі жоғары квадратты тікбұрышты шекті элемент қолданылады. А қосымшасында қажетті түсіндірмелері бар MATLAB [1-4] тіліндегі BASIC_NDS_MKE_8_uzlov есептеу бағдарламасының мәтіні (листинг) берілген. Б қосымшасында элементтің қаттылық матрицасын есептеуге арналған BDB0 бағдарламасының листингі келтірілген. BASIC_NDS_MKE_8_uzlov бағдарламасын қолданудың қысқаша нұқсаулығы да берілген. Құжаттың соңында қолданылған әдебиеттер тізімі көрсетілген.

1 Қолдану саласы

1.1 Қазақстан Республикасының жалпы пайдаланымдағы автомобиль жолдарының желісіне таралады және көп қабатты жалпы пайдаланымдағы автомобиль жолдарын жобалау мәселелерін шешуге арналған.

1.2 Жалпы пайдаланымдағы автомобиль жолдарына арналған жол төсемелерінің құрылымдарын жобалауда, жобалау және пайдалану сатысында жол төсемелерін есептеуде (ҚР ЕЖ 3.03-103-2014 және ҚР ЕЖ 3.03-104-2014 сәйкес), сондай-ақ автомобиль жолдарына қатысты инженерлік-экономикалық міндеттерді шешуде осы ұсынымдарды басшылыққа алу қажет.

2 Нормативтік сілтемелер

Осы ұсынымдарды қолдану үшін келесі сілтемелік нормативтік құжаттар қажет:

ҚР ЕЖ 3.03-103-2014 "Қатты жол төсемдерін жобалау"

ҚР ЕЖ 3.03-104-2014 "Қатты емес типті жол төсемдерін жобалау"

Ескерту – Осы ұсынымдарды қолдану кезінде ағымдағы жылдың күйі бойынша жасалған "Стандарттау бойынша нормативтік құжаттар" көрсеткіші бойынша, және де ағымдағы жылда жарияланған тиісті ақпараттық көрсеткіштер бойынша сілтемелік стандарттардың әрекет етуін тексеру қажет. Егер сілтемелік құжат ауыстырылған (өзгертілген) болса, онда осы ұсынымдарды пайдалану кезінде ауыстырылған (өзгертілген) стандартты нұсқау етіп қолдану қажет. Егер сілтемелік құжат ауыстырылмай күші жойылған болса, онда оған сілтеме берілген ереже осы сілтемені қозғамай қатысты қолданылады.

3 Терминдер мен анықтамалар

Осы Ұсынымдарда тиісті анықтамалары бар келесі терминдер қолданылады:

3.1 Жол төсемесі: Автокөлік жүруге және оның салмағының әсерін жер төсемесіне беру үшін жасалған жүру бөлігінің көп қабатты құрылымы.

3.2 Жер төсемесі: Жол төсемесін, сондай-ақ жол қозғалысын ұйымдастырудың техникалық құралдарын орналастыруға және автомобиль жолдарын жайластыру үшін негіз болып қызмет атқаратын құрылымдық элемент.

3.3 Деформация: Дененің сызықтық өлшемдерінің бастапқы өлшемдеріне қарағанда өзгеруімен сипатталатын салыстырмалы шама.

3.4 Кернеу: алаң бірлігіне келетін қалыпты немесе жанама жүктемемен анықталынатын салыстырмалы шама.

3.5 Қатқыл жол төсемесі: Асфальтбетондардың әртүрлі түрлерінен, битуммен, цементпен, әктаспен, кешенді және басқа да тұтқырғыштармен, сондай-ақ әлсіз байланысты түйірлі материалдармен нығайтылған материалдардан салынған қабатты жол төсемесі.

3.6 Құрылымдық қабат: Біркелкі материалдардан тұратын және іргелес қабаттардардан материал түрімен, беріктігімен және құрамымен ерекшелінетін жол төсемесінің әрбір қабаты.

3.7 Жол жамылғысы: Жол негізіне салынатын, көлік құралдарынан жүктемені тура қабылдайтын және белгіленген пайдалану талаптарын қамтамасыз етуге арналған және жол негізін атмосфералық факторлардың ықпалынан қорғайтын жол төсемесінің бір немесе көп қабатты жоғарғы бөлігі.

3.8 Жол негізі: Жамылғымен бірге көлік құралдарының жүктемесін қабылдайтын және оны қосымша қабаттарға немесе тікелей жер төсемесінің топырағына үлестіруге арналған жол төсемесінің салмақ түсетін төменгі қабаты.

3.9 Шекті элементтер әдісі: математикалық физиканың дифференциалдық теңдеулерін шешудің сандық әдісі.

4 BASIC_NDS_MKE_8_UZLOV бағдарламасын пайдалану әдістемесі

4.1 Міндеттің жалпы қойылымы

Жалпы пайдаланымдағы автомобиль жолының математикалық үлгісін құру үшін екі қабатты асфальтбетон жол төсемесінің құрылымы мен 8% цемент қоспасы, шағыл тасты қоспа және қиыршық тасты-құмды қоспасы бар шағыл тасты-құмды қоспадан тұратын жол төсемесінің негізін қарастырамыз. Жол төсемесінің негізі жеңіл топырақты негізге орналастырылады (1-сурет).

Көпқабатты құрылымның тепетеңдігі туралы есепті шешуге арналған сандық әдістің - шекті элементтер әдісінің мүмкіндіктері құрылымдық қабаттардың геометриялық өлшемдері мен құрылымда қолданылатын материалдардың физика-механикалық қасиеттерін айтарлықтай жеңіл тағайындауға жағдай жасайды.

Сондықтан 1-кестеде келтірілген математикалық үлгіде қарастырылған құрылымдық қабаттардың геометриялық өлшемдері (1-сурет) мен құрылымдық элементтердің матриалдарының физика-механикалық қасиеттері (оның ішінде серпімділік модульдері) [6] сәйкес тағайындалған.

1 – жамылғы, ұсақ түйірлі тығыз асфальтбетон; 2 – түйірлері ірі кеуек асфальтбетон; 3 – 8% шағыл тасты-құмды қоспа; 4 – шағыл тасты қоспа; 5 – қиыршық тасты-құмды қоспа; 6 – жеңіл балшық топырақ

1-сурет - Жол құрылымының сұлбалық көрінісі

Құрылымдық элемент материалдарының физика-механикалық қасиеттерінің мәндері А қосымшасының 108-110 жолдарында келтірілген.

1-кесте – Құрылымдық қабаттардың сипаттамасы

Құрылымдық қабат материалы	Серпімділік модулі Е, МПа	Пуассон коэффициенті,	Қабаттың биіктігі, м
1 Ұсақ түйірлі тығыз асфальтбетон	3200	0,18	0,05
2 Ірі түйірлері кеуек асфальтбетон	2000	0,27	0,10
3 8% шағыл тасты-құмды қоспа	800	0,30	0,40
4 Шағыл тасты қоспа	275	0,30	0,45
5 Қиыршық тасты-құмды қоспа	180	0,30	1,00
6 Жеңіл балшық топырақ	50	0,35	-

Шекті элементтер әдісімен жол құрылымының кернеулі-деформациялық күйі туралы есепті шешу алгоритмін компьютерлік жүзеге асыру үшін MATLAB жүйесінің алгоритм тіліндегі BASIC_NDS_MKE_8_uzlov есептеу бағдарламасы әзірленді.

Зерттелетін сала сегіз түйінді квадраттық тікбұрышты шекті элементтерге бөлінеді (2-сурет). Мұндай соңғы элементтегі түйіндердің жергілікті нөмірі 1-ден бастап 8-ге дейін сол жақтан төменгі түйіннен басталады және сағат тіліне қарама-қарсы бағытта жүзеге асырылады.

2-сурет – Сегіз түйінді тікбұрышты элемент

Шекті-элемент торын құру үшін және координаттық өстері бағытындағы және айнымалы адымдар қолданылды (3-сурет).

Әрбір көлденең қатарда (3-сурет) 60 элементтен орналасқан (egor=60), ал элементтердің көлденең қатарларының саны 23 (ever=23). Осыдан, зерттелетін аймақ 1380 шекті элементтерге бөлінген, ал түйіндер саны np=(3*ever+2)*egor+2*ever+1=4307 тең (А қосымшасының 12-15 жолдары).

3-сурет – Зерттелетін аймақты шекті-элементтік бөлудің сұлбалық көрінісі

Түйіндердің жалпыланған нөмірлері түйіндердің шеткі сол жақтағы тік қатарынан басталады және жоғарыдан төменге қарап жүзеге асырылады және солдан оңға қарай жылжыйды (3-сурет).

Сонымен қатар, бөліну тұрақсыз болады, себебі сегіз түйінді элементтерде орталық түйін болмайды. Бұл жағдайда тік қатарда 47 түйін, ал келесісінде – 24 түйін болады.

Автомобильдің бір дөңгелегі арқылы түйіндерде тік жүктеме түрінде берілетін автомобиль салмағының ықпалындағы шекті-элементтік тордағы түйіндік нүктелердің жылжу векторының компоненттерін анықтау қажет (4-сурет).

L және H – зерттелетін аймақтың ені мен биіктігі; Р – көлік құралының дөңгелегінен түсетін жүктеме; u_i, v_i - х және у координаттық өстер бойымен бағытталған шекаралардағы нүктелердің жылжу компоненттері; I – VI – жол төсемесі мен топырақты негіздің құрылымдық қабаттары.

4-сурет – Есептің есептік сұлбасы

4.2 Координаттар ауқымдарын қалыптастыру

BASIC_NDS_MKE_8_uzlov бағдарламасында dx және dy ауқымдары берілген, олар координаттар бойынша адым мәндерінен тұрады. Мұндай ауқымдардың болуы координаттар бойынша айнымалы адымдарды беруге мүмкіндік береді (А қосымшасының 44-63 жолдары).

Бағдарламаның координаттарды қалыптастыру бөлігінде (А қосымшасының 68-91 жолдары) ағымдағы шекті элемент түйіндерінің жалпыланған нөмірі 1,2,..., 8 жергілікті нөмірлерге сәйкес келетін n1, n2, …, n8 арқылы белгіленген, 2-суретте келтірілген. BASIC_NDS_MKE_8_uzlov бағдарламасында n1, n2, …, n8 есептеу алгоритмдері бойынша есептелінетін түйіндердің жалпыланған нөмірлерінің нақты мәндері түйіндердің ағымдағы көлденең қатарларындағы элемент нөмірімен (m=1:egor) және түйіндердің ағымдағы тік қатарындағы нөмірімен (n=1:ever) анықталады.

4.3 Шекаралық шарттар

4-суретте зерттелетін аймақтың бүйірлерінде және төменгі шекараларында жылжу векторының тиісті компоненттерінің нөлдік мәндері түрінде беріледі (u=0 и u=v=0). Кейбір түйіндердің шекаралық шарттарын белгілі жылжулар түрінде беру үшін жылжулары белгілі түйіндерінің нөмірлерінен тұратын арнайы Mz ауқымын құру қажет (А қосымшасының 23-43 жолдары).

Жүктемелер түріндегі шекаралық шарттар бір дөңгелекке келетін автомобиль салмағымен анықталатын қадалған тік күштер түрінде беріледі. Жүктемелердегі шекаралық шарттарды есептеу әдістемесі А қосымшасында келтірілген (112-166 жолдар). Есептеудің қорытынды нәтижесі жүйенің (1) тепе-теңдік теңдеулерінің оң жақ бөлігінің векторы арқылы анықталады.

4.4 Тепе-теңдік теңдеуін анықтайтын ара қатынастар

Сызықтық үшбұрышты шекті элементтерді қолданған секілді квадратты элементтерді қолданғанда да жүйенің тепе-теңдік теңдеуі матрицалық пішінде келесіндей көрініске ие [7]:

(1)

мұнда [K] –жүйенің қаттылығының (2*np,2*np) өлшемді матрицасы; – түйіндік жылжулардың 2*np өлшемді векторы;

–түйіндік жүктемелердің 2*np өлшемді векторы;

np – шекті-элементтік тордағы түйіндердің жалпы саны.

Жүктемелердегі шекаралық шарттарға сәйкес келетін векторының компоненттерінің мәндерін есептеу тәртібі А қосымшасының 150-165 жолдарында келтірілген. Бұл вектордың қалған компоненттері түйіндердің тепе-теңдігі жағдайларына қарай нөлге тең.

(1) теңдеуіндегі жүйенің қаттылық матрицасы мынадай жиынтық түрінде анықталады

мұнда – шекті-элементтік тордағы элементтердің жалпы саны; –ағымдағы "" элементінің қаттылық матрицасы.

Элементтің қаттылық матрицасы келесі көлемдік интегралдың көмегімен есептелінеді

(2)

мұнда серпімділік матрицасын есептеу алгоритмі А қосымшасының 204-207 жолдарында келітірілген.

(2) өрнекте градиенттер матрицасы мүсін функциясынан координаттары бойынша дифференциялау арқылы анықталады.

Мүсін функциясының сегіз түйінді квадратты элементтерде орын алатын координаттарынан бейсызықты тәуелділік жағдайында мүсін функциясыны координаттар бойынша алынған туындылары да координаттарына тәуелді болады және интегралын есептеу орасан зор математикалық амалдарды қажет ететін өзіндік есепке айналады.

4.5 Сегіз түйіні бар төртбұрышты элементтің градиенттер матрицасын есептеу

(2) ара қатынасында көлемді интегралды сандық әдіспен есептеу үшін координаттардың жергілікті жүйесіне өту ұсынылады [7] (5-сурет). Мұнда интегралдаудың айнымалыларын ауыстыруды төмендегі ара қатынас арқылы жүзеге асыруға болады:

(3)

мұнда – элементтің бірлік қалыңдығы, ал – Якобидың координаттарды түрлендіру матрицасының анықтауышы. Осылайша (2) қатынасындағы көлемді интеграл төмендегідей түрге ие болады:

(4)

Интегралдаудың және айнымалыларына өту интегралдаудың шектерін анықтауды жеңілдетеді және Гаусс-Лежандр квадратурасы арқылы интегралдарын есептеудің біріңғай алгоритмін әзірлеуге мүмкіндік береді.

5-сурет - Сызықтық тікбұрышты элементке арналған координаттар жүйесі

4.6 Деформациялар мен кернеулерді анықтау алгоритмі

Өрнекті келесі түрде жазамыз (5):

(8) өрнегінде және координаттары сызықтық тікбұрышқа жатады (5-сурет) және оларға сегізбұрышты тікбұрыштардағы және нүктелері сәйкес келетін болады.

Деформациялар тензорының қалған компоненттері де осыған ұқсас анықталады. (5) формуладан компоненті үшін төмендегіні аламыз:

немесе, (4) және екенін ескере отырып, келесіні аламыз:

Деформациялар векторының компоненттерін және кернеулер векторын есептеудің келтірілген алгоритмін сандық орындау бағдарламасы осы құжаттың А қосымшасының 412-527 жолдарында келтірілген. Одан басқа, А қосымшасында тік қималарға тән нүктелердің көлденең жылжуларының (А қосымшасының 297-312 жолдары), түрлі тереңдіктерде орналасқан көлденең қималар нүктелерінің тік жылжуларының (А қосымшасының 316-408 жолдары) және деформациялар мен кернеулер векторлары компоненттері мәндерінің (А қосымшасының 529-602 жолдары) кестелерін түзу бағдарламалары берілген.

4.7 BASIC_NDS_MKE_8_UZLOV бағдарламасын пайдаланудың қысқаша әдістемесі

Келтірілген нормативтік құжаттың мақсаты MATLAB тілінде есептеу бағдарламасын жасау болып табылады.

Есптеу бағдарламасы зерттеудегі жол құрылымының нүктелеріндегі (4-сурет) жылжулар, деформациялар мен кернеулердің компоненттерінің мәндерін анықтауға арналған.

Әлемдік тәжірибеде математикалық физика есептерін, соның ішінде деформацияланатын қатты денелердің керенеулі деформацияланған күйіне арналған есептерді шығаруға қабілетті басқа да бағдарламалар жиынтықтары бар. Олардың ішінде ең белгілісі ANSYS бағдарламар кешені. Бірақ олардың бәрі коммерциялық бағдарламалар болғандықтан олардың алгоритмдері жасырын болады да, қарапайым қолданушыларға оларды өз есептеріне пайдалану, оларды қажетінше өзгерту мүмкін емес. Ұсынылып отырған BASIC_NDS_MKE_8_uzlov бағдарламасы ашық, оның алгоритмі толық ашылып жазылған, бағдарлама жылдам бейімделеді және оны әрқашанда жақсартып отыруға мүмкіндік бар. Мысалы, оған жеке блоктар қосуға болады және жол жамылғысындағы төменгі температуралық жарықшақтардың пайда болу шарттары туралы, жол төсемінің құрылым қабаттарындағы шаршау жарықшақтарының пайда баолу шарттары туралы есептерді шығаруға болады.

Нақты бір есепті шығаратын кезде BASIC_NDS_MKE_8_uzlov бағдарламасын іске қосу үшін директориясында BASIC_NDS_MKE_8_uzlov бағдарламасының мәтіні (листинг) орналасқан MATLAB бағдарламалар кешенінің бір үлгісі орналасуы шарт. Компьютердің жұмыс столында MATLAB бағдарламалар кешенінің таңбасы орналастырылады (6-сурет).

6-сурет – MATLAB бағдарламар кешенінің таңбасы

Осы құжаттың А қосымшасында BASIC_NDS_MKE_8_uzlov бағламасының листингі келтірілген. Мәтінді қолдану ыңғайлы болу үшін листингтің қатарлары нөмірленген, және оны MATLAB бағдарламалар кешенінің директориясына отырғызу алдында алып тастау керек.

MATLAB бағдарламалар кешенін қолдануға шақыру жарлықты екі рет шерту арқылы жасаланады (6-сурет ).

Жұмыс столында MATLAB – тың басты беті пайда болады (7-сурет).

7-сурет – MATLAB бағдарламасының басты бетінің көрінісі

MATLAB-тың 7-суреттегі басты бетінде BASIC_NDS_MKE_8_uzlov бағдарламасының бір бөлігі келтірілген. Егер бұл бөлік көрінбесе, ал басты бетте басқа бағдарламаның бөлігі көрініп тұрса, онда BASIC_NDS_MKE_8_uzlov бағдарламасын іздеу үшін басты беттің жоғарғы жағында OPEN түймесін басу керек. OPEN түймесін басқаннан кейін MATLAB диреториясындағы бағдарламалар тізімі шығады (8-сурет).

8-сурет - Директориядағы бағдарламалар тізімінің көрінісі

BASIC_NDS_MKE_8_uzlov бағдарламасын басты бетке шақыру үшін тізімдегі оның атын шертсе болды (7-сурет).

9-сурет – Қабаттардың материалдарының физика-механикалық сипаттамаларын енгізуге қажет 109-110 қатарлар

Келтірілген бағдарламада автомобиль жолының физика-механикалық сипаттамалары қойылған. Оларды қажетінше өзгерту үшін 109-110 қатарлардағы олардың мәндерін енгізу керек (9-сурет). Физика-механикалық сипаттамалар енгізілгеннен кейін жоғарғы көлбеу бетке әсер етуші тік күштің мәнін беруді 150 қатарда жүргіземіз және олардың әсер ету нүктелерінің нөмірлерін береміз (10-сурет).

10-сурет – Бетке әсер етуші күштің шамасы мен ол әсер етуші түйіннің нөмірін беру қатары

Есепті шығарудың келесі қадамында қолданушының талабына байланысты графикалық нәтижелер алу керек. Төменде тұйіндердің тік бағытта әртүрлі тереңдіктердегі жылжуларын анықтау тапсырмасын ұйымдастыру мысалдары келтірілген (11-17 - суреттер).

11-сурет – Түйіндердің h= 5 см тереңдіктегі тік жылжулар кестесін алу

12-сурет – Түйіндердің h= 15 см тереңдіктегі тік жылжулар кестесін алу

13-сурет – Түйіндердің h= 25 см тереңдіктегі тік жылжулар кестесін алу

14-сурет– Түйіндердің h= 45 см тереңдіктегі тік жылжулар кестесін алу

15-сурет – Кернеулердің SigmaX компонентінің кестесін алу

16-сурет – Кернеулердің SigmaY компонентінің кестесін алу

17-сурет – Кернеулердің TauXY компонентінің кестесін алу

Бағдарламаның басты беттегі барлық мәтінін қарап шығып, қажетті бастапқы мәндерді енгізіп, қажет болса бағдарламаға өзгертулер енгізгеннен кейін оны есептеуге қосуға болады. Ол үшін Run түймесін шертсе болды (18-сурет).

18-сурет– Құралдар бетіндегі бағдарламаны іске қосудың Run түймесі

Бағдарламада тапсырушының қалауы бойынша есепті шығару нәтижелері туралы қажетті мәліметті алу командалары бар. Мұндай мәлімет кесте түрінде, кесте түрінде, немесе мәтін түрінде т.б. болуы мүмкін. Бұл тапсырушының қалауына байланысты.

Төменде, 13-18 - суреттерде, мысал ретінде шешілген есептің нәтижелері зерттеу аймағының қажетті қималарының нүктелеріндегі тік және көлбеу жылжулар өрістерінің қалыптасу суреттері, және кернеулердің компоненттерінің де сол қималардағы өзгеру суреттері келтірілген. Мысал есепте зерттеудегі аймақтың жоғарғы беті автомобильдің жалғыз дөңгелегінің тік бағыттағы әсеріне ұшыраған.

4.8 BASIC_NDS_MKE_8_uzlov бағдарламасы бойынша есептеу үлгілері

Төменде BASIC_NDS_MKE_8_uzlov бағдарласын компьютерде жүзеге асыру нәтижелері жылжулар мен кернеулер кестелері түрінде келтірілген.

1 – тік симметрия өсінің сол жағынан 30 см қашықтықта орналасқан тік қима нүктелерінің көлденең ығысуы; 2 – тік симметрия өсі нүктелерінің көлденең ығысуы; 3 – тік симеттерия өсінің оң жағынан 30 см қашықтықта орналасқан тік қима нүктелерінің көлденең ығысуы

19-сурет – Тік қима нүктелерінің көлденең ығысуы

20-сурет - 1– h=5 см; 2–h=15 см; 3–h=25 см; 4–h=45 см; 5–h=80 см тереңдіктерінде орналасқан көлденең қима нүктелерінің тік ығысуы

21-сурет - h тереңдіктерде кернеулердің үлесуі

22-сурет - h тереңдіктерде кернеулердің үлесуі

23-сурет - h тереңдіктерде кернеулердің үлесуі

А қосымшасы (міндетті) BASIC_NDS_MKE_8_uzlov бағдарламасының бастапқы коды

1	%BASIC_NDS_MKE_8_uzlov
2	% бағдарламысыны базалық нұсқасы
3
4	% СОҢҒЫ ЭЛЕМЕНТТЕР ӘДІСІМЕН 8 ТҮЙІНДІ КВАДРАТТЫ ЭЛЕМЕНТТІҢ
5	% КӨМЕГІМЕН ҚҚС ТУРАЛЫ ТЕСТІЛІК ТАПСЫРМАСЫН % ШЕШУ
6
7	% Әзірлеменің басталуы 8.07.2017 ж
8	% аяқталуы 1.09.2018 ж
9
10	tic
11
12	ever=23; egor=60; % соңғы-элементтік бөліну элементтеріндегі
13	% тік (ever) және
14	% көлденең (egor) қатарлардағы соңғы элементтер саны
15	np=(3*ever+2)*egor+2*ever+1; % түйіндердің жалпы саны
16	x=zeros(np); y=zeros(np);
17
18	% ЭЛЕМЕНТТЕРДІ НӨМІРЛЕУ ТӘРТІБІ – СОЛДАН ОҢҒЫ (ҚАТАРЛАР % ТӨМЕННЕН БАСТАЛАДЫ)
19	% ТҮЙІНДЕРДІ НӨМІРЛЕУ ТӘРТІБІ - СОЛДАН ОҢҒЫ. ТҮЙІНДЕР % ҚАТАРЫ ДА
20	% ТӨМЕННЕН БАСТАЛАДЫ. ТҮЙІНДЕРДІҢ АРАЛЫҚ ҚАТАРЛАРЫ, % ЭЛЕМЕНТТЕР ОРТАСЫ АРҚЫЛЫ ӨТЕТІН
21	% egor+1 ТҮЙІНДЕРДЕН ТҰРАДЫ
22
23	% ЖЫЛЖУЛАРДАҒЫ ШЕКТІК ЖАҒДАЙЛАРДЫ БЕРУ ҮШІН БЕРІЛГЕН МӘНДЕРІ
24	% БАР ТҮЙІНДЕР НӨМІРЛЕРІНІҢ АУҚЫМДАРЫН ҚАЛЫПТАСТЫРУ
25
26	Mz=zeros(400);
27	np=ever*(3*egor+2)+2*egor+1;
28	for i=1:(2*egor+1)
29	Mz(i)=i;
30	end
31	k=2*egor+2;
32	for i=1:ever
33	Mz(k)=(2*egor+1)+1+(i-1)*(3*egor+2);
34	Mz(k+1)=Mz(k)+egor;
35	Mz(k+2)=Mz(k+1)+1;
36	Mz(k+3)=Mz(k+2)+2*egor;
37	k=k+4;
38	end
39	for i=1:(2*egor+1)
40	Mz(k)=Mz(i)+np;
41	k=k+1;
42	end
43
44	% МЕТРЛАРДАҒЫ КООРДИНАТ ӨСТЕРІ БОЙЫНША АУҚЫМДЫ АДЫМДАР
45	format long
46	dx=zeros(2*egor);
47	dy=zeros(2*ever);
48
49	dx=[
50	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
51	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
52	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
53	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
54	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
55	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
56	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
57	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
58	0.250e-01 0.250e-01 0.250e-01 0.250e-01];
59
60	dy=[0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 ...
61	0.1000 0.1000 0.1000 0.1000 0.1000 0.0750 0.0750 ...
62	0.0625 0.0625 0.0500 0.0500 0.0500 0.0500 0.0250 ...
63	0.0250 0.0250];
64
65	yy=[0.000 0.200 0.400 0.600 0.800 1.000 1.200 1.400 1.600 1.800 2.000 2.200 ...
66	2.400 2.550 2.700 2.825 2.950 3.050 3.150 3.250 3.350 3.400 3.450 3.500];
67
68	% КООРДИНАТТАРДЫ ҚАЛЫПТАСТЫРУ БЛОГЫ
69	% х өсіне оң бағыт – солдан оңға
70	% у өсіне оң бағыт – төменнен жоғарыға
71	k=1;
72	NN=zeros(ever);
73	for n=1:ever
74	for m=1:egor
75	n1=2*m-1+(3*egor+2)*(n-1); n2=n1+1; n3=n1+2;
76	n4=n3+2*egor-m+1; n5=n3+(3*egor+2);
77	n6=n5-1; n7=n5-2; n8=n4-1;
78
79	x(n1)=(m-1)*2*dx(m); x(n2)=x(n1)+dx(m); x(n3)=x(n1)+2*dx(m);
80	x(n4)=x(n3); x(n5)=x(n3); x(n6)=x(n2); x(n7)=x(n1); x(n8)=x(n1);
81
82	y(n1)=yy(n); y(n2)=y(n1); y(n3)=y(n1); y(n4)=y(n3)+dy(n);
83	y(n5)=y(n3)+2*dy(n); y(n6)=y(n5); y(n7)=y(n5); y(n8)=y(n4);
84	if(m==36)
85	NN(2*k-1)=n3;
86	NN(2*k)=n4;
87	k=k+1;
88	end
89	end
90	end
91	NN(1:2*ever+1);
92
93	% БАСТАПҚЫ МӘЛІМЕТТЕР ЖӘНЕ МАТРИЦАЛАРДЫҢ НӨЛДЕНУІ
94	np2=2*np;
95	K=zeros(np2,np2);
96	U=zeros(np2+1);
97	U0=zeros(np2+1);
98	F=zeros(np2);
99
100	% ҚҮРЫЛЫМДЫҚ ҚАБАТ МАТЕРИАЛДАРЫНЫҢ АТАУЫ (ЖОҒАРЫДАН ТӨМЕНГЕ)
101	% 1. ҰТ ТЫҒЫЗ АБ
102	% 2. ІТ КЕУЕК АБ
103	% 3. 8% ШАҒЫЛ ТАСТЫ-ҚОСПА
104	% 4. ШАҒЫЛ ТАСТЫ ҚОСПА
105	% 5. ҚИЫРШЫҚ ТАСТЫ-ҚҰМДЫ ҚОСПА
106	% 6. ЖЕҢІЛ БАЛШЫҚ ТОПЫРАҚ
107
108	% ҚАБАТТАР МАТЕРИАЛДАРЫНЫҢ ФИЗИКА-МЕХАНИКАЛЫҚ % СИПАТТАМАЛАРЫ (ЖОҒАРЫДАН ТӨМЕНГЕ)
109	E=[6933e+05 4142e+05 1000e+05 250e+05 180e+05 58e+05]; % Серп. модулі кг/м^2
110	nu=[0.18 0.27 0.30 0.30 0.30 0.35]; % ПУАССОН коэффициенті
111
112	% ЖҮКТЕМЕЛЕРДЕГІ ШЕКТІК ЖАҒДАЙЛАР
113
114	% Дара дөңгелектен түсетін жүктеме
115
116	% ИНТЕРНЕТТЕН АЛЫНҒАН МӘЛІМЕТТЕР
117	% 1 Па=0.102 кгс/м^2
118	% 0.6 МПа=6.12 кгс/см^2=60000.0 кгс/м^2
119
120	% ЭЛЕМЕНТТІҢ БҰРЫШТЫҚ ТҮЙІНДЕРІНЕ ҮЛЕСТІРІЛГЕН ДАРА % ДӨҢГЕЛЕКТЕН ТҮСЕТІН ЖҮКТЕМЕ
121	% ЕСЕПТІК СҰЛБАНЫҢ ОРТАСЫНА ҚОЙЫЛАДЫ. % Аралық түйіндер
122	% ескерілмейді
123
124	% ДӨҢГЕЛЕК ЖҮКТЕМЕСІ НОРМАТИВТІК ҚАРҚЫНДЫЛЫҚ % ЕСЕБІНЕН ЕСЕПТЕЛЕДІ
125	% ЖҮКТЕМЕНІҢ q=0.6 MPa
126	% ЕГЕР 1 Па=0.102 кгс/м^2 (Интернеттен) ОНДА
127	% 0.6 МПа=6.12 кгс/см^2 = 60000.0 кгс/м^2.
128	% Осылайша, егер жамылғының 1 метр бойынша 60000 кг салмақ тура келсе, онда
129	% 0.05 м бір арақашықтыққа – 3000 кг салмақ .
130	% Егер дөңгелектің енін 40 см деп, ал көлденең бағыттағы
131	% бір соңғы элементтің енін – 5 см деп қабылдасақ, онда көлік жүктемесі
132	% 8 соңғы элемент шегінде 9 түйін бойынша үлеседі
133	% Бұрыштық түйіндер арасындғы арақашықтық - 0.05 м (5 cм – соңғы элементтің ені).
134	% Сонымен қатар келесіні ескеру қажет
135	% екі шеткі түйіндерге (1-ші және 9-түйіндер) жартылай жүктеме түседі.
136	% Мұндай жағдайда түйінді нүктелердегі жүктемелер
137	% келесіндей үлестіріледі:
138	% P1=P9=1500 кг, P2=P3= … =P8=3000 кг
139
140	% БІР ТҮЙІНГЕ КЕЛЕТІН ЖҮКТЕМЕНІ ЕСЕПТЕУ
141	% ТҮЙІНДЕР АРАСЫНДАҒЫ АРАҚАШЫҚТЫҚ 0.025 М (аралық түйіндер ескеріледі)
142	% 1 М-ГЕ 60000 КГ САЛМАҚ ТУРА КЕЛЕДІ, АЛ 0.025 М БІР АРАҚАШЫҚТЫҚҚА
143	% 1500 КГ САЛМАҚ ТУРА КЕЛЕДІ
144
145	% НӘТИЖЕЛЕР КГ-ДА ЖӘНЕ МЕТРЛЕРДЕ АЛЫНАДЫ
146
147	% ЭПЮРАНЫҢ ОРТАЛЫҚ БӨЛІГІНЕ ТҮСІРІЛГЕН
148	% ТІК КҮШТЕРДІҢ ШАМАСЫ
149
150	P=0.3e+04; % кг
151
152	% КҮШ ТҮСІРІЛГЕН БҰРЫШТЫҚ ТҮЙІНДЕР НӨМІРЛЕРІ
153
154	nF(1:9)=[4239 4241 4243 4245 4247 4249 4251 4253 4255];
155	for i=1:9
156	j=nF(i);
157	xP(i)=x(j);
158	end
159
160	% БҰРЫШТЫҚ ТҮЙІНДЕРГЕ КҮШ ТҮСІРУ
161	F(nF(1)+np)=-P/2;
162	F(nF(9)+np)=-P/2;
163	for i=1:7
164	F(nF(1+i)+np)=-P;
165	end
166
167	% ЭЛЕМЕНТ ҚАТТЫЛЫҒЫНЫҢ [ke] МАТРИЦАСЫН ЖӘНЕ
168	% ЖҮЙЕ ҚАТТЫЛЫҒЫНЫҢ [K] ҒАЛАМДЫҚ МАТРИЦАСЫН ҚАЛЫПТАСТЫРУ
169
170	mpe=[-0.577350 0.577350]; % ИНТЕГРАЛДАУ НҮКТЕЛЕРЕНІҢ КООРДИНАТТАРЫ
171	mH=[1 1]; % САЛМАҚТЫҚ КОЭФФИЦИЕНТТЕР
172	K=zeros(np2,np2);
173	for n=1:ever
174	for m=1:egor
175	n1=2*m-1+(3*egor+2)*(n-1); n2=n1+1; n3=n1+2;
176	n4=n3+2*egor-m+1; n5=n3+(3*egor+2);
177	n6=n5-1; n7=n5-2; n8=n4-1;
178
179	X1=x(n1); Y1=y(n1);
180	X2=x(n3); Y2=y(n3);
181	X3=x(n5); Y3=y(n5);
182	X4=x(n7); Y4=y(n7);
183
184	% ЖАЗЫҚ ДЕФОРМАЦИЯ ЖАҒДАЙЫ ҮШІН СЕРПІМДІЛІК МАТРИЦАСЫН % ҚАЛЫПТАСТЫРУ
185	if((n>=1)&(n<=7))
186	e1=E(6); nu1=nu(6); % ЖЕҢІЛ БАЛШЫҚ ТОПЫРАҚ
187	end
188	if((n>=8)&(n<=12))
189	e1=E(5); nu1=nu(5); % ҚИЫРШЫҚ ТАСТЫ-ҚҰМДЫ ҚОСПА
190	end
191	if((n>=13)&(n<=16))
192	e1=E(4); nu1=nu(4); % ШАҒЫЛ ТАСТЫ ҚОСПА
193	end
194	if((n>=17)&(n<=20))
195	e1=E(3); nu1=nu(3); % 8% ШАҒЫЛ ТАТЫ-ҚҰМДЫ ҚОСПА
196	end
197	if((n>=21)&(n<=22))
198	e1=E(2); nu1=nu(2); % ІТ КЕУЕК АБ
199	end
200	if(n==23)
201	e1=E(1); nu1=nu(1); % ҰТ ТЫҒЫЗ АБ
202	end
203
204	D=zeros(3,3);
205	D(3,3)=e1/(2*(1+nu1)); D(2,2)=2*D(3,3)*(1-nu1)/(1-2*nu1);
206	D(1,1)=D(2,2); D(1,2)=2*D(3,3)*nu1/(1-2*nu1);
207	D(2,1)=D(1,2);
208
209	dJ=abs((X1-X2)*(Y2-Y4))/4; % якобиан
210
211	% ЭЛЕМЕНТ ҚАТТЫЛЫҒЫНЫҢ[Ke] МАТРИЦАСЫНДАҒЫ ҚОС ИНТЕГРАЛДЫ % САНДЫҚ ИНТЕГРАЛДАУ ҮШІН
212	% ГАУСС КВАДРАТУРАСЫН ПАЙДАЛАНУ
213
214	% BtDBnds1(psi,eta,D,X1,X2,Y1,Y3) – ішкі есептеу бағдарламасы
215	% элемент қаттылығының матрицасының
216	Ke=zeros(16,16);
217	psi=mpe(1); eta=mpe(1);
218	k11=BtDBnds1(psi,eta,D,X1,X2,Y1,Y3);
219	psi=mpe(1); eta=mpe(2);
220	k12=BtDBnds1(psi,eta,D,X1,X2,Y1,Y3);
221	psi=mpe(2); eta=mpe(1);
222	k21=BtDBnds1(psi,eta,D,X1,X2,Y1,Y3);
223	psi=mpe(2); eta=mpe(2);
224	k22=BtDBnds1(psi,eta,D,X1,X2,Y1,Y3);
225	Ke=dJ*(mH(1)*(mH(1)*k11+mH(2)*k12)+mH(2)*(mH(1)*k21+mH(2)*k22));
226
227	ne8=[n1 n2 n3 n4 n5 n6 n7 n8];
228	for i=1:8
229	for j=1:8
230	K(ne8(i),ne8(j))=K(ne8(i),ne8(j))+Ke(i,j);
231	K(ne8(i),ne8(j)+np)=K(ne8(i),ne8(j)+np)+Ke(i,j+8);
232	K(ne8(i)+np,ne8(j))=K(ne8(i)+np,ne8(j))+Ke(i+8,j);
233	K(ne8(i)+np,ne8(j)+np)=K(ne8(i)+np,ne8(j)+np)+Ke(i+8,j+8);
234	end
235	end
236	end
237	end
238
239	disp('К матрицасы қалыптастырылды ')
240
241	% САТЖ ТҮРЛЕНДІРУ
242
243	np2=2*np;
244	i=1;
245	for n=1:np2
246	if(n==Mz(i))
247	F(n)=K(n,n)*U(n);
248	for m=1:np2
249	if(m~=n)
250	K(n,m)=0.0;
251	end
252	end
253	i=i+1;
254	end
255	end
256
257	i=1;
258	for n=1:np2
259	if(n==Mz(i))
260	for m=1:np2
261	if(m~=n)
262	F(m)=F(m)-K(m,n)*U(n);
263	K(m,n)=0.0;
264	end
265	end
266	i=i+1;
267	end
268	end
269
270	% САТЖ ШЕШІМІ
271	U=K\F;
272	disp(' [K]{U}={F} жүйесі шешілді ')
273
274	% Сипатты тік қималардағы түйіндер нөмірлері
275	nom12=[
276	49 146 231 328 413 510 595 692 ...
277	777 874 959 1056 1141 1238 1323 1420 ...
278	1505 1602 1687 1784 1869 1966 2051 2148 ...
279	2233 2330 2415 2512 2597 2694 2779 2876 ...
280	2961 3058 3143 3240 3325 3422 3507 3604 ...
281	3689 3786 3871 3968 4053 4150 4235];
282	nom15=[
283	61 152 243 334 425 516 607 698 ...
284	789 880 971 1062 1153 1244 1335 1426 ...
285	1517 1608 1699 1790 1881 1972 2063 2154 ...
286	2245 2336 2427 2518 2609 2700 2791 2882 ...
287	2973 3064 3155 3246 3337 3428 3519 3610 ...
288	3701 3792 3883 3974 4065 4156 4247];
289	nom18=[
290	73 158 255 340 437 522 619 704 ...
291	801 886 983 1068 1165 1250 1347 1432 ...
292	1529 1614 1711 1796 1893 1978 2075 2160 ...
293	2257 2342 2439 2524 2621 2706 2803 2888 ...
294	2985 3070 3167 3252 3349 3434 3531 3616 ...
295	3713 3798 3895 3980 4077 4162 4259];
296
297	% ОРТАША ҚИМАДАҒЫ ТҮЙІНДЕРДІҢ КӨЛДЕНЕҢ ЖЫЛЖУЛАР КЕСТЕЛЕРІ
298	Y=zeros(2*ever+1);
299	X150=zeros(2*ever+1);
300	X120=zeros(2*ever+1);
301	X180=zeros(2*ever+1);
302	for i=1:2*ever+1
303	j120=nom12(i);
304	j150=nom15(i);
305	j180=nom18(i);
306	Y(i)=y(j150);
307	X120(i)=U(j120)-0.0001; % 30 СМ ТІК СИММЕТРИЯ ӨСІНЕН СОЛ ЖАҚТА
308	X150(i)=U(j150)*10; % ТІК СИММЕТРИЯ ӨСІНДЕ
309	X180(i)=U(j180)+0.0001; % 30 СМ ТІК СИММЕТРИЯ ӨСІНЕН % ОҢ ЖАҚТА
310	end
311	hPlot=plot(X150,Y,'b-',X120,Y,'r-',X180,Y,'k-.');grid on;xlabel('Ux, m');ylabel('y, m')
312	set(hPlot,'LineWidth',2);
313
314	figure
315
316	% ТҮЙІНДЕРДІҢ КӨЛДЕНЕҢ ҚАТАРЛАРЫНЫҢ АУҚЫМДАРЫН ҚАЛЫПТАСТЫРУ % Nv(I,m)
317
318	% 2*egor+1 ТЕҢ ӘРБІР ТОЛЫҚ КӨЛДЕНЕҢ ҚАТАРДАҒЫ ТҮЙІНДЕРДІҢ % ЖАЛПЫ САНЫ
319	% ОСЫНДАЙ КӨЛДЕНЕҢ ҚАТАРЛАРДЫҢ САНЫ ever+1 ТЕҢ
320	% ЭЛЕМЕНТТЕР ОРТАСЫНАН ӨТЕТІН ҚИМАЛАР % ҚАРАСТЫРЫЛМАЙДЫ.
321	% ОЛАР ӘРБІР КӨЛДЕНЕҢ ДЕҢГЕЙДЕ (ТЕРҢДІКТЕ) ІЗДЕЛЕТІН % ШАМАЛАРДЫҢ КЕСТЕЛЕРІН ҚҰРУ ҮШІН
322	% ПАЙДАЛАНЫЛАДЫ
323
324	Nv=zeros(ever+1,2*egor+1);
325	for n=1:(ever+1)
326	for m=1:(2*egor+1)
327	I=2*n-1;
328	Nv(I,m)=m+(3*egor+2)*(n-1);
329	end
330	end
331	for n=1:ever
332	for m=1:(egor+1)
333	I=2*n;
334	Nv(I,2*m-1)=m+(2*egor+1)*n+(egor+1)*(n-1);
335	end
336	end
337
338	format short
339	Set=zeros(2*ever+1,4);
340	for n=1:2*ever+1
341	i=Nv(n,1);
342	Set(n,1)=n;
343	Set(n,2)=i;
344	Set(n,3)=y(i);
345	Set(n,4)=3.50-y(i);
346	end
347	for m=1:47
348	% fprintf('%6u %6u %10.4f %10.4f\n',Set(m,1),Set(m,2),Set(m,3),Set(m,4));
349	end
350
351	% һ ТЕРЕҢДІКТЕРІНДЕГІ КӨЛДЕНЕҢ ҚИМАЛАРДЫҢ ТІК ЖЫЛЖУЛАРЫНЫҢ
352	% КЕСТЕЛЕРІ
353	X=zeros(2*egor+1);
354	Y=zeros(2*egor+1);
355	for i=1:(2*egor+1)
356	j=Nv(45,i); % ТЕРЕҢДІК h=5 СМ
357	X(i)=x(j);
358	Y(i)=U(j+np);
359	end
360	hPlot=plot(X,Y,'-');grid on;xlabel('X, m');ylabel('Uy, m')
361	set(hPlot,'LineWidth',2);
362	hold on
363
364	X=zeros(2*egor+1);
365	Y=zeros(2*egor+1);
366	for i=1:(2*egor+1)
367	j=Nv(41,i); % ТЕРЕҢДІК h=15 СМ
368	X(i)=x(j);
369	Y(i)=U(j+np);
370	end
371	hPlot=plot(X,Y,'-');grid on;xlabel('X, m');ylabel('Uy, m')
372	set(hPlot,'LineWidth',2);
373	hold on
374
375	X=zeros(2*egor+1);
376	Y=zeros(2*egor+1);
377	for i=1:(2*egor+1)
378	j=Nv(39,i); % ТЕРЕҢДІК h=25 СМ
379	X(i)=x(j);
380	Y(i)=U(j+np);
381	Ux325(i)=U(j);
382	Uy325(i)=Y(i);
383	end
384	hPlot=plot(X,Y,'--');grid on;xlabel('X, m');ylabel('Uy, m')
385	set(hPlot,'LineWidth',2);
386	hold on
387
388	for i=1:(2*egor+1)
389	j=Nv(35,i); % ТЕРЕҢДІК h=45 СМ
390	X(i)=x(j);
391	Y(i)=U(j+np);
392	X345(i)=X(i);
393	Ux305(i)=U(j);
394	Uy305(i)=Y(i);
395	end
396	hPlot=plot(X,Y,'--');grid on;xlabel('X, m');ylabel('Uy, m')
397	set(hPlot,'LineWidth',2);
398	hold on
399
400	for i=1:(2*egor+1)
401	j=Nv(29,i); % ТЕРЕҢДІК h=80 СМ
402	X(i)=x(j);
403	Y(i)=U(j+np);
404	Ux270(i)=U(j);
405	Uy270(i)=Y(i);
406	end
407	hPlot=plot(X,Y,'-.');grid on;xlabel('X, m');ylabel('Uy, m')
408	set(hPlot,'LineWidth',2);
409
410	figure
411
412	% ЭЛЕМЕНТТЕРДІҢ ОРТАЛЫҚ НҮКТЕЛЕРІНДЕГІ ДЕФОРМАЦИЯЛАР
413	% КОМПОНЕНТТЕРІН ЕСЕПТЕУ (eps=0.0, eta=0.0)
414
415	epsX0=zeros(ever*egor);
416	epsY0=zeros(ever*egor);
417	gamXY0=zeros(ever*egor);
418
419	for n=1:ever
420	for m=1:egor
421	n1=2*m-1+(3*egor+2)*(n-1); n2=n1+1; n3=n1+2;
422	n4=n3+2*egor-m+1; n5=n3+(3*egor+2);
423	n6=n5-1; n7=n5-2; n8=n4-1;
424
425	dJ=abs((x(n1)-x(n3))*(y(n3)-y(n7)))/4; % якобиан
426	i=m+egor*(n-1);
427	epsX0(i)=(U(n4)-U(n8))/(x(n4)-x(n8));
428	epsY0(i)=(-U(n2+np)+U(n6+np))/(y(n6)-y(n2));
429	gamXY0(i)=-(-U(n2)+U(n6))/(x(n4)-x(n8))+(U(n4+np)-U(n8+np))/(y(n6)-y(n2));
430	end
431	end
432
433	% ЖОЛМА-ЖОЛ ЭЛЕМЕНТТЕР НӨМІРІ
434	NEv=zeros(ever,egor);
435	for n=1:ever
436	for m=1:egor
437	i=m+egor*(n-1);
438	NEv(n,m)=i;
439	end
440	end
441
442	% ДЕФОРМАЦИЯЛАР МЕН КЕРНЕУЛЕР ЕСЕПТЕЛІНЕТІН ОРТАЛЫҚ % КЕРНЕУДІҢ НӨМІРІ
443	% NEv(n,m) АУЫҚЫМЫНДА ОРНАЛАСҚАН ЭЛЕМЕНТ НӨМІРІНЕ СӘЙКЕС КЕЛЕДІ
444
445	X=zeros(egor); Y25=zeros(egor); Y200=zeros(egor); Y400=zeros(egor);
446	for m=1:egor
447	n=Nv(1,2*m);
448	X(m)=x(n);
449	i25=NEv(23,m);
450	i200=NEv(20,m);
451	i400=NEv(18,m);
452	Y25(m)=epsX0(i25); % 2,5 cм
453	Y200(m)=epsX0(i200); % 20,0 cм
454	Y400(m)=epsX0(i400); % 40,0 см
455	end
456	hPlot=plot(X,Y25,'-',X,Y200,'--',X,Y400,'-.');...
457	grid on;xlabel('X, m');ylabel('epsX')
458	set(hPlot,'LineWidth',2);
459	figure
460
461	X=zeros(egor); Y25=zeros(egor); Y200=zeros(egor); Y400=zeros(egor);
462	for m=1:egor
463	n=Nv(1,2*m);
464	X(m)=x(n);
465	i25=NEv(23,m);
466	i200=NEv(20,m);
467	i400=NEv(18,m);
468	Y25(m)=epsY0(i25); % 2,5 cм
469	Y200(m)=epsY0(i200); % 20,0 cм
470	Y400(m)=epsY0(i400); % 40,0 см
471	end
472	hPlot=plot(X,Y25,'-',X,Y200,'--',X,Y400,'-.');...
473	grid on;xlabel('X, m');ylabel('epsY')
474	set(hPlot,'LineWidth',2);
475	figure
476
477	% ЭЛЕМЕНТТЕРДІҢ ОРТА НҮКТЕЛЕРІНДЕГІ Sigma КЕРНЕУЛЕР
478	% КОМПОНЕНТТЕРІН ЕСЕПТЕУ (eps=0.0, eta=0.0)
479	SigX=zeros(ever); SigY=zeros(ever);
480	k=0;
481	for n=1:ever
482	for m=1:egor
483	n1=2*m-1+(3*egor+2)*(n-1); n2=n1+1; n3=n1+2;
484	n4=n3+2*egor-m+1; n5=n3+(3*egor+2);
485	n6=n5-1; n7=n5-2; n8=n4-1;
486
487	% ЖАЗЫҚ ДЕФОРМАЦИЯ ЖАҒДАЙЫ ҮШІН СЕРПІМДІЛІК МАТРИЦАСЫН % ҚАЛЫПТАСТЫРУ
488	if((n>=1)&(n<=7))
489	e1=E(6); nu1=nu(6); % ЖЕҢІЛ БАЛШЫҚ ТОПЫРАҚ
490	end
491	if((n>=8)&(n<=12))
492	e1=E(5); nu1=nu(5); % ҚИЫРШЫҚ ТАСТЫ-ҚҰМДЫ ҚОСПА
493	end
494	if((n>=13)&(n<=16))
495	e1=E(4); nu1=nu(4); % ШАҒЫЛ ТАСТЫ ҚОСПА
496	end
497	if((n>=17)&(n<=20))
498	e1=E(3); nu1=nu(3); % 8% ШАҒЫЛ ТАСТЫ-ҚҰМДЫ ҚОСПА
499	end
500	if((n>=21)&(n<=22))
501	e1=E(2); nu1=nu(2); % ІТ КЕУЕК АБ
502	end
503	if(n==23)
504	e1=E(1); nu1=nu(1); % ҰТ ТЫҒЫЗ АБ
505	end
506
507	D=zeros(3,3);
508	D(3,3)=e1/(2*(1+nu1)); D(2,2)=2*D(3,3)*(1-nu1)/(1-2*nu1);
509	D(1,1)=D(2,2); D(1,2)=2*D(3,3)*nu1/(1-2*nu1);
510	D(2,1)=D(1,2);
511
512	X1=x(n1); X2=x(n3); Y2=y(n3); Y4=y(n7);
513	dJ=abs((X1-X2)*(Y2-Y4))/4; % Якобиан
514	G=e1/(2*(1+nu1));
515
516	i=m+egor*(n-1);
517	% В МПа
518	SigmaX0(i)=(D(1,1)*epsX0(i)+D(1,2)*epsY0(i))/100000;
519	SigmaY0(i)=(D(2,1)*epsX0(i)+D(2,2)*epsY0(i))/100000;
520	TauXY0(i)=D(3,3)*gamXY0(i)/100000;
521	if(m==30)
522	k=k+1;
523	SigX(n)=SigmaX0(i);
524	SigY(n)=SigmaY0(i);
525	end
526	end
527	end
528
529	% КЕРНЕУЛЕР КОМПОНЕНТТЕРІНІҢ КЕСТЕСІН ҚҰРУ SigmaX0, SigmaY0 И
530	% TauXY0 ЭЛЕМЕНТТЕРДІҢ ОРТА ТҮЙІНДЕРІНДЕ
531
532	X=zeros(egor); Y23=zeros(egor); Y22=zeros(egor); Y20=zeros(egor);
533	Y18=zeros(egor); Y14=zeros(egor); Y12=zeros(egor);
534	for m=1:egor
535	n=Nv(1,2*m);
536	X(m)=x(n);
537	i23=NEv(23,m); % h=2,5 см
538	i22=NEv(22,m); % h=7,5 см
539	i20=NEv(20,m); % h=20,0 см
540	i18=NEv(18,m); % h=40,0 см
541	i14=NEv(14,m); % h=87,5 см
542	i12=NEv(12,m); % h=120,0 см
543	Y23(m)=SigmaX0(i23); % 2,5 cм
544	Y22(m)=SigmaX0(i22); % 7,5 cм
545	Y20(m)=SigmaX0(i20); % 20,0 cм
546	Y18(m)=SigmaX0(i18); % 40,0 cм
547	Y14(m)=SigmaX0(i14); % 87,5 cм
548	Y12(m)=SigmaX0(i12); % 120,0 cм
549	end
550	hPlot=plot(X,Y22,'-',X,Y20,'--',X,Y18,'-.',X,Y14,'-r');...
551	grid on;xlabel('X, m');ylabel('SigmaX, МПа')
552	set(hPlot,'LineWidth',2);
553	figure
554
555	X=zeros(egor); Y23=zeros(egor); Y22=zeros(egor); Y20=zeros(egor);
556	Y18=zeros(egor); Y14=zeros(egor); Y12=zeros(egor);
557	for m=1:egor
558	n=Nv(1,2*m);
559	X(m)=x(n);
560	i23=NEv(23,m); % h=2,5 см
561	i22=NEv(22,m); % h=7,5 см
562	i20=NEv(20,m); % h=20,0 см
563	i18=NEv(18,m); % h=40,0 см
564	i14=NEv(14,m); % h=87,5 см
565	i12=NEv(12,m); % h=120,0 см
566	Y23(m)=SigmaY0(i23); % 2,5 cм
567	Y22(m)=SigmaY0(i22); % 7,5 cм
568	Y20(m)=SigmaY0(i20); % 20,0 cм
569	Y18(m)=SigmaY0(i18); % 40,0 cм
570	Y14(m)=SigmaY0(i14); % 87,5 cм
571	Y12(m)=SigmaY0(i12); % 120,0 cм
572	end
573	hPlot=plot(X,Y22,'-',X,Y20,'--',X,Y18,'-.',X,Y14,'-r');...
574	grid on;xlabel('X, m');ylabel('SigmaY, МПа')
575	set(hPlot,'LineWidth',2);
576	figure
577
578	X=zeros(egor); Y23=zeros(egor); Y22=zeros(egor); Y20=zeros(egor);
579	Y18=zeros(egor); Y14=zeros(egor); Y12=zeros(egor);
580	for m=1:egor
581	n=Nv(1,2*m);
582	X(m)=x(n);
583	i23=NEv(23,m); % h=2,5 см
584	i22=NEv(22,m); % h=7,5 см
585	i21=NEv(21,m); % h=12,5 см
586	i20=NEv(20,m); % h=20,0 см
587	i18=NEv(18,m); % h=40,0 см
588	i16=NEv(16,m); % h=61,25 см
589	i14=NEv(14,m); % h=87,5 см
590	i12=NEv(12,m); % h=120,0 см
591	Y23(m)=TauXY0(i23); % 2,5 cм
592	Y22(m)=TauXY0(i22); % 7,5 cм
593	Y21(m)=TauXY0(i21); % 12,5 cм
594	Y20(m)=TauXY0(i20); % 20,0 cм
595	Y18(m)=TauXY0(i18); % 40,0 cм
596	Y16(m)=TauXY0(i16); % 61,25 cм
597	Y14(m)=TauXY0(i14); % 87,5 cм
598	Y12(m)=TauXY0(i12); % 120,0 cм
599	end
600	hPlot=plot(X,Y22,'-',X,Y20,'--',X,Y14,'-.');...
601	grid on;xlabel('X, m');ylabel('TauXY, МПа')
602	set(hPlot,'LineWidth',2);
603
604	toc

Б қосымшасы (міндетті) BtDBnds кіші бағдарламасының бастапқы коды

1	function[BDB]=BtDBnds(psi,eta,D,X1,X2,Y1,Y3);
2	A=(-2/(X1-X2));
3	B1=A*(1-eta)*(2*psi+eta)/4; B2= A*(-psi*(1-eta));
4	B3= A*(1-eta)*(2*psi-eta)/4; B4= A*(1-eta^2)/2;
5	B5= A*(1+eta)*(2*psi+eta)/4; B6= A*(-psi*(1+eta));
6	B7= A*(1+eta)*(2*psi-eta)/4; B8= A*(-(1-eta^2)/2);
7	B=(-2/(Y1-Y3));
8	C1=B*(1-psi)*(psi+2*eta)/4; C2=-B*(1-psi^2)/2;
9	C3=B*(1+psi)*(-psi+2*eta)/4; C4=-B*eta*(1+psi);
10	C5=B*(1+psi)*(psi+2*eta)/4; C6=B*(1-psi^2)/2;
11	C7=B*(-(1-psi)*(psi-2*eta)/4); C8=-B*eta*(1-psi);
12
13
14	B=[B1 B2 B3 B4 B5 B6 B7 B8 0 0 0 0 0 0 0 0;
15	0 0 0 0 0 0 0 0 C1 C2 C3 C4 C5 C6 C7 C8;
16	C1 C2 C3 C4 C5 C6 C7 C8 B1 B2 B3 B4 B5 B6 B7 B8];
17	Bt=B';
18	BtD=Bt*D;
19	BDB=BtD*B;

Библиография

[1] Мартынов Н.Н., Иванов А.П. MATLAB 5.Х. Есептеу, визуалдау, бағдарламалау. – М.:КУДИЦ-ОБРАЗ, 2000.-336 б.

[2] Коткин Г.Л., Черкасский В.С. MATLAB қолданып физикалық үдерістерді компьютерлік моделдеу: Оқу құралы/Новосиб. Новосибирск ун-т., 2001. 173 б.

[3] Потемкин В.Г. MATLAB 5.Х инженерлік және ғылыми есептеулер жүйесі 2 томда. –М.:ДИАЛОГ-МИФИ, 1999. 670 б.

[4] Тынкевич М.А. Сандық әдістер.-Кемерово.:КузГТУ.1977.-122 б.

[5] Сегерлинд Л. Шекті элементтер әдісін қолдану. М: Мир, 1979.–392 б.

[6] ҚР СТ 1293-2004 Қатты емес типті жол төсемелерінің серпімділік модулін анықтау әдістері және олардың жіктелуі.

Р РК 218-149-2018 МЕТОДИЧЕСКИЕ РЕКОМЕНДАЦИИ ПО ПРАКТИЧЕСКОМУ ИСПОЛЬЗОВАНИЮ ПРОГРАММЫ ПО ОПРЕДЕЛЕНИЮ НАПРЯЖЕННО-ДЕФОРМИРОВАННОГО СОСТОЯНИЯ МНОГОСЛОЙНОЙ КОНСТРУКЦИИ ДОРОЖНОЙ ОДЕЖДЫ

Приказ Председателя Комитета автомобильных дорог Министерства по инвестициям и развитию Республики Казахстан от 21 декабря 2018 года № 122

Предисловие

1	РАЗРАБОТАНЫ И ВНЕСЕНЫ	Акционерным обществом "Казахстанский дорожный научно-исследовательский институт" (АО "КаздорНИИ")
2	УТВЕРЖДЕНЫ И ВВЕДЕНЫ В ДЕЙСТВИЕ	Приказом Председателя Комитета автомобильных дорог Министерства по инвестициям и развитию Республики Казахстан № 122 от 21 декабря 2018 года
3	СОГЛАСОВАНЫ	Акционерным обществом "НК "ҚазАвтоЖол" № 03/14-2-2591-И от 12 ноября 2018 года
4	СРОК ПЕРВОЙ ПРОВЕРКИ	2023 год
	ПЕРИОДИЧНОСТЬ ПРОВЕРКИ	5 лет
5	ВВЕДЕНЫ ВПЕРВЫЕ

Документ доступен к просмотру в информационно-правовой системе нормативно-правовых актов Республики Казахстан "Әдiлет"

Настоящие рекомендации не могут быть полностью или частично воспроизведены, тиражированы и распространены без разрешения Комитета автомобильных дорог Министерства по инвестициям и развитию Республики Казахстан

Содержание

	Введение
1	Область применения
2	Нормативные ссылки
3	Термины и определения
4	Методика использования программы BASIC_NDS_MKE_8_uzlov
4.1	Общая постановка задачи
4.2	Формирование массивов координат
4.3	Граничные условия
4.4	Определяющие соотношения уравнения равновесия
4.5	Вычисление матрицы градиентов четырехугольного элемента с восемью узлами
4.6	Алгоритм определения деформации и напряжений
4.7	Краткая инструкция по использованию программы BASIC_NDS_MKE_8_uzlov
4.8	Примеры расчетов по программе BASIC_NDS_MKE_8_uzlov
Приложение А (обязательное) Исходный код программы BASIC_NDS_MKE_8_uzlov Приложение Б (обязательное) Исходный код подпрограммы BtDBnds Библиография

Введение

Настоящие методические рекомендации разработаны согласно плану работ, принятому в рамках выполнения темы Б01.02 "Разработка теоретических решений и программ для определения напряженно-деформированного состояния многослойной конструкции дорожной одежды".

Документ содержит постановку типовой задачи по определению напряженно-деформированного состояния многослойной дорожной конструкции, статически нагруженной весом автомобиля, и основные формулы из алгоритма расчета методом конечных элементов. В методе применяется высокоточный квадратичный прямоугольный конечный элемент с восемью узлами. В приложении А приведен текст (листинг) расчетной программы BASIC_NDS_MKE_8_uzlovна языке MATLAB [1-4], снабженный необходимыми комментариями. В приложении Б приводится листинг подпрограммы BDB0, предназначенный для вычисления матрицы жесткости элемента. Приводится также краткая инструкция по использованию программы BASIC_NDS_MKE_8_uzlov. В конце документа приведен список использованной литературы.

1 Область применения

1.1 Настоящие рекомендации распространяются на сеть автомобильных дорог общего пользования Республики Казахстан и предназначены для решения вопросов, связанных с проектированием многослойных автомобильных дорог общего пользования.

1.2 Рекомендациями следует руководствоваться при проектировании конструкций дорожных одежд для автомобильных дорог общего пользования, для расчета дорожных одежд на стадиях проектирования и эксплуатации, а также при решении инженерно-экономических задач применительно к автомобильным дорогам.

2 Нормативные ссылки

Для применения настоящих рекомендаций необходимы следующие ссылочные нормативные документы:

СП РК 3.03-103-2014 "Проектирование жестких дорожных одежд"

СП РК 3.03-104-2014 "Проектирование нежестких дорожных одежд"

Примечание - При пользовании настоящими рекомендациями целесообразно проверить действие ссылочных документов по ежегодно издаваемому информационному указателю "Нормативные документы по стандартизации", составленному по состоянию на текущий год и соответствующим ежемесячно издаваемым информационным указателям, опубликованным в текущем году. Если ссылочный документ заменен (изменен), то при пользовании настоящими рекомендациями следует руководствоваться замененным (измененным) стандартом. Если ссылочный документ отменен без замены, то положение, в котором дана ссылка на него, применяется в части, не затрагивающей это ссылку.

3 Термины и определения

В настоящих рекомендациях применяются следующие термины с соответствующими определениями:

3.1 Дорожная одежда: Многослойная конструкция в пределах проезжей части автомобильной дороги, воспринимающая нагрузку от автотранспортного средства и передающая ее на грунт.

3.2 Земляное полотно: Конструктивный элемент, служащий основанием для размещения дорожной одежды, а также технических средств организации дорожного движения и обустройства автомобильной дороги.

3.3 Деформация: Относительная величина, характеризующая изменение линейных размеров тела по отношению к его первоначальным размерам.

3.4 Напряжение: Относительная величина, определяемая нормальной или касательной нагрузкой, приходящейся на единицу площади.

3.5 Дорожная одежда нежесткая: Дорожная одежда со слоями, устроенными из разного вида асфальтобетонов, из материалов и грунтов укрепленных битумом, цементом, известью, комплексными и другими вяжущими, а также из слабосвязных зернистых материалов

3.6 Конструктивный слой: Каждый слой дорожной одежды, состоящий из однородных материалов и отличающийся от соседних слоев видом материалов, его прочностью и составом. Учитывается при расчете прочности дорожной одежды.

3.7 Покрытие дорожное: Одно- или многослойная верхняя часть дорожной одежды, устраиваемая на дорожном основании, непосредственно воспринимающая нагрузки от транспортных средств и предназначенная для обеспечения заданных эксплуатационных требований и защиты дорожного основания от воздействия атмосферных факторов.

3.8 Основание дорожное: Нижний несущий слой дорожной одежды, воспринимающий нагрузки от транспортных средств совместно с покрытием и предназначенный для ее распределения на дополнительные слои или непосредственно на грунт земляного полотна.

3.9 Метод конечных элементов: Численный метод решения дифференциальных уравнений математической физики.

4 Методика использования программы BASIC_NDS_MKE_8_uzlov

4.1 Общая постановка задачи

Для построения математической модели автомобильной дороги общего пользования рассмотрим конструкцию, состоящую из двухслойной асфальтобетонной дорожной одежды и основания дорожной одежды, состоящего из щебеночно-песчаной смеси с 8 % цементной добавкой, щебеночной смеси и гравийно-песчаной смеси. Основание дорожной одежды устраивается на грунтовом основании из легкого суглинка (рисунок 1).

Возможности численного метода – метода конечных элементов, применяемого для решения задачи о равновесии многослойной конструкции, позволяют сравнительно легко назначать требуемые значения геометрических размеров конструктивных слоев и физико-механических характеристик материалов, используемых в конструкции.

Поэтому, приведенные в таблице 1 геометрические размеры конструктивных слоев в рассматриваемой математической модели (рисунок 1) и физико-механические свойства материалов конструктивных элементов (в том числе модули упругости) назначены согласно [6].

1 – покрытие, асфальтобетон мелкозернистый плотный; 2 – асфальтобетон крупнозернистый пористый; 3 – щебеночно-песчаная смесь 8 %; 4 – щебеночная смесь; 5 –гравийно-песчаная смесь; 6 – грунт суглинок легкий

Рисунок 1 – Схематический вид дорожной конструкции

Значения физико-механических свойств материалов конструктивных элементов вводятся в 108-110 строках приложения А.

Таблица 1 – Характеристики конструктивных слоев

Материал конструктивного слоя	Модуль упругости Е, МПа	Коэффициент Пуассона,	Высота слоя, м
1 Асфальтобетон мелкозернистый плотный	3200	0,18	0,05
2 Асфальтобетон крупнозернистый пористый	2000	0,27	0,10
3 Щебеночно-песчаная смесь 8%	800	0,30	0,40
4 Щебеночная смесь	275	0,30	0,45
5 Гравийно-песчаная смесь	180	0,30	1,00
6 Грунт суглинок легкий	50	0,35	-

Для компьютерной реализации алгоритма решения задачи о напряженно-деформированном состоянии дорожной конструкции методом конечных элементов разработана расчетная программа BASIC_NDS_MKE_8_uzlovна алгоритмическом языке системы MATLAB.

Исследуемая область разбивается на квадратичные прямоугольные конечные элементы с восемью узлами (рисунок 2). Местная нумерация узлов с 1 по 8 в таком конечном элементе начинается с левого нижнего узла, и осуществляется в направлении против часовой стрелки.

Рисунок 2 – Прямоугольный элемент с восемью узлами

Для создания конечно-элементной сетки применялись переменные шаги и в направлении координатных осей и (рисунок 3).

В каждом горизонтальном ряде (рисунок 3) расположено по 60 элементов (egor=60), а количество горизонтальных рядов элементов – 23 (ever=23). Следовательно, исследуемая область разбита на 1380 конечных элементов, а число узлов равно np=(3*ever+2)*egor+2*ever+1=4307 (строки 12-15 приложения А).

Рисунок 3 – Схематический вид конечно-элементной разбивки исследуемой области

Глобальная нумерация узлов начинается с крайнего левого вертикального ряда узлов, и осуществляется сверху вниз и перемещается слева направо (рисунок 3).

При этом разбивка получается нерегулярной, т.к. в восьмиузловых элементах центральный узел отсутствует. В данном случае, в одном вертикальном ряде имеются 47 узлов, а в следующем – 24 узла.

Требуется определить компоненты вектора перемещений в узловых точках конечно-элементной сетки под действием веса автомобиля, передаваемого через одно колесо и приложенного в виде вертикальных нагрузок в соответствующих узлах (рисунок 4).

L и H – ширина и высота исследуемой области соответственно; P – нагрузка от колеса транспортного средства; u_i, v_i – компоненты перемещения точек на границах в направлении координатных осей x и y; I – VI конструктивные слои дорожной одежды и грунтового основания

Рисунок 4 – Расчетная схема задачи

4.2 Формирование массивов координат

В программе BASIC_NDS_MKE_8_uzlovзаданы массивы dx и dy, содержащие значения шагов по координатам. Наличие таких массивов позволяет задавать переменные шаги по координатам (строки 44-63 приложения А).

В блоке формирования координат программы (строки 68-91 приложения А) глобальные номера узлов текущего конечного элемента обозначены через n1, n2, …, n8, соответствующие локальным номерам 1,2, …, 8, приведенным на рисунке 2. Конкретные значения глобальных номеров узлов, вычисляемые по алгоритмам вычисления n1, n2, …, n8 в программе BASIC_NDS_MKE_8_uzlovопределяются номером элемента в текущем горизонтальном ряде узлов (m=1:egor) и номером элемента в текущем вертикальном ряде узлов (n=1:ever).

4.3 Граничные условия

Из рисунка 4 видно, что граничные условия в перемещениях по бокам и на нижней границе исследуемой области задаются в виде нулевых значений соответствующих компонентов вектора перемещений (u=0 и u=v=0). Чтобы задавать граничные условия в виде известных перемещений некоторых узлов также нужно создать специальный массив Mz, содержащий номера узлов задаваемых известных перемещений (строки 23-43 приложения А).

Граничные условия в нагрузках задаются в виде сосредоточенных вертикальных усилий, рассчитанных от веса автомобиля, приходящегося на одно колесо. Методика расчета граничных условий в нагрузках приведена в приложении А (строки 112-166). Конечный результат расчета представляет собой вектор правой части уравнения равновесия системы (1).

4.4 Определяющие соотношения уравнения равновесия

Уравнение равновесия системы в случае применения квадратичных элементов, как и в случае линейных треугольных конечных элементов, в матричной форме имеет вид [7]:

(1)

где [K] – матрица жесткости системы размерности (2*np,2*np);

– вектор узловых перемещений размерности 2*np;

– вектор узловых нагрузок размерности 2*np;

np – общее число узлов в конечно-элементной сетке.

Порядок вычисления значений компонентов вектора , соответствующие граничным условиям в нагрузках приведены в строках 150-165приложения А. Остальные компоненты этого вектора, в силу условий равновесия узлов, равны нулю.

В уравнении (1) матрица жесткости системы определяется в виде суммы

где – общее число элементов в конечно-элементной сетке; – матрица жесткости текущего элемента "".

Матрица жесткости элемента вычисляется с помощью объемного интеграла (2):

(2)

где алгоритм вычисления матрицы упругости приведены в строках 204-207 приложения А.

В выражении (2) так называемая матрица градиентов формируется путем дифференцирования функций формы по координатам .

В случае нелинейной зависимости функций формы от координат , что имеет место в случае квадратичных элементов с восемью узлами, производные от функций формы по координатам сами являются функциями координат , и вычисление интеграла представляет собой самостоятельную задачу с громоздкими математическими выкладками.

4.5 Вычисление матрицы градиентов четырехугольного элемента с восемью узлами

Для численного интегрирования объемного интеграла в соотношении (2) предлагается [7] переходить к локальной системе координат (рисунок 5). Здесь замену переменных интегрирования можно сделать с помощью соотношения (3):

(3)

где – единичная толщина элемента, а – определитель матрицы преобразования координат Якоби. Тогда объемный интеграл в (2) приводится к виду (4):

(4)

Переход к переменным интегрирования и упрощает пределы интегрирования, что позволяет разработать единый алгоритм вычисления интегралов через квадратуру Гаусса-Лежандра.

Рисунок 5 – Система координат для линейного прямоугольного элемента

4.6 Алгоритм определения деформации и напряжений

Выражение перепишем в следующем виде (5):

(5)

где

,

, ,

,

, (6)

,

, .

Теперь из выражения (5) компоненту тензора деформации определим следующим образом:

(7)

или учитывая равенства (5)

(8)

В выражении (8) координаты и относятся к линейному прямоугольнику (рисунок 5), и им при рассмотрении восьмиугольных прямоугольников будут соответствовать и .

Аналогично определяются и остальные компоненты тензора деформации. Из формулы (5) для компоненты имеем:

или, учитывая (4) и , получим:

(9)

Деформация сдвига определяется аналогичным образом.

Для определения компонентов напряжений воспользуемся формулами вычисления напряжений [7]:

, и , (10)

где, для плоской деформации ;

- модуль упругости;

- коэффициент Пуассона.

Программа численной реализации приведенного алгоритма вычисления компонентов вектора деформации и вектора напряжений приведена в строках 412-527 приложения А настоящего документа. Кроме этого, в приложении А предусмотрены фрагменты программ построения графиков горизонтальных перемещений точек характерных вертикальных сечений (строки 297-312 приложения А), вертикальных перемещений точек горизонтальных сечений, расположенных на различных глубинах (строки 316-408 приложения А) и значений компонент векторов деформаций и напряжений (строки 441-475, 529-602 приложения А).

4.7 Краткая инструкция по использованию программы

BASIC_NDS_MKE_8_uzlov

Целью создания настоящего нормативного документа является разработка расчетной программы на языке MATLAB.

Расчетная программа предназначена для вычисления значений компонентов перемещений, деформаций и напряжений в точках исследуемой дорожной конструкции (рисунок 4).

В мировой практике существуют и другие пакеты программ, способных решать задачу о напряженно-деформированном состоянии в твердом деформируемом теле. Наиболее известным среди них является программный комплекс ANSYS. Однако, все они являются коммерческими, алгоритм их реализации скрытый, и простому пользователю использовать их для расчетов своих задач, тем более их видоизменить просто невозможно. Настоящий программный продукт BASIC_NDS_MKE_8_uzlovявляется открытым, алгоритм его реализации расписан открытым текстом, программа легко адаптируется и всегда есть возможность ее усовершенствовать. Например, ее можно дополнить отдельными блоками и решить задачи об условиях образования низкотемпературных трещин в дорожном покрытии, усталостных трещин в конструктивных слоях дорожной одежды.

Для запуска программы BASIC_NDS_MKE_8_uzlovдля решения конкретной задачи,в первую очередь необходимо, чтобы в компьютере была установлена одна из версии программного комплекса MATLAB,и в директорию этого комплекса был внесен текст (листинг) программы BASIC_NDS_MKE_8_uzlov. На рабочий стол компьютера выносится ярлык программного комплекса MATLAB (рисунок 6).

Рисунок 6 – Ярлык программного комплекса MATLAB

В приложении А настоящего документа приведен исходный код программы BASIC_NDS_MKE_8_uzlov. Для облегчения обращения к тексту листинг снабжен построчной нумерацией, которую нужно удалить перед установкой программы в директорию программного комплекса MATLAB.

Вызов программного комплекса MATLAB для расчетов осуществляется двойным щелчком на ярлыке (рисунок 6).

На рабочем столе появится главная страница MATLAB (рисунок 7).

Рисунок 7 – Вид главной страницы программы наMATLAB

На рисунке 7 на главной странице MATLAB приведен фрагмент программы BASIC_NDS_MKE_8_uzlov. Если же этот фрагмент будет отсутствовать, а на главной странице присутствует фрагмент другой программы, то поиск необходимой программы BASIC_NDS_MKE_8_uzlov будет произведен простым нажатием на кнопку OPEN, расположенной в верхней части главной страницы. После нажатия кнопки OPEN получим список программ, находящихся в директории MATLAB (рисунок 8).

Рисунок 8 – Вид списка программ в директории

Простым нажатием на название нужной программы BASIC_NDS_MKE_8_uzlovиз приведенного списка вызывается эта программа на главную страницу (рисунок 7).

Рисунок 9 – Строка 109-110 для введения физико-механических характеристик материалов слоев

По умолчанию в программе установлены физико-механические характеристики автомобильной дороги.Для изменения их значений данные необходимо ввести в строки 109-110 (рисунок 9). После введения физико-механических характеристик вводим значение прилагаемой вертикальной силы на поверхность в строку 150,а так же вводим номера поверхностных узлов куда необходимо приложить нагрузку (рисунок10).

Рисунок 10 – Строка для введения значения прилагаемой силы на поверхность и номеров узлов приложения сил

Следующим этапом при решении задачи является получение графических результатов по необходимым для пользователя требованиям. Ниже приведены примеры запроса графиков вертикальных перемещений узлов на разных глубинах (рисунки 11-17).

Рисунок 11 – Получение графиков вертикальных перемещений узлов на глубине h= 5 см

Рисунок 12 – Получение графиков вертикальных перемещений узлов на глубине h= 15 см

Рисунок 13 – Получение графиков вертикальных перемещений узлов на глубине h= 25 см

Рисунок 14 – Получение графиков вертикальных перемещений узлов на глубине h= 45 см

Рисунок 15 – Получение графиков компонент напряжений SigmaX

Рисунок 16 – Получение графиков компонент напряжений SigmaY

Рисунок 17 – Получение графиков компонент напряжений TauXY

Просмотрев вес текст программы на главной странице, и установив все необходимые исходные данные, а при необходимости, вносив изменения в тело программы можно будет запустить программу на выполнение. Для этого достаточно нажать на кнопкуRun (рисунок 18).

Рисунок 18 – Кнопка Run на панели инструментов для запуска программы

В тексте программы присутствуют команды, согласно которым на экран будет выводиться необходимая информация о результатах решения задачи. Такая информация может иметь графическую форму, может быть в виде таблиц, или в виде текстового файла. Все это зависит от воли заказчика.

Ниже, на рисунках 13-18, в качестве примера приведены результаты решения тестовой задачи, где приведены картины формирования поля горизонтальных и вертикальных перемещений точек характерных сечений исследуемой области и компоненты напряжений в характерных горизонтальных сечениях. В тестовой задаче дорожная конструкция нагружена распределенной вертикальной нагрузкой от одиночного колеса автомобиля, действующей в центре поверхности дорожного покрытия.

4.8 Примеры расчетов по программе BASIC_NDS_MKE_8_uzlov

Ниже приводятся результаты компьютерной реализации программы BASIC_NDS_MKE_8_uzlovв виде графиков перемещений и напряжений.

1–горизонтальные смещения точек вертикального сечения, отстоящего от вертикальной оси симметрии слева на 30 см; 2–горизонтальные смещения точек вертикальной оси симметрии; 3–горизонтальные смещения точек вертикального сечения, отстоящего от вертикальной оси симметрии справа на 30 см

Рисунок 19 – Горизонтальные смещения точек вертикального сечения

Рисунок 20 – Вертикальные смещения точек горизонтальных сечений, расположенных на глубинах: 1– h=5 см; 2–h=15 см; 3–h=25 см; 4–h=45 см; 5–h=80 см.

Рисунок 21 – Распределение напряжений на глубинах h

Рисунок 22 – Распределение напряжений на глубинах h

Рисунок 23 – Распределение напряжений на глубинах h

Приложение А

(обязательное)

Исходный код программы BASIC_NDS_MKE_8_uzlov

1	% БАЗОВАЯ ВЕРСИЯ ПРОГРАММЫ BASIC_NDS_MKE_8_uzlov
2	% РЕШЕНИЕ ТЕСТОВОЙ ЗАДАЧИ О НДС С ПОМОЩЬЮ
3	%ЧЕТЫРЕХУГОЛЬНОГО
4	% КВАДРАТИЧНОГО ЭЛЕМЕНТА С 8-Ю УЗЛАМИ МЕТОДОМ КОНЕЧНЫХ
5	% ЭЛЕМЕНТОВ
6
7	% Начало разработки 8.07.2017 г
8	% Окончание 1.09.2018 г
9
10	Tic
11
12	ever=23; egor=60; % количество конечных элементов в вертикальном (ever) и
13	% горизонтальном (egor) рядах элементов в конечно-элементной
14	% разбивке
15	np=(3*ever+2)*egor+2*ever+1; % общее число узлов
16	x=zeros(np); y=zeros(np);
17
18	% ПОРЯДОК НУМЕРАЦИИ ЭЛЕМЕНТОВ - СЛЕВА НАПРАВО (РЯДЫ % НАЧИНАЮТСЯ СНИЗУ)
19	% ПОРЯДОК НУМЕРАЦИИ УЗЛОВ - СЛЕВА НАПРАВО. РЯДЫ УЗЛОВ ТАКЖЕ % НАЧИНАЮТСЯ
20	% СНИЗУ. ПРОМЕЖУТОЧНЫЕ РЯДЫ УЗЛОВ, ПРОХОДЯЩИЕ ЧЕРЕЗ СЕРЕДИНУ % ЭЛЕМЕНТОВ
21	% СОДЕРЖАТ egor+1 УЗЛОВ
22
23	% ФОРМИРОВАНИЕ МАССИВА НОМЕРОВ УЗЛОВ С ЗАДАННЫМИ ЗНАЧЕНИЯМИ
24	% ПЕРЕМЕЩЕНИЙ ДЛЯ ЗАДАНИЯ ГРАНИЧНЫХ УСЛОВИЙ В ПЕРЕМЕЩЕНИЯХ
25
26	Mz=zeros(400);
27	np=ever*(3*egor+2)+2*egor+1;
28	for i=1:(2*egor+1)
29	Mz(i)=i;
30	End
31	k=2*egor+2;
32	for i=1:ever
33	Mz(k)=(2*egor+1)+1+(i-1)*(3*egor+2);
34	Mz(k+1)=Mz(k)+egor;
35	Mz(k+2)=Mz(k+1)+1;
36	Mz(k+3)=Mz(k+2)+2*egor;
37	k=k+4;
38	End
39	for i=1:(2*egor+1)
40	Mz(k)=Mz(i)+np;
41	k=k+1;
42	End
43
44	% МАССИВЫ ШАГОВ ПО КООРДИНАТНЫМ ОСЯМ В МЕТРАХ
45	format long
46	dx=zeros(2*egor);
47	dy=zeros(2*ever);
48
49	dx=[
50	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
51	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
52	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
53	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
54	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
55	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
56	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
57	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
58	0.250e-01 0.250e-01 0.250e-01 0.250e-01];
59
60	dy=[0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 ...
61	0.1000 0.1000 0.1000 0.1000 0.1000 0.0750 0.0750 ...
62	0.0625 0.0625 0.0500 0.0500 0.0500 0.0500 0.0250 ...
63	0.0250 0.0250];
64
65	yy=[0.000 0.200 0.400 0.600 0.800 1.000 1.200 1.400 1.600 1.800 2.000 2.200 ...
66	2.400 2.550 2.700 2.825 2.950 3.050 3.150 3.250 3.350 3.400 3.450 3.500];
67
68	% БЛОК ФОРМИРОВАНИЯ КООРДИНАТ
69	% Положительное направление оси х - слева направо
70	% Положительное направление оси у - снизу вверх
71	k=1;
72	NN=zeros(ever);
73	for n=1:ever
74	for m=1:egor
75	n1=2*m-1+(3*egor+2)*(n-1); n2=n1+1; n3=n1+2;
76	n4=n3+2*egor-m+1; n5=n3+(3*egor+2);
77	n6=n5-1; n7=n5-2; n8=n4-1;
78
79	x(n1)=(m-1)*2*dx(m); x(n2)=x(n1)+dx(m); x(n3)=x(n1)+2*dx(m);
80	x(n4)=x(n3); x(n5)=x(n3); x(n6)=x(n2); x(n7)=x(n1); x(n8)=x(n1);
81
82	y(n1)=yy(n); y(n2)=y(n1); y(n3)=y(n1); y(n4)=y(n3)+dy(n);
83	y(n5)=y(n3)+2*dy(n); y(n6)=y(n5); y(n7)=y(n5); y(n8)=y(n4);
84	if(m==36)
85	NN(2*k-1)=n3;
86	NN(2*k)=n4;
87	k=k+1;
88	end
89	End
90	End
91	NN(1:2*ever+1);
92
93	% ИСХОДНЫЕ ДАННЫЕ И ОБНУЛЕНИЕ МАТРИЦ
94	np2=2*np;
95	K=zeros(np2,np2);
96	U=zeros(np2+1);
97	U0=zeros(np2+1);
98	F=zeros(np2);
99
100	% НАИМЕНОВАНИЕ МАТЕРИАЛОВ КОНСТРУКТИВНЫХ СЛОЕВ (СВЕРХУ ВНИЗ)
101	% 1. АБ МЗ ПЛОТНЫЙ
102	% 2. АБ КЗ ПОРИСТЫЙ
103	% 3. ЩЕБЕНОЧНО-ПЕСЧАНАЯ СМЕСЬ 8%
104	% 4. ЩЕБЕНОЧНАЯ СМЕСЬ
105	% 5. ГРАВИЙНО-ПЕСЧАНАЯ СМЕСЬ
106	% 6. ГРУНТ СУГЛИНОК ЛЕГКИЙ
107
108	% ФИЗИКО-МЕХАНИЧЕСКИЕ ХАРАКТЕРИСТИКИ МАТЕРИАЛОВ СЛОЕВ(СВЕРХУ % ВНИЗ)
109	E=[6933e+05 4142e+05 1000e+05 250e+05 180e+05 58e+05]; % Модуль упр.кг/м^2
110	nu=[0.18 0.27 0.30 0.30 0.30 0.35]; % КОЭФФИЦИЕНТЫ ПУАССОНА
111
112	% ГРАНИЧНЫЕ УСЛОВИЯ В НАГРУЗКАХ
113
114	% Нагрузка от одиночного колеса
115
116	% ДАННЫЕ ИЗ ИНТЕРНЕТА
117	% 1 Па=0.102 кгс/м^2
118	% 0.6 МПа=6.12 кгс/см^2=60000.0 кгс/м^2
119
120	% РАСПРЕДЕЛЕННАЯ ПО УГЛОВЫМ УЗЛАМ ЭЛЕМЕНТА НАГРУЗКА ОТ % ОДИНОЧНОГО КОЛЕСА
121	% ПРИКЛАДЫВАЕТСЯ ПО ЦЕНТРУ РАСЧЕТНОЙ СХЕМЫ. При этом промежуточные % узлы
122	% игнорируются
123
124	% НАГРУЗКА НА КОЛЕСО РАССЧИТЫВАЕТСЯ ИЗ РАСЧЕТА НОРМАТИВНОЙ % ИНТЕНСИВНОСТИ
125	% НАГРУЗКИ q=0.6 MPa
126	% ЕСЛИ 1 Па=0.102 кгс/м^2 (из Интернета) ТО
127	% 0.6 МПа=6.12 кгс/см^2 = 60000.0 кгс/м^2.
128	% Таким образом, если на 1 погонный метр полотна приходится 60000 кг веса, то
129	% на одно расстояние 0.05 м – 3000 кг веса.
130	% Если ширину колеса принять равной 40 см, а ширину одного конечного
131	% элемента в горизонтальном направлении – 5 см, то транспортная нагрузка
132	% будет распределена по 9 узлам в пределах 8 конечных элементов.
133	% Расстояние между угловыми узлами - 0.05 м (5 cм - ширина конечного элемента).
134	% При этом надо учесть, что
135	% два крайних узла (1-ый и 9-ый узлы) будут нагружены наполовину.
136	% В таком случае, нагрузки в узловых точках будут распределены
137	% следующим образом:
138	% P1=P9=1500 кг, P2=P3= … =P8=3000 кг
139
140	% РАСЧЕТ НАГРУЗКИ НА УЗЕЛ
141	% РАССТОЯНИЕ МЕЖДУ УЗЛАМИ 0.025 М (учитываются и промежуточные узлы)
142	% НА 1 М ПРИХОДИТСЯ 60000 КГ ВЕСА, А НА ОДНО РАССТОЯНИЕ 0.025 М
143	% ПРИХОДИТСЯ 1500 КГ ВЕСА
144
145	% РЕЗУЛЬТАТЫ БУДУТ ПОЛУЧЕНЫ В КГ И МЕТРАХ
146
147	% ВЕЛИЧИНЫ ВЕРТИКАЛЬНЫХ СИЛ, ПРИЛОЖЕННЫХ ПО УГЛОВЫМ УЗЛАМ
148	% В ЦЕНТРАЛЬНОЙ ЧАСТИ ЭПЮРЫ
149
150	P=0.3e+04; % кг
151
152	% НОМЕРА УГЛОВЫХ УЗЛОВ ПРИЛОЖЕНИЯ СИЛ
153
154	nF(1:9)=[4239 4241 4243 4245 4247 4249 4251 4253 4255];
155	for i=1:9
156	j=nF(i);
157	xP(i)=x(j);
158	end
159
160	% ПРИЛОЖЕНИЕ СИЛ К УГЛОВЫМ УЗЛАМ
161	F(nF(1)+np)=-P/2;
162	F(nF(9)+np)=-P/2;
163	for i=1:7
164	F(nF(1+i)+np)=-P;
165	end
166
167	% ФОРМИРОВАНИЕ МАТРИЦЫ ЖЕСТКОСТИ ЭЛЕМЕНТА [ke] И
168	% ГЛОБАЛЬНОЙ МАТРИЦЫ ЖЕСТКОСТИ СИСТЕМЫ [K]
169
170	mpe=[-0.577350 0.577350]; % КООРДИНАТЫ ТОЧЕК ИНТЕГРИРОВАНИЯ
171	mH=[1 1]; % ВЕСОВЫЕ КОЭФФИЦИЕНТЫ
172	K=zeros(np2,np2);
173	for n=1:ever
174	for m=1:egor
175	n1=2*m-1+(3*egor+2)*(n-1); n2=n1+1; n3=n1+2;
176	n4=n3+2*egor-m+1; n5=n3+(3*egor+2);
177	n6=n5-1; n7=n5-2; n8=n4-1;
178
179	X1=x(n1); Y1=y(n1);
180	X2=x(n3); Y2=y(n3);
181	X3=x(n5); Y3=y(n5);
182	X4=x(n7); Y4=y(n7);
183
184	% ФОРМИРОВАНИЕ МАТРИЦЫ УПРУГОСТИ ДЛЯ СЛУЧАЯ ПЛОСКОЙ % ДЕФОРМАЦИИ
185	if((n>=1)&(n<=7))
186	e1=E(6); nu1=nu(6); % ГРУНТ СУГЛИНОК ЛЕГКИЙ
187	end
188	if((n>=8)&(n<=12))
189	e1=E(5); nu1=nu(5); % ГРАВИЙНО-ПЕСЧАНАЯ СМЕСЬ
190	end
191	if((n>=13)&(n<=16))
192	e1=E(4); nu1=nu(4); % ЩЕБЕНОЧНАЯ СМЕСЬ
193	end
194	if((n>=17)&(n<=20))
195	e1=E(3); nu1=nu(3); % ЩЕБЕНОЧНО-ПЕСЧАНАЯ СМЕСЬ 8%
196	end
197	if((n>=21)&(n<=22))
198	e1=E(2); nu1=nu(2); % АБ КЗ ПОРИСТЫЙ
199	end
200	if(n==23)
201	e1=E(1); nu1=nu(1); % АБ МЗ ПЛОТНЫЙ
202	end
203
204	D=zeros(3,3);
205	D(3,3)=e1/(2*(1+nu1)); D(2,2)=2*D(3,3)*(1-nu1)/(1-2*nu1);
206	D(1,1)=D(2,2); D(1,2)=2*D(3,3)*nu1/(1-2*nu1);
207	D(2,1)=D(1,2);
208
209	dJ=abs((X1-X2)*(Y2-Y4))/4; % якобиан
210
211	% ИСПОЛЬЗОВАНИЕ КВАДРАТУРЫ ГАУССА ДЛЯ ЧИСЛЕННОГО % ИНТЕГРИРОВАНИЯ
212	% ДВОЙНОГО ИНТЕГРАЛА В МАТРИЦЕ ЖЕСТКОСТИ ЭЛЕМЕНТА [Ke]
213
214	% BtDBnds1(psi,eta,D,X1,X2,Y1,Y3) - подпрограмма вычисления
215	% матрицы жесткости элемента
216	Ke=zeros(16,16);
217	psi=mpe(1); eta=mpe(1);
218	k11=BtDBnds1(psi,eta,D,X1,X2,Y1,Y3);
219	psi=mpe(1); eta=mpe(2);
220	k12=BtDBnds1(psi,eta,D,X1,X2,Y1,Y3);
221	psi=mpe(2); eta=mpe(1);
222	k21=BtDBnds1(psi,eta,D,X1,X2,Y1,Y3);
223	psi=mpe(2); eta=mpe(2);
224	k22=BtDBnds1(psi,eta,D,X1,X2,Y1,Y3);
225	Ke=dJ*(mH(1)*(mH(1)*k11+mH(2)*k12)+mH(2)*(mH(1)*k21+mH(2)*k22));
226
227	ne8=[n1 n2 n3 n4 n5 n6 n7 n8];
228	for i=1:8
229	for j=1:8
230	K(ne8(i),ne8(j))=K(ne8(i),ne8(j))+Ke(i,j);
231	K(ne8(i),ne8(j)+np)=K(ne8(i),ne8(j)+np)+Ke(i,j+8);
232	K(ne8(i)+np,ne8(j))=K(ne8(i)+np,ne8(j))+Ke(i+8,j);
233	K(ne8(i)+np,ne8(j)+np)=K(ne8(i)+np,ne8(j)+np)+Ke(i+8,j+8);
234	end
235	end
236	end
237	end
238
239	disp('Матрица К сформирована')
240
241	% ПРЕОБРАЗОВАНИЕ СЛАУ
242
243	np2=2*np;
244	i=1;
245	for n=1:np2
246	if(n==Mz(i))
247	F(n)=K(n,n)*U(n);
248	for m=1:np2
249	if(m~=n)
250	K(n,m)=0.0;
251	end
252	end
253	i=i+1;
254	end
255	end
256
257	i=1;
258	for n=1:np2
259	if(n==Mz(i))
260	for m=1:np2
261	if(m~=n)
262	F(m)=F(m)-K(m,n)*U(n);
263	K(m,n)=0.0;
264	end
265	end
266	i=i+1;
267	end
268	end
269
270	% РЕШЕНИЕ СЛАУ
271	U=K\F;
272	disp('Система [K]{U}={F} решена')
273
274	% Номера узлов характерных вертикальных сечений
275	nom12=[
276	49 146 231 328 413 510 595 692 ...
277	777 874 959 1056 1141 1238 1323 1420 ...
278	1505 1602 1687 1784 1869 1966 2051 2148 ...
279	2233 2330 2415 2512 2597 2694 2779 2876 ...
280	2961 3058 3143 3240 3325 3422 3507 3604 ...
281	3689 3786 3871 3968 4053 4150 4235];
282	nom15=[
283	61 152 243 334 425 516 607 698 ...
284	789 880 971 1062 1153 1244 1335 1426 ...
285	1517 1608 1699 1790 1881 1972 2063 2154 ...
286	2245 2336 2427 2518 2609 2700 2791 2882 ...
287	2973 3064 3155 3246 3337 3428 3519 3610 ...
288	3701 3792 3883 3974 4065 4156 4247];
289	nom18=[
290	73 158 255 340 437 522 619 704 ...
291	801 886 983 1068 1165 1250 1347 1432 ...
292	1529 1614 1711 1796 1893 1978 2075 2160 ...
293	2257 2342 2439 2524 2621 2706 2803 2888 ...
294	2985 3070 3167 3252 3349 3434 3531 3616 ...
295	3713 3798 3895 3980 4077 4162 4259];
296
297	% ГРАФИКИ ГОРИЗОНТАЛЬНЫХ ПЕРЕМЕЩЕНИЙ УЗЛОВ В СРЕДНЕМ СЕЧЕНИИ
298	Y=zeros(2*ever+1);
299	X150=zeros(2*ever+1);
300	X120=zeros(2*ever+1);
301	X180=zeros(2*ever+1);
302	for i=1:2*ever+1
303	j120=nom12(i);
304	j150=nom15(i);
305	j180=nom18(i);
306	Y(i)=y(j150);
307	X120(i)=U(j120)-0.0001; % СЛЕВА ОТ ВЕРТИКАЛЬНОЙ ОСИ СИММЕТРИИ НА 30 СМ
308	X150(i)=U(j150)*10; % НА ВЕРТИКАЛЬНОЙ ОСИ СИММЕТРИИ
309	X180(i)=U(j180)+0.0001; % СПРАВА ОТ ВЕРТИКАЛЬНОЙ ОСИ СИММЕТРИИ НА % 30 СМ
310	end
311	hPlot=plot(X150,Y,'b-',X120,Y,'r-',X180,Y,'k-.');grid on;xlabel('Ux, m');ylabel('y, m')
312	set(hPlot,'LineWidth',2);
313
314	figure
315
316	% ФОРМИРОВАНИЕ МАССИВА НОМЕРОВ ГОРИЗОНТАЛЬНЫХ РЯДОВ УЗЛОВ % Nv(I,m)
317
318	% ОБЩЕЕ ЧИСЛО УЗЛОВ В КАЖДОМ ПОЛНОМ ГОРИЗОНТАЛЬНОМ РЯДУ РАВНО % 2*egor+1,
319	% А ЧИСЛО ТАКИХ ГОРИЗОНТАЛЬНЫХ РЯДОВ РАВНО ever+1
320	% СЕЧЕНИЯ, ПРОХОДЯЩИЕ ЧЕРЕЗ СЕРЕДИНУ ЭЛЕМЕНТОВ, НЕ % РАССМАТРИВАЮТСЯ.
321	% ОНИ ИСПОЛЬЗУЕТСЯ ДЛЯ ПОСТРОЕНИЯ ГРАФИКОВ ИСКОМЫХ ВЕЛИЧИН НА % КАЖДОМ
322	% ГОРИЗОНТАЛЬНОМ УРОВНЕ (ГЛУБИНЕ)
323
324	Nv=zeros(ever+1,2*egor+1);
325	for n=1:(ever+1)
326	for m=1:(2*egor+1)
327	I=2*n-1;
328	Nv(I,m)=m+(3*egor+2)*(n-1);
329	end
330	end
331	for n=1:ever
332	for m=1:(egor+1)
333	I=2*n;
334	Nv(I,2*m-1)=m+(2*egor+1)*n+(egor+1)*(n-1);
335	end
336	end
337
338	format short
339	Set=zeros(2*ever+1,4);
340	for n=1:2*ever+1
341	i=Nv(n,1);
342	Set(n,1)=n;
343	Set(n,2)=i;
344	Set(n,3)=y(i);
345	Set(n,4)=3.50-y(i);
346	end
347	for m=1:47
348	% fprintf('%6u %6u %10.4f %10.4f\n',Set(m,1),Set(m,2),Set(m,3),Set(m,4));
349	end
350
351	% ГРАФИКИ ВЕРТИКАЛЬНЫХ ПЕРЕМЕЩЕНИЙ УЗЛОВ
352	% ГООРИЗОНТАЛЬНЫХ СЕЧЕНИЙ НА ГЛУБИНАХ h
353	X=zeros(2*egor+1);
354	Y=zeros(2*egor+1);
355	for i=1:(2*egor+1)
356	j=Nv(45,i); % ГЛУБИНА h=5 СМ
357	X(i)=x(j);
358	Y(i)=U(j+np);
359	end
360	hPlot=plot(X,Y,'-');grid on;xlabel('X, m');ylabel('Uy, m')
361	set(hPlot,'LineWidth',2);
362	hold on
363
364	X=zeros(2*egor+1);
365	Y=zeros(2*egor+1);
366	for i=1:(2*egor+1)
367	j=Nv(41,i); % ГЛУБИНА h=15 СМ
368	X(i)=x(j);
369	Y(i)=U(j+np);
370	end
371	hPlot=plot(X,Y,'-');grid on;xlabel('X, m');ylabel('Uy, m')
372	set(hPlot,'LineWidth',2);
373	hold on
374
375	X=zeros(2*egor+1);
376	Y=zeros(2*egor+1);
377	for i=1:(2*egor+1)
378	j=Nv(39,i); % ГЛУБИНА h=25 СМ
379	X(i)=x(j);
380	Y(i)=U(j+np);
381	Ux325(i)=U(j);
382	Uy325(i)=Y(i);
383	end
384	hPlot=plot(X,Y,'--');grid on;xlabel('X, m');ylabel('Uy, m')
385	set(hPlot,'LineWidth',2);
386	hold on
387
388	for i=1:(2*egor+1)
389	j=Nv(35,i); % ГЛУБИНА h=45 СМ
390	X(i)=x(j);
391	Y(i)=U(j+np);
392	X345(i)=X(i);
393	Ux305(i)=U(j);
394	Uy305(i)=Y(i);
395	end
396	hPlot=plot(X,Y,'--');grid on;xlabel('X, m');ylabel('Uy, m')
397	set(hPlot,'LineWidth',2);
398	hold on
399
400	for i=1:(2*egor+1)
401	j=Nv(29,i); % ГЛУБИНА h=80 СМ
402	X(i)=x(j);
403	Y(i)=U(j+np);
404	Ux270(i)=U(j);
405	Uy270(i)=Y(i);
406	end
407	hPlot=plot(X,Y,'-.');grid on;xlabel('X, m');ylabel('Uy, m')
408	set(hPlot,'LineWidth',2);
409
410	figure
411
412	% ВЫЧИСЛЕНИЕ КОМПОНЕНТОВ ДЕФОРМАЦИИ В СРЕДИННЫХ ТОЧКАХ
413	% ЭЛЕМЕНТОВ (eps=0.0, eta=0.0)
414
415	epsX0=zeros(ever*egor);
416	epsY0=zeros(ever*egor);
417	gamXY0=zeros(ever*egor);
418
419	for n=1:ever
420	for m=1:egor
421	n1=2*m-1+(3*egor+2)*(n-1); n2=n1+1; n3=n1+2;
422	n4=n3+2*egor-m+1; n5=n3+(3*egor+2);
423	n6=n5-1; n7=n5-2; n8=n4-1;
424
425	dJ=abs((x(n1)-x(n3))*(y(n3)-y(n7)))/4; % якобиан
426	i=m+egor*(n-1);
427	epsX0(i)=(U(n4)-U(n8))/(x(n4)-x(n8));
428	epsY0(i)=(-U(n2+np)+U(n6+np))/(y(n6)-y(n2));
429	gamXY0(i)=-(-U(n2)+U(n6))/(x(n4)-x(n8))+(U(n4+np)-U(n8+np))/(y(n6)-y(n2));
430	end
431	end
432
433	% НОМЕРА ЭЛЕМЕНТОВ ПОСТРОЧНО
434	NEv=zeros(ever,egor);
435	for n=1:ever
436	for m=1:egor
437	i=m+egor*(n-1);
438	NEv(n,m)=i;
439	end
440	end
441
442	% НОМЕР СРЕДИННОГО УЗЛА, ДЛЯ КОТОРОГО ВЫЧИСЛЯЮТСЯ ДЕФОРМАЦИИ % И НАПРЯЖЕНИЯ,
443	% СОВПАДАЕТ С НОМЕРОМ ЭЛЕМЕНТА, РАСПОЛОЖЕННОГО В МАССИВЕ NEv(n,m)
444
445	X=zeros(egor); Y25=zeros(egor); Y200=zeros(egor); Y400=zeros(egor);
446	for m=1:egor
447	n=Nv(1,2*m);
448	X(m)=x(n);
449	i25=NEv(23,m);
450	i200=NEv(20,m);
451	i400=NEv(18,m);
452	Y25(m)=epsX0(i25); % 2,5 cм
453	Y200(m)=epsX0(i200); % 20,0 cм
454	Y400(m)=epsX0(i400); % 40,0 см
455	end
456	hPlot=plot(X,Y25,'-',X,Y200,'--',X,Y400,'-.');...
457	grid on;xlabel('X, m');ylabel('epsX')
458	set(hPlot,'LineWidth',2);
459	figure
460
461	X=zeros(egor); Y25=zeros(egor); Y200=zeros(egor); Y400=zeros(egor);
462	for m=1:egor
463	n=Nv(1,2*m);
464	X(m)=x(n);
465	i25=NEv(23,m);
466	i200=NEv(20,m);
467	i400=NEv(18,m);
468	Y25(m)=epsY0(i25); % 2,5 cм
469	Y200(m)=epsY0(i200); % 20,0 cм
470	Y400(m)=epsY0(i400); % 40,0 см
471	end
472	hPlot=plot(X,Y25,'-',X,Y200,'--',X,Y400,'-.');...
473	grid on;xlabel('X, m');ylabel('epsY')
474	set(hPlot,'LineWidth',2);
475	figure
476
477	% ВЫЧИСЛЕНИЕ КОМПОНЕНТОВ НАПРЯЖЕНИЙ Sigma В СРЕДИННЫХ ТОЧКАХ
478	% ЭЛЕМЕНТОВ (eps=0.0, eta=0.0)
479	SigX=zeros(ever); SigY=zeros(ever);
480	k=0;
481	for n=1:ever
482	for m=1:egor
483	n1=2*m-1+(3*egor+2)*(n-1); n2=n1+1; n3=n1+2;
484	n4=n3+2*egor-m+1; n5=n3+(3*egor+2);
485	n6=n5-1; n7=n5-2; n8=n4-1;
486
487	% ФОРМИРОВАНИЕ МАТРИЦЫ УПРУГОСТИ ДЛЯ СЛУЧАЯ ПЛОСКОЙ % ДЕФОРМАЦИИ
488	if((n>=1)&(n<=7))
489	e1=E(6); nu1=nu(6); % ГРУНТ СУГЛИНОК ЛЕГКИЙ
490	end
491	if((n>=8)&(n<=12))
492	e1=E(5); nu1=nu(5); % ГРАВИЙНО-ПЕСЧАНАЯ СМЕСЬ
493	end
494	if((n>=13)&(n<=16))
495	e1=E(4); nu1=nu(4); % ЩЕБЕНОЧНАЯ СМЕСЬ
496	end
497	if((n>=17)&(n<=20))
498	e1=E(3); nu1=nu(3); % ЩЕБЕНОЧНО-ПЕСЧАНАЯ СМЕСЬ 8%
499	end
500	if((n>=21)&(n<=22))
501	e1=E(2); nu1=nu(2); % АБ КЗ ПОРИСТЫЙ
502	end
503	if(n==23)
504	e1=E(1); nu1=nu(1); % АБ МЗ ПЛОТНЫЙ
505	end
506
507	D=zeros(3,3);
508	D(3,3)=e1/(2*(1+nu1)); D(2,2)=2*D(3,3)*(1-nu1)/(1-2*nu1);
509	D(1,1)=D(2,2); D(1,2)=2*D(3,3)*nu1/(1-2*nu1);
510	D(2,1)=D(1,2);
511
512	X1=x(n1); X2=x(n3); Y2=y(n3); Y4=y(n7);
513	dJ=abs((X1-X2)*(Y2-Y4))/4; % Якобиан
514	G=e1/(2*(1+nu1));
515
516	i=m+egor*(n-1);
517	% В МПа
518	SigmaX0(i)=(D(1,1)*epsX0(i)+D(1,2)*epsY0(i))/100000;
519	SigmaY0(i)=(D(2,1)*epsX0(i)+D(2,2)*epsY0(i))/100000;
520	TauXY0(i)=D(3,3)*gamXY0(i)/100000;
521	if(m==30)
522	k=k+1;
523	SigX(n)=SigmaX0(i);
524	SigY(n)=SigmaY0(i);
525	end
526	end
527	end
528
529	% ПОСТРОЕНИЕ ГРАФИКОВ КОМПОНЕНТ НАПРЯЖЕНИЙ SigmaX0, SigmaY0 И
530	% TauXY0 В СРЕДИННЫХ УЗЛАХ ЭЛЕМЕНТОВ
531
532	X=zeros(egor); Y23=zeros(egor); Y22=zeros(egor); Y20=zeros(egor);
533	Y18=zeros(egor); Y14=zeros(egor); Y12=zeros(egor);
534	for m=1:egor
535	n=Nv(1,2*m);
536	X(m)=x(n);
537	i23=NEv(23,m); % h=2,5 см
538	i22=NEv(22,m); % h=7,5 см
539	i20=NEv(20,m); % h=20,0 см
540	i18=NEv(18,m); % h=40,0 см
541	i14=NEv(14,m); % h=87,5 см
542	i12=NEv(12,m); % h=120,0 см
543	Y23(m)=SigmaX0(i23); % 2,5 cм
544	Y22(m)=SigmaX0(i22); % 7,5 cм
545	Y20(m)=SigmaX0(i20); % 20,0 cм
546	Y18(m)=SigmaX0(i18); % 40,0 cм
547	Y14(m)=SigmaX0(i14); % 87,5 cм
548	Y12(m)=SigmaX0(i12); % 120,0 cм
549	end
550	hPlot=plot(X,Y22,'-',X,Y20,'--',X,Y18,'-.',X,Y14,'-r');...
551	grid on;xlabel('X, m');ylabel('SigmaX, МПа')
552	set(hPlot,'LineWidth',2);
553	figure
554
555	X=zeros(egor); Y23=zeros(egor); Y22=zeros(egor); Y20=zeros(egor);
556	Y18=zeros(egor); Y14=zeros(egor); Y12=zeros(egor);
557	for m=1:egor
558	n=Nv(1,2*m);
559	X(m)=x(n);
560	i23=NEv(23,m); % h=2,5 см
561	i22=NEv(22,m); % h=7,5 см
562	i20=NEv(20,m); % h=20,0 см
563	i18=NEv(18,m); % h=40,0 см
564	i14=NEv(14,m); % h=87,5 см
565	i12=NEv(12,m); % h=120,0 см
566	Y23(m)=SigmaY0(i23); % 2,5 cм
567	Y22(m)=SigmaY0(i22); % 7,5 cм
568	Y20(m)=SigmaY0(i20); % 20,0 cм
569	Y18(m)=SigmaY0(i18); % 40,0 cм
570	Y14(m)=SigmaY0(i14); % 87,5 cм
571	Y12(m)=SigmaY0(i12); % 120,0 cм
572	end
573	hPlot=plot(X,Y22,'-',X,Y20,'--',X,Y18,'-.',X,Y14,'-r');...
574	grid on;xlabel('X, m');ylabel('SigmaY, МПа')
575	set(hPlot,'LineWidth',2);
576	figure
577
578	X=zeros(egor); Y23=zeros(egor); Y22=zeros(egor); Y20=zeros(egor);
579	Y18=zeros(egor); Y14=zeros(egor); Y12=zeros(egor);
580	for m=1:egor
581	n=Nv(1,2*m);
582	X(m)=x(n);
583	i23=NEv(23,m); % h=2,5 см
584	i22=NEv(22,m); % h=7,5 см
585	i21=NEv(21,m); % h=12,5 см
586	i20=NEv(20,m); % h=20,0 см
587	i18=NEv(18,m); % h=40,0 см
588	i16=NEv(16,m); % h=61,25 см
589	i14=NEv(14,m); % h=87,5 см
590	i12=NEv(12,m); % h=120,0 см
591	Y23(m)=TauXY0(i23); % 2,5 cм
592	Y22(m)=TauXY0(i22); % 7,5 cм
593	Y21(m)=TauXY0(i21); % 12,5 cм
594	Y20(m)=TauXY0(i20); % 20,0 cм
595	Y18(m)=TauXY0(i18); % 40,0 cм
596	Y16(m)=TauXY0(i16); % 61,25 cм
597	Y14(m)=TauXY0(i14); % 87,5 cм
598	Y12(m)=TauXY0(i12); % 120,0 cм
599	end
600	hPlot=plot(X,Y22,'-',X,Y20,'--',X,Y14,'-.');...
601	grid on;xlabel('X, m');ylabel('TauXY, МПа')
602	set(hPlot,'LineWidth',2);
603
604	toc

Приложение Б

(обязательное)

Исходный код подпрограммы BtDBnds

1	function[BDB]=BtDBnds(psi,eta,D,X1,X2,Y1,Y3);
2	A=(-2/(X1-X2));
3	B1=A*(1-eta)*(2*psi+eta)/4; B2= A*(-psi*(1-eta));
4	B3= A*(1-eta)*(2*psi-eta)/4; B4= A*(1-eta^2)/2;
5	B5= A*(1+eta)*(2*psi+eta)/4; B6= A*(-psi*(1+eta));
6	B7= A*(1+eta)*(2*psi-eta)/4; B8= A*(-(1-eta^2)/2);
7	B=(-2/(Y1-Y3));
8	C1=B*(1-psi)*(psi+2*eta)/4; C2=-B*(1-psi^2)/2;
9	C3=B*(1+psi)*(-psi+2*eta)/4; C4=-B*eta*(1+psi);
10	C5=B*(1+psi)*(psi+2*eta)/4; C6=B*(1-psi^2)/2;
11	C7=B*(-(1-psi)*(psi-2*eta)/4); C8=-B*eta*(1-psi);
12
13
14	B=[B1 B2 B3 B4 B5 B6 B7 B8 0 0 0 0 0 0 0 0;
15	0 0 0 0 0 0 0 0 C1 C2 C3 C4 C5 C6 C7 C8;
16	C1 C2 C3 C4 C5 C6 C7 C8 B1 B2 B3 B4 B5 B6 B7 B8];
17	Bt=B';
18	BtD=Bt*D;
19	BDB=BtD*B;

Библиография

[1]Мартынов Н.Н., Иванов А.П. MATLAB 5.Х. Вычисления, визуализация, программирование. – М.:КУДИЦ-ОБРАЗ, 2000.-336 с.

[2] Коткин Г.Л., Черкасский В.С. Компьютерное моделирование физических процессов с использованием MATLAB:Учебное пособие/Новосиб. ун-т. Новосибирск, 2001. 173 с.

[3] Потемкин В.Г. Система инженерных и научных расчетов MATLAB 5.Х. В 2-х томах. –М.:ДИАЛОГ-МИФИ, 1999. 670 с.

[4] Тынкевич М.А. Численные методы.-Кемерово.:КузГТУ.1977.-122 с.

[5] Сегерлинд Л. Применение метода конечных элементов. М: Мир, 1979.– 392 с.

[6] СТ РК 1293-2004 Методы определения модуля упругости дорожных одежд нежесткого типа и их классификация.