УДК 536.242:539.37 DOI 10.18698/0536-1044-2016-6-13-20
Численное моделирование температурно-структурного и напряженного состояний в процессе закалки железнодорожного рельса
А.М. Покровский, Ю.В. Воронов, Д.Н. Третьяков
МГТУ им. Н.Э. Баумана, 105005, Москва, Российская Федерация, 2-я Бауманская ул., д. 5, стр. 1
Numerical Simulation of Thermal-Structural and Stress States in the Process of Hardening Railway Rails
A.M. Pokrovskiy, Y.V. Voronov, D.N. Tretyakov
BMSTU, 105005, Moscow, Russian Federation, 2nd Baumanskaya St., Bldg. 5, Block 1 IQl e-mail: [email protected], [email protected], [email protected]
Изготовление высокопрочных железнодорожных рельсов — чрезвычайно важная проблема. В связи с этим разработка методов численного анализа температурно-структурного и напряженного состояния рельсов в процессе закалки, способствующих рационализации технологических режимов их изготовления, является актуальной задачей. Создана математическая модель, позволяющая описывать температурные поля, распределение структуры и термических напряжений в рельсе в течение всего процесса закалки. В основу решения нелинейной нестационарной задачи теплопроводности и термоупругопластичности положен метод конечных элементов. Для описания теплообмена использованы граничные условия третьего рода. Моделирование превращения аустенита в феррито-карбид в изотермических условиях проведено на основе уравнения Авраами. Переход от изотермической кинетики распада аусте-нита к неизотермическим условиям осуществлен согласно теории изокинетических реакций с привлечением правила аддитивности. Представлены результаты расчета температур, структур и напряжений в железнодорожном рельсе для различных моментов закалки. Показано, что головка рельса Р65 после закалки имеет структуру пластинчатого феррито-карбида. Мартенсит при закалке в масло присутствует только в структуре шейки и вблизи пера подошвы. Разработанные программные средства могут быть полезны для прогнозирования прочности рельса при эксплуатации.
Ключевые слова: железнодорожные рельсы, закалка, нелинейная нестационарная задача теплопроводности, кинетика структурных превращений, метод конечных элементов, термоупругопластичность, остаточные напряжения.
The manufacture of high-strength railroad rails is an extremely important issue. It necessitates the development of such methods for numerical analysis of thermal-structural and stress states of rails in the process of hardening that can be used to rationalize rail manufacturing processes. A mathematical model is created that can describe temperature fields, distribution patterns and thermal stresses in a rail during the entire hardening process. The finite element method is used as the basis for solving the nonlinear nonstationary problem of heat conduction and thermo-elastic-plasticity. Boundary conditions of the third kind are used to describe the heat transfer. Modelling of the transformation of austenite into ferrite-carbide in isothermal conditions is carried out using
the Avrami equation. The transition from the isothermal kinetics of austenite decomposition to the non-isothermal conditions is described by the theory of isokinetic reactions applying the additivity rule. The calculation results of the temperatures, structures and stresses in a railway rail at different stages of hardening are presented. It is shown that the head of the R65 rail after hardening has the structure of lamellar ferrite-carbide. When quenched in oil, martensite is present only in the structure of the neck and near the rail foot blade. The software developed by the authors can be used to predict the strength of the rail during operation.
Keywords: railroad rails, hardening, nonlinear nonstationary problem of heat conduction, kinetics of structural transformations, finite element method, thermos-elastic-plasticity, residual stresses.
Длина железнодорожных путей в России составляет около 10 % мировой протяженности. По отечественным железным дорогам перевозят порядка 40 % грузооборота этого вида транспорта в мире. Отечественные железные дороги отличает высокая интенсивность движения как пассажирских, так и товарных поездов. Загруженность железных дорог в России существенно выше, чем в США и европейских странах, где значительная часть перевозок по суше осуществляется автомобильным транспортом. Высокая интенсивность движения, постоянно возрастающие скорости поездов и массы перевозимых грузов предъявляют особые требования к качеству рельсов, являющихся основным элементом железнодорожного пути.
Для придания рельсам требуемых эксплуатационных свойств в процессе изготовления их подвергают термической обработке. При этом очень важно выбрать рациональные режимы термообработки, так как некорректные режимы могут привести к короблению рельсов, а также к возникновению и росту трещиноподобных дефектов вследствие высокого уровня остаточных термонапряжений. Экспериментальные методы исследования температурно-структурного и напряженного состояния рельсов при термической обработке малоэффективны, так как не позволяют проследить всю кинетику формирования структуры и остаточных напряжений, особенно во внутренних областях рельса.
Цель работы — поиск численных методов компьютерного моделирования, позволяющих отслеживать температурно-структурное и напряженное состояние рельса в течение всего процесса термической обработки.
В настоящей работе созданы программные средства и компьютерные методики для решения указанной задачи, которую можно разбить на три части: теплопроводность, моделирование структурного состава и расчет напряжений.
Поскольку температурные поля определяют кинетику структурных превращений и формирования термических напряжений, а теплофизиче-ские и физико-механические свойства зависят от температуры и структуры, все три задачи оказываются взаимосвязанными и решаются совместно. В расчете использован шаговый метод, при котором на каждом шаге по времени последовательно решают указанные задачи. В основу задач теплопроводности и расчета напряжений положен метод конечных элементов (МКЭ).
В качестве объекта исследования выбран железнодорожный рельс тяжелого типа Р65 [1] длиной 25 м, изготовленный из рельсовой стали 76ХФ (рис. 1). Основные размеры поперечного сечения рельса Р65, мм, приведены ниже:
Высота рельса Н....................................................180
Высота шейки к....................................................105
Ширина головки Ь..................................................75
Ширина подошвы В............................................150
Толщина шейки е....................................................18
Высота пера подошвы т....................................11,2
В качестве термической обработки использована объемная закалка, заключающаяся в сквозном нагреве рельса в печи до температуры 860 °С и последующем охлаждении в масляной ванне. Поскольку длина рельса во много раз превышает размеры поперечного сечения, задачу нелинейной нестационарной теплопроводности решали в плоской двумерной постановке.
Для изотропного тела при переменных теп-лофизических коэффициентах эту задачу можно описать дифференциальным уравнением [2]
дТ д f, дТ У д f дТ
. дУ.
ср— =—I À— I + — дт дх I дх J ду
\
+ qv, (1)
где с — коэффициент теплоемкости; р — плотность; Т(х, у, т) — температура, х, у — коорди-
В
Рис. 1. Поперечное сечение рельса: 1 — головка; 2 — шейка; 3 — подошва
наты, х — время; X — коэффициент теплопроводности; qv — мощность удельных источников энерговыделения.
Для описания условий теплообмена использованы граничные условия третьего рода [2]:
X ^ ^ = к [7с(х) - Гп(х)], (2)
где п — нормаль к поверхности; к — суммарный коэффициент теплоотдачи, учитывающий теплообмен конвекцией и излучением; Тс и Тп — температура окружающей среды и поверхности соответственно; индекс «п» относится к значениям на поверхности.
Краевая задача (1)-(2) проинтегрирована при начальном условии
Т (х, / ,0) = То, (3)
где То — начальная температура.
Решение задачи нестационарной теплопроводности для изотропного тела с использованием МКЭ при переменных теплофизических коэффициентах сводится к минимизации функционала, описывающего краевую задачу (1)-(3). Для ансамбля конечных элементов это приводит к следующему матричному уравнению [3]:
[С] М + [К ]{Т} = (7}, от
(4)
вектор-столбец температур в узлах конечно-элементной сетки; {7} — вектор-столбец тепловой нагрузки в узлах.
Формирование матриц [С],[К] и вектора {7} осуществляли, согласно МКЭ, посредством суммирования соответствующих компонентов матриц теплоемкости и теплопроводности конечных элементов [С]е, [К]е и вектора узловой нагрузки элементов (7 }е.
Для аппроксимации производной по времени в уравнении (4) применяли безусловно устойчивую конечно-разностную схему Кран-ка — Никольсона. В расчетах использовали треугольный осесимметричный симплекс-элемент, а для вычисления матриц теплоемкости, теплопроводности и вектора узловой нагрузки — формулы, приведенные в работе [4].
Моделирование структурного состава стали проводили на основании информации, снятой с изотермической диаграммы (ИТД) превращений переохлажденного аустенита. В расчете использовали ИТД стали 85ХФ [5], близкой по химическому составу к рельсовой стали 76ХФ. Эта сталь претерпевает при охлаждении два вида превращений — феррито-карбидное и мартенситное. В процессе численного моделирования на каждом шаге по времени в каждом конечном элементе вычислялся вектор удельных долей аустенита, феррито-карбида и мартенсита соответственно — {V} =^а, Vфк, Vм}. При этом плавная кривая изменения температуры в каждом конечном элементе заменялась ступенчатой, т. е. принималось, что на каждом п-м шаге по времени температура мгновенно меняется с Тп-1 на Тп и остается постоянной на данном шаге.
Для описания изотермического распада аустенита в феррито-карбид использовали уравнение Авраами [6]
Vфк(х) = 1 - ехр (- К хп), (5)
где Vфк — объемная доля феррито-карбида; К, п — зависящие от температуры эмпирические коэффициенты, определяемые по ИТД превращений переохлажденного аустенита.
Зная из ИТД стали для каждой температуры время начала Хн и конца Хк феррито-карбидного превращения, коэффициенты К и п, зависящие от температуры, можно определить по формулам [7]
где [С], [К ] — глобальные матрицы теплоемкости и теплопроводности соответственно; {Т} —
п(Т) =
2,66 ^ Хк/Х!
-; К(Т) = 0,01005хп(Т).
Согласно методу наименьших квадратов [8], изотермические диаграммы стали 85ХФ для феррито-карбидной области описывали следующими выражениями:
^ тн = 8,38 • 10-6\Г - 500|2,27 + 0,997;
lg т к = 2,12 •Ю-6 \Т - 500|2 57 + 1,82 при 370 °С<г < 680 °С; lg Тн = - 6,5435 • 10-5\Т - 200|1,87 + 2,49;
^ т к = -3,45 •Ю-5 |Т-20^2,07 + 3,82
при 220 °С < г < 370 °С.
Для перехода к неизотермической кинетике превращения применена теория изокинетиче-ских реакций [6], согласно которой объемную долю феррито-карбида на и-м шаге по времени определяют по уравнению (5) для времени ти + Дти, где ти — время, необходимое для достижения накопленной к моменту ти_1 степени превращения Уф^Т1 при температуре Ти. Тогда объемная доля феррито-карбида на и-м шаге
Уфк(Ти)=1_
( 1п(1 _ Уф";!) '
ехр
_К (Тя)
фк
К (Тя)
+Дти
Цу =Р^фк
Дтя
"; цу = р^м
Дти
случае приращение тензора полной деформации
Д£у = Де) + Дер + 8) Дет , (6)
где Дее) и Дер — приращение упругой и пластической деформации соответственно; 8) — символ Кронекера; ДеТ — приращение свободной деформации, учитывающей температурные и структурные изменения объема.
Согласно методу дополнительных деформаций, решение задачи термоупругопластичности сведено к последовательному решению задачи термоупругости. При этом два последних слагаемых в уравнении (6) объединены в одно:
ДЕ;) = Дее + Де0;,
где Де0 = Дер + 8)ДеТ — дополнительная деформация.
При определении приращений пластических деформаций принималось существование пластического потенциала, который для неизотермической теории течения в случае нестационарного структурного состава и использования критерия пластичности Хубера — Мизеса можно представить в виде [10]
¥р = (3/25)5) )1/2 _ /т (цр,Т,{У}) = 0, (7)
где Б) — девиатор напряжений; Цр — параметр Удквиста при пластичности,
Мартенситное превращение относится к атермическим превращениям, степень распада которых определяется только температурой и не зависит от времени [9], поэтому в расчете удельной доли мартенсита использована чисто температурная эмпирическая зависимость [4].
Тепловыделения при феррито-карбидном цфк и мартенситном ц" превращениях учитывали посредством включения в уравнение теплопроводности мощности удельных источников энерговыделения:
фк = , ДУифк , ДУГ
цр = | ^ ер
где !фк и Ьм — удельная теплота феррито-карбидного и мартенситного превращений соответственно; ДУифк и ДУим — изменение на и-м шаге по времени удельной доли феррито-карбида и мартенсита соответственно.
Расчет напряжений проведен посредством решения задачи термоупругопластичности для материала с нестационарной структурой [4]. В основу решения положен шаговый метод дополнительных (начальных) деформаций. В этом
(^ёр — интенсивность приращений пластических деформаций).
Вводя для упрощения обобщенный параметр О, характеризующий температурное и структурное состояние стали и учитывая, что уменьшаемое в формуле (7) представляет собой интенсивность напряжений а;, получим
а; = /т(Цр , О). (8)
Выбор условия пластичности в виде соотношения (8) равносилен гипотезе о том, что при данных температуре и структуре интенсивность напряжений является функцией параметра Удквиста, не зависящей от типа напряженного состояния. Функцию / можно получить из мгновенных кривых растяжения, представив их в виде уравнения а = /т(£р, О), в котором учтено, что для одноосного растяжения а; = а и Цр =Ер, где е р — накопленная пластическая деформация. Используя модель упругопластической среды со степенным упрочнением для описания кривых растяжения
отдельных фаз, мгновенную кривую растяжения гетерогенной структуры можно представить в виде
СТ = Бе при £<£т = (СТт.а Уа + Ът.фк Уфк + + ^т.м Ум)! Б;
т,° С
а = ат
У а + а
/ \шфк
е
т.фк
VAk +
+ ат
Шм
Ум при е >ет,
где о — напряжение; Е — модуль упругости; е — деформация; ат.а, а т.фк, а т.м и Ша, Шфк, Шм — пределы текучести и показатели упрочнения, зависящие от температуры, для аустенита, фер-рито-карбида и мартенсита соответственно;
ет.а = ат.а/Е; ет.фк = ат.фк/Е; ет.м = ат.м/Е.
Приращения пластической деформации можно рассчитать, имея мгновенную кривую растяжения и приращение интенсивности напряжений dat для каждого конечного элемента на данном шаге по времени [4]:
dep = 3/2(1/Ек " 1/E) ( dat -f dQ 1 ^, I dQ J at
где Ек = dfr/dep — касательный модуль.
По описанной методике в среде Visual Fortran создана авторская программа для расчета температур, структур и напряжений при закалке рельса. Все теплофизические коэффициенты и физико-механические характеристики приняты по работе [4]. В силу симметрии поперечного сечения рельса относительно вертикальной оси y рассматривали только правую половину сечения, которую разбивали на 1 546 треугольных симплекс-элементов. Результаты расчета показали, что охлаждение рельса до нормальной температуры происходит за 12 мин.
На рис. 2-4 приведены зависимости температуры T, объемной доли феррито-карбида Уфк и остаточных напряжений ax в головке рельса от расстояния до поверхности h. Как видно из рис. 2, в начале процесса охлаждения в масле температура рельса быстро снижается. За полторы минуты температура поверхности уменьшается до 360 °С, а внутренних областей — до 480 °С. Следует отметить, что интенсивное снижение температуры происходит и вблизи тонкой шейки, так как она быстро охлаждается.
700 600 500 400 300
V-
/ / 2 * "ч _ " ч
/ у ч \
/ / / _____ / чз ч \
0
18 27 36 h, мм
Рис. 2. Зависимость температуры в головке рельса Т от расстояния до поверхности И при различном времени охлаждения: 1 — Тохл = 44 с; 2 — Тохл = 70 с; 3 — Тохл = 90 с
Рфк,%
80 60 40 20
0
\ \ \ \ / 7" / /
\ \ \ \ 3 \ 1 / / / /
\ \ \ / / /
1 \ 2 V ч / / /
ч *----- - —
18
27 36 h, мм
Рис. 3. Формирование феррито-карбидной структуры в головке рельса при различном
времени охлаждения: 1 — Тохл = 70 с; 2 — Тохл = 90 с; 3 — Тохл = 106 с
а„ МПа
ч ч / / * ><у ч \ \ \ —___ _
J8
27.-
S36-----h,
200 0
-200 -400 -600 -800
Рис. 4. Зависимость остаточных (1) и временных (2, 3) нормальных напряжений ъх в головке рельса от расстояния до поверхности и: 2 — Тохл = 44 с; 3 — Тохл = 70 с
Из рис. 3 следует, что при времени охлаждения Тохл = 70 мин в центральных областях сохраняется структура чистого аустенита, а поверхности содержит порядка 25 % феррито-карбида. При Тохл = 90 мин на поверхности го-
Рис. 5. Распределение остаточных нормальных напряжений ох по сечению головки (а) и подошвы (б) рельса (значения указаны в МПа)
ловки и вблизи шейки структура преобразуется в чистый феррито-карбид. При Тохл = 106 мин только в центральных областях в структуре сохраняется аустенит.
Как видно из рис. 4, в начале охлаждения на поверхности головки и вблизи шейки напряжения растягивающие, а во внутренних областях сжимающие. По мере охлаждения знаки напряжений меняются и уже на поверхности напряжения становятся сжимающими, а внутри растягивающими.
Изолинии остаточных осевых напряжений ох в поперечном сечении головки и подошвы рельса приведены на рис. 5. Отсюда следует, что на внешних поверхностях головки и подошвы сжимающие напряжения достигают значений от -500 до -600 МПа. Во внутренних областях как головки, так и подошвы фор-
мируются зоны с повышенным уровнем растягивающих напряжений, достигающих 290 и 240 МПа соответственно. Проведенный расчет показал, что в шейке остаточные напряжения ох практически равны нулю, а напряжения Оу меняются по толщине от 200 МПа внутри до -370 МПа на поверхности. Общий уровень остаточных напряжений Оу ниже, чем ох. Максимальные растягивающие напряжения, расположенные во внутренних областях головки и подошвы составляют 250 и 210 МПа соответственно.
Выводы
1. При закалке в масле головка рельса имеет структуру пластинчатого феррито-карбида. Мартенсит присутствует только в структуре
шейки и вблизи пера подошвы, причем его доля достигает 65 и 95 % соответственно.
2. Уровень наиболее опасных с точки зрения хрупкой прочности растягивающих остаточных термических напряжения ох выше, чем оу. Максимальные сжимающие напряжения ох возникают на рабочей поверхности головки рельса,
Литература
а Оу — на боковой поверхности. Их значения близки и составляют 640...660 МПа.
3. Разработанные программные средства можно использовать для рационализации режимов термической обработки железнодорожных рельсов.
[1] ГОСТ P 51685-2013. Рельсы железнодорожные. Москва, Изд-во стандартов, 2001. 23 с.
[2] Цветков Ф.Ф., Григорьев Б.А. Тепломассобмен. Москва, Издательский дом МЭИ, 2006.
550 с.
[3] Zienkiewicz O.C., Taylor R.L., Fox D.D. The finite element method for solid and structural
mechanics. New York, Elsevier, 2014. 657 p.
[4] Вафин Р.К., Покровский А.М., Лешковцев В.Г. Прочность термообрабатываемых про-
катных валков. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2004. 264 с.
[5] Попов А.А., Попова Л.Е. Справочник термиста: Изотермические и термокинетиче-
ские диаграммы распада переохлажденного аустенита. Москва, Машгиз, 1961. 430 с.
[6] Christian J.W. The Theory of Transformations in Metals and Alloys. Pt. I, II Oxford, Per-
gamon Press, 2002. 1200 p.
[7] Покровский А.М., Рыжиков А.В. Математическое моделирование температурного и
фазово-структурного состояний при наплаве биметаллического прокатного валка. Машиностроение и инженерное образование, 2016, № 1, с. 42-51.
[8] Деммель Д. Вычислительная линейная алгебра: теория и приложения. Москва, Мир,
2001. 429 с.
[9] Гуляев А.П., Гуляев А.А. Металловедение. Москва, ИД Альянс, 2011. 644 с.
[10] Покровский А.М. Оценка ресурса прокатных валков с учетом остаточных напряжений от термической обработки. Производство проката, 2005, № 9, с. 26-31.
References
[1] GOST P 51685-2013. Rel'sy zheleznodorozhnye [State Standard R 51685-2013. Railway rails].
Moscow, Standartinform publ., 2001. 23 p.
[2] Tsvetkov F.F., Grigor'ev B.A. Teplomassobmen [Heat and Mass Transfer]. Moscow, MPEI
publ., 2006. 550 p.
[3] Zienkiewicz O.C., Taylor R.L., Fox D.D. The finite element method for solid and structural
mechanics. New York, Elsevier, 2014. 657 p.
[4] Vafin R.K., Pokrovskii A.M., Leshkovtsev V.G. Prochnost' termoobrabatyvaemykh pro-
katnykh valkov [The strength of heat treatable mill rolls]. Moscow, Bauman Press, 2004. 264 p.
[5] Popov A.A., Popova L.E. Spravochnik termista: Izotermicheskie i termokineticheskie dia-
grammy raspada pereokhlazhdennogo austenite [Directory treater: Isothermal and thermo-kinetic decay diagram of supercooled austenite]. Moscow, Mashgiz publ., 1961. 430 p.
[6] Christian J.W. The Theory of Transformations in Metals and Alloys. Pt. I, II Oxford, Per-
gamon Press, 2002. 1200 p.
[7] Pokrovskii A.M., Ryzhikov A.V. Matematicheskoe modelirovanie temperaturnogo i fazovo-
strukturnogo sostoianii pri naplave bimetallicheskogo prokatnogo valka [Mathematical modeling of temperature and phase-structural states during surfacing bimetallic rolling roll]. Mashinostroenie i inzhenernoe obrazovanie [Mechanical Engineering and Engineering Education]. 2016, no. 1, pp. 42-51.
[8] Demmel' D. Vychislitel'naia lineinaia algebra: teoriia i prilozheniia [Computational linear al-
gebra: theory and applications]. Moscow, Mir publ., 2001. 429 p.
[9] Guliaev A.P., Guliaev A.A. Metallovedenie [Metallography]. Moscow, ID Al'ians publ., 2011.
644 p.
[10] Pokrovskii A.M. Otsenka resursa prokatnykh valkov s uchetom ostatochnykh napriazhenii ot termicheskoi obrabotki [Evaluation of service life of rolls taking into account residual stresses due to heat treatment]. Proizvodstvo prokata [Rolled Products Manufacturing]. 2005, no. 9, pp. 26-31.
Информация об авторах
ПОКРОВСКИЙ Алексей Михайлович (Москва) — доктор технических наук, профессор кафедры «Прикладная механика». МГТУ им. Н.Э. Баумана (105005, Москва, Российская Федерация, 2-я Бауманская ул., д. 5, стр. 1, e-mail: [email protected]).
ВОРОНОВ Юрий Викторович (Москва) — аспирант кафедры «Прикладная механика». МГТУ им. Н.Э. Баумана (105005, Москва, Российская Федерация, 2-я Бауманская ул., д. 5, стр. 1, e-mail: [email protected]).
ТРЕТЬЯКОВ Денис Николаевич (Москва) — инженер кафедры «Прикладная механика». МГТУ им. Н.Э. Баумана (105005, Москва, Российская Федерация, 2-я Бауманская ул., д. 5, стр. 1, e-mail: [email protected]).
Статья поступила в редакцию 25.03.2016 Information about the authors
POKROVSKIY Aleksey Mikhailovich (Moscow) — Doctor of Science (Eng.), Professor, Department of Applied Mechanics. Bauman Moscow State Technical University (105005, Moscow, Russian Federation, 2nd Baumanskaya St., Bldg. 5, Block 1, e-mail: [email protected]).
VORONOV Yuriy Viktorovich (Moscow) — Postgraduate, Department of Applied Mechanics. Bauman Moscow State Technical University (105005, Moscow, Russian Federation, 2nd Baumanskaya St., Bldg. 5, Block 1, e-mail: [email protected]).
TRETYAKOV Denis Nikolaevich (Moscow) — Engineer, Department of Applied Mechanics. Bauman Moscow State Technical University (105005, Moscow, Russian Federation, 2nd Baumanskaya St., Bldg. 5, Block 1, e-mail: [email protected]).
В Издательстве МГТУ им. Н.Э. Баумана вышел в свет учебник Р.З. Кавтарадзе
«Локальный теплообмен в поршневых двигателях»
Учебник посвящен исследованию локального теплообмена в поршневых двигателях. Значительная его часть написана на основе результатов, полученных в МГТУ им. Н.Э. Баумана. Ряд вопросов в теории поршневых двигателей рассматривается впервые. В данный учебник, написанный на основе учебного пособия с тем же названием (1-е изд. — 2001 г., 2-е изд. — 2007 г.), включены новые материалы, отражающие достижения последних лет в этой области науки. Содержание учебника соответствует курсу лекций, который автор читает в МГТУ им. Н.Э. Баумана.
Для магистрантов, аспирантов, научных и инженерно-технических работников, занимающихся созданием перспективных двигателей, а также исследованием и доводкой существующих моделей.
По вопросам приобретения обращайтесь:
105005, Москва, 2-я Бауманская ул., д. 5, стр. 1. Теп.: +7 499 263-60-45, факс: +7 499 261-45-97; [email protected]; www.baumanpress.ru