книги / Теория инженерного эксперимента
..pdf11.7. Моделирование на ЭВМ
Современная вычислительная машина становится од ним из средств, используемых при проведении эксперимен тов. Сложные технические системы, не поддающиеся даже самым разработанным формам математического анализа, часто можно легко промоделировать с помощью последо вательности простых взаимосвязанных уравнений, пол ностью описывающих схему поведения. Затем эти систе мы уравнений решаются при различных значениях вход ных параметров, что позволяет исследовать область зна чений зависимой переменной и соотношение между вход ными и выходными величинами.
Рассмотрим простой пример. Геологические исследо вания показали, что бассейны рек приобретают ту или иную форму в результате огромного множества неболь ших случайных воздействий за огромный промежуток времени. В совокупности множество этих случайных со бытий привело к образованию бассейнов, имеющих раз личное число притоков разной длины, ширины и формы. Такое разнообразие характерно для бассейнов всех рек земного шара. Проверить эту теорию экспериментальным путем невозможно, но возможна ее проверка путем моде лирования на цифровой вычислительной машине с по мощью так называемого метода случайного блуждания.
Результаты моделирования обнаруживают поразитель ное совпадение с теорией. При моделировании бассейны рек были описаны количественно, при этом длина, раз меры и распределения вероятностей почти полностью совпадали с реально наблюдаемыми в природе [5].
Трудно представить себе все области возможного при менения моделирования в геологии и биологии, задачах движения транспорта, теплотехнике, механике и общест венных науках. Моделирование требует применения совер шенно других (более простых) математических методов, которые в настоящее время играют важную роль в инже нерной деятельности. Наших студентов уже можно учить тому, ^как решать с помощью вычислительных машин такие задачи, с которыми не в состоянии справиться даже доктора наук, если они будут применять формальные мате матические методы. Применение моделирования встречает
значительное сопротивление со стороны тех, кто привык к старым методам и затратил в прошлом много труда на их изучение и освоение.
Моделирование дает возможность исследователю пол ностью провести эксперимент на вычислительной машине. Хотя в данном случае не существует такой проблемы, как точность прибора, ошибка появляется в виде неоп ределенностей физических свойств и констант, исполь зуемых в программе для моделирования. Моделирование и параметрические исследования могут дополняться ана лизом размерностей и проведением факторного экспери мента (гл. 6), позволяющего значительно сократить число вариантов, исследуемых с помощью вычислительной ма шины.
По-видимому, наиболее подходящим случаем приме нения моделирования в экспериментальной работе яв ляется оценка ошибки результата методами, рассмотрен ными в гл. 3, когда программа эксперимента является очень сложной. Если работу всей измерительной системы удается промоделировать на ЭВМ, то можно ввести раз личные ошибки и, тщательно проверив их комбинации, построить полную область ошибок, позволяющую оценить наихудшие условия, ожидаемые при проведении экспери мента. Однако предстоит еще много поработать над ана лизом такого рода.
Моделирование завершает круг вопросов, рассматри ваемых в данной книге. В некоторых отношениях модели рование ближе к эксперименту, чем к анализу. Можно ожидать, что в будущем оно будет играть важную роль при проведении самых сложных экспериментов.
ЗАДАЧИ
11.1 Составьте программу на языке ФОРТРАН для ввода данных из табл. 2.1 и записи интерполированных значений, вычисляемых по формуле Лагранжа.
11.2. Составьте программу на языке ФОРТРАН, обе спечивающую выбор N значений X и Y и построение по значениям X и Y (ошибку содержит переменная Y) пря мой наименьших квадратов с помощью методов, рассмо тренных в гл. 9.
11.3.Составьте программу на языке ФОРТРАН, обе* спечивающую подбор квадратного уравнения для данных из примера 1 0 .1 .
11.4.Доставьте программу на языке ФОРТРАН, что бы проверить, будут ли N точек с координатами X и Y лежать на прямой при построении графика в линейных, полулогарифмических и логарифмических координатах. Выберите количественный критерий оценки «прямоли нейности», обеспечивающий выполнение вычислений на ЭВМ.
11.5.Составьте программу, обеспечивающую: а) вы числение параметров пуассоновского распределения, опи санного в гл. 8 , для N точек; б) проверку соответствия данных пуассоновскому распределению с помощью кри терия х2 (без записи значений х2 и вероятностей). Составь те логическую схему таким образом, чтобы число точек
вкаждой ячейке было не меньше пяти.
1 1.6 . Составьте программу для анализа латинского квадрата 3 x 3 , аналогичного анализу, выполненному в примере 6 .2 .
ЛИ ТЕРАТУРА
1. H a m m i n g R. W., Numerical Methods for Scientists and Engi neers, McGraw-Hill, 1962.
2. M c C r a c k e n D. D., A Guide to FORTRAN Programming. Wiley, N.Y., 1961.
3.R a 1 s t о n A., W i 1 f |W ., Mathematical Methods for Digital Computers, Wiley, N.Y., 1960.
4. |
S c h e n c k |
H., |
Jr., |
FORTRAN Methods in Heat Flow, |
Ronald |
5. |
Press, N.Y., 1963. |
Jr., Simulation of the Evolution of Drainage- |
|||
S c h e n c k |
H., |
||||
|
basin Networks with a Digital Computer, J . Geophys. Res., 68, 5739— |
||||
6. |
5745 (October 15, |
1963). |
|
||
S h e r m a n |
P. M., |
Programming and Coding Digital Computers, |
|||
7. |
Wiley, N.Y., 1963. |
|
N.Y., |
||
W e i s s E., |
Programming the IBM 1620, McGraw-Hill, |
||||
|
1965. |
|
|
|
|
23— 168
Приложение А
Три простых эксперимента
Хотя невозможно дать даже беглый обзор множества экспериментов, проводимых в высших учебных заведе ниях, при переподготовке инженеров и в промышлен ности, все же имеет смысл детально изучить три простых эксперимента с точки зрения их планирования, методики проведения и анализа данных. Наша задача состоит здесь не в том, чтобы выработать определенную последователь ность действий или форму для составления отчета, а в том, чтобы показать применение материала различных разделов данной книги при проведении экспериментов, известных большинству инженеров.
Пример А.1. Измерение вязкости с помощью рота ционного вискозиметра Стормера. Вискозиметр Стормера схематически изображен на фиг. А.1. При обозначе ниях, принятых на этой схеме, вязкость находится из формулы1
_ |
д е I |
(Ь2 — a 2) k 2g |
(А.1 |
) |
||
" — |
1 / Д т |
4 n a 2b2S ( / + |
е) |
|||
|
|
|||||
где I — толщина слоя жидкости, а |
е — дополнительная |
|||||
фиктивная толщина, обусловленная трением о дно вра щающегося цилиндра.
Величины a, b, k, S и I измеряются с помощью штан генциркуля или линейки. Будем считать, что максималь
но |
возможная неопределенность составляет |
0,076 см |
|||
для |
всех величин. |
Значение е |
принимается |
равным |
|
6 , |
1 см (предполагается, что эта величина известна точно, |
||||
хотя, |
по-видимому, |
некоторая |
неопределенность все же |
||
1 V e n n a r d J. К.-, Elementary Fluid Mechanics, Wiley, 1949, pp. 227—229.
имеется). Отношение Д0/(1/Дт) будет определено далее в процессе эксперимента. Записав формулу (А.1) в виде
__ еде
^ 1/Д т *
(А.2 )
найдем ошибку результата С и примем решение относи тельно целесообразности повторного измерения размер-
Ф и г. А .1 . С хем а вискозим етра |
С торм ера. |
И сп ользую тся о б о зн а |
чения, приняты е |
в ф орм ул е |
(А .1 ). |
ных величин с целью уменьшения неопределенности. Со ставим таблицу величин, необходимых для определе ния С:
Величина
|
а |
\ |
Ь |
\ k |
5 |
1 |
е |
1+ е |
Значение, |
5,00 |
|
6,05 |
1,55 |
|
|
|
13,72 |
см |
|
57,15 |
7,60 |
6,10 |
||||
Н еопреде |
± 0 ,0 7 6 |
|
|
|
|
|
|
|
ленность, см |
|
± 0 ,0 7 6 |
± 0 ,0 7 6 |
± 0 ,0 7 6 |
± 0 ,0 7 6 |
0 |
± 0 ,0 7 6 |
23*
Будем |
считать, что ускорение |
силы тяжести |
g ~ |
|
= 981 |
см/сек2 — точная |
величина. |
|
|
Чтобы упростить анализ ошибки результата, за |
||||
пишем выражение для С в следующем виде: |
|
|||
|
с = ° |
( • ? — |
? - ) • |
<А - 3 > |
где D — k2g/4nS(l + е). После подстановки числовых значений этих величин получаем D = 0,237. Теперь при меним формулу (3.15) к каждому из неопределенных членов формулы (А.З). Мы видим, что
дС ___ |
2D |
_ |
« 0,237 |
0,00379. |
|
да |
а » |
~ ~ |
Z '(5 ,0 )3 |
||
|
Тогда (dC/dafwl = (—0,00379)2(0,076)2 = 8,33* 10“8. Это значение характеризует неопределенность, которую вно сит величина а в С2.
Аналогично
дС _ |
2 D |
2-0,237 |
0,00214 |
|
db |
‘ Ь3 |
(6,05)» |
||
|
и(dC/mfwl = (0,00214)2 • (0,076)2 = 2,66-10"®. Исследуем теперь неопределенность, вносимую k. Из
формулы |
(А.1) |
|
|
|
|
|
дС |
b2 —а2 |
[ 4 n S « + |
e) ] 2k. |
|
|
dk |
а2Ь2 |
|
||
Член (6 |
2 — а2 )/а2 6 |
2 равен |
0,0127. |
После |
подстановки |
значений g, S, (I + е) и k |
получаем дС/dk |
= 0,00388 и |
|||
(" S ')* w*= (°>0°388)2 • (0,076)2= 8,43 -10"8.
Таким образом, неопределенность, вносимая k, примерно такого же порядка, что и неопределенность, вносимая а.
Далее найдем неопределенность, вносимую S. Снача ла вычислим производную
дС |
Ь2 — а2 |
k2g |
w 1 |
dS |
a2b2 |
4 n ( l + e) |
S* |
Подставив соответствующие числовые значения, найдем, что
Влияние S пренебрежимо мало, и можно показать, что то же самое справедливо и в отношении (/ + е). Таким обра зом, соотношение (3.15) можно записать в виде
откуда |
wc =(8,33.10-® +2,66.10-® |
+ 8,43* 10- 8 ) 1 / 2= |
= 4,22-10-4; С = 3 ,0 Ы 0 _3, поэтому |
неопределенность |
|
результата, обусловленная выполнением измерений, со ставляет (4,22-10г4 )/(3,0Ы0-®), или 14%. [Уменьшения этой довольно большой неопределенности можно достиг нуть двумя способами. Можно применить более точные измерительные приборы и методы измерений и тем самым уменьшить неопределенность, составляющую 0,076 см, либо многократно измерять a, b a k, чтобы повысить об щую точность измерений. Для снижения неопределенности до 5% потребуется девятикратное измерение этихХтрех важных величин, разумеется, при условии, что неопре деленность полностью обусловлена наличием случайных
ошибок, а |
не вызывается систематической ошибкой. |
! Теперь |
необходимо найти Д0/(1/Дт) — тангенс угла |
наклона кривой зависимости времени падения груза от величины, обратной массе, при фиксированных значениях всех остальных величин (включая S и /). Теория предска зывает, что график будет иметь вид прямой, а отсюда сле дует, что анализ должен выполняться методом наимень ших квадратов. К сожалению, в распоряжении исследова теля нет достаточного количества дробных значений мас сы, чтобы значения /и- 1 изменялись непрерывно и чтобы можно было применить графический 'метод Асковица (см. разд. 9.2). Вязкость очень сильно зависит от темпе ратуры. Поэтому было бы.весьма неразумно изменять массу в возрастающем или убывающем порядке, посколь ку в этом случае будет также несколько увеличиваться или уменьшаться температура масла. С помощью какого-
либо игрового метода план эксперимента легко рандоми зировать (разд. 6 .2 ).
Типичный набор данных (записанных не в порядке
выбора) |
имеет вид |
|
|
|
|
|
Масса т, |
г |
10 |
20 |
30 |
40 |
50 |
Время падения при |
2,8 |
1,7 |
1,2 |
1,0 |
0,7 |
|
S = 5 7 ,15 см, сек |
||||||
(Масса)'1, г"1 |
0,100 |
0,050 |
0,033 |
0,025 |
0,020 |
|
Хотя совершенно одинаковые интервалы между значения ми тГ1 невозможны, все же их, конечно, можно выб рать более удачно при условии, что точность измерения времени одинакова на всей кривой (при использовании уравнений наименьших квадратов это условие должно выполняться). Если же окажется, что неопределенность при измерении интервалов времени больше в области ма лых интервалов, то такой выбор значений массы представ ляется удачным, так как на графике зависимости времени
от тг1 большее число |
точек попадает на участок малых |
||
интервалов. |
прямую, |
удовлетворяющую следующим |
|
|Построим |
|||
данным: |
|
|
|
m - i, X |
е, у |
т - 1 0 , X Y |
(т-1) 2, Х2 |
0,1 |
2,8 |
0,28 |
0,01 |
0,05 |
1,7 |
0,085 |
0,0025 |
0,033 |
1,2 |
0,0395 |
0,00109 |
0,025 |
1,0 |
0,025 |
0,000625 |
0,020 |
0,7 |
0,014 |
0,000400 |
1тг*е=0,228 |
10= 7,4 |
1т->0«О,4435 |
2 (т _1)2= 0 ,014615 |
Здесь использованы обозначения, принятые в форму
лах (9.6) и (9.7). |
находим V — величину |
||
По формулам (9.6) и (9.7) |
|||
отрезка, отсекаемого прямой |
на |
координатной оси 0 , |
|
и а' — угловой коэффициент |
прямой: |
||
и, |
0,0146-7,4 — 0,228-0,4435 _ _ ЛООО |
||
0 — |
5-0,0146 — 0,228® |
— v .ooo, |
|
|
5-0,4435 - 0,228-7,4 _ |
ПЛ с |
|
“ — |
5-0,0146 — 0,228* |
|
|
Соответствующие точки и сама прямая показаны на фиг. А.2.
Теперь можно рассмотреть точки, отклоняющиеся от прямой. Явно сомнительная точка находится в области наименьших значений времени, и ее нельзя исключать до тех пор, пока в области малых значений 0 не будет по лучено большее число точек. В таких простых экспери ментах легче снять несколько повторных отсчетов, чем отбрасывать сомнительные точки. Заметим, что прямая
Ф и г . А.2. Прямая, построенная методом наименьших квадратов, для зависимости времени 0 паде ния груза с высоты S от 1/т.
на фиг. А.2 пересекает ось 0 в точке 0,333 сек, поэтому возможна проверка с помощью экстраполяции. Вопрос заключается в том, какой должна быть величина отрезка, отсекаемого прямой на оси 0, и каким образом можно предсказать его величину. Изучение схемы на фиг. А.1 показывает, что при п г 1 -*■ 0 (при бесконечном увеличе нии т) тормозящий эффект трения жидкости о стенки вискозиметра становится пренебрежимо малым и время падения определяется с помощью обычного уравнения для движения с ускорением под действием силы тяжести
S = - Ç - + V 0Q, |
(А.4) |
где S =57,15 см; g =981 см/сек*; V0 — начальная ско' рость движения в момент попадания на хронометрируе мый участок.
Расстояние от начала движения до начала хронометри руемого участка выбирается таким образом, чтобы груз любой заданной массы достиг постоянной или оконча тельной скорости до момента входа на регистрируемый участок. Допустим, что некоторый груз при падении про ходит расстояние 57,15 см, а на следующем участке дли ной 57,15 см проводится хронометрирование (фактически начальный участок составил около 28,5 см). На преодоле ние первого отрезка длиной 57,15 см требуется (2S/g)lf* = = 0,342 сек. Для прохождения всего расстояния, равно го 104,3 см, требуется 0,482 сек. Таким образом, для про хождения второго отрезка длиной 57,15 см требуется 0,482 — 0,342 = 0,14 сек. Так как груз не проходит полностью начальный отрезок длиной 57,15 см, истинное время падения на хронометрируемом участке длиной 57,15 см заключено в интервале от 0,342 сек (время дви жения из состояния покоя) до 0,14 сек (время движения при старте с отметки 57,15 см). Поэтому полученное путем экстраполяции значение, равное 0,333 сек, правдо подобно, хотя, возможно, является несколько завышен ным.
Теперь следует оценить неопределенность или ошибку в определении углового коэффициента прямой, так как для определения вязкости по формуле (А.1 ) необходимо знать величину углового коэффициента. Угловой коэффи циент является функцией Д0 и 1/Дт, но масса т известна точно, так как используются стандартные значения мас сы. В конечном счете А0 = 2,8 сек и абсолютная макси мальная ошибка в измерении времени при тГ1 = 0 со ставляет 0,333—0,140 = 0,193 сек (хотя, по-видимому, она в два раза меньше). Таким образом, максимально воз можная ошибка составляет 0,193/2,8, т. е. около 7%; наиболее вероятно, что эта ошибка равна 3,5%.
Наконец, по формуле (А.2 ) можно определить вяз кость:
р=24,5-(3,01 • 10'3 )=7,374-10' 2 е/см -сек.
С учетом неопределенности измерений в отчете об испыта ниях следует указать результат с точностью до двух зна чащих цифр: 0,074 г/см• сек. Фактически при максималь ной ошибке в определении С порядка 14% и ошибке в