ҚР Ұ 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 – зерттелетін аймақтың ені мен биіктігі; Р – көлік құралының дөңгелегінен түсетін жүктеме; ui, vi - х және у координаттық өстер бойымен бағытталған шекаралардағы нүктелердің жылжу компоненттері; 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 – нагрузка от колеса транспортного средства; ui, vi – компоненты перемещения точек на границах в направлении координатных осей 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 Методы определения модуля упругости дорожных одежд нежесткого типа и их классификация.