Код:
> # Нужно вывести зависимость критерия Нуссельта при течении воды в трубе от Re и геометрических параметров турбулизатора потока типа "витая лента", таких как шаг закрутки ленты например S/d и ее толщины 1 и 2 мм
> # Уравнение должно быть вида Nu = A•Re^n •(S/d)^k•(толщина)^m
> # Вот данные: зависимости Nu=f(Re)
> #
> # Re 26575,0 22145,8 17716,6 14763,9 10334,7 7381,9
> # Nu при 1 и s/d=1 142,3 121,8 114,8 105,8 88,4 74,4
> # Nu при 1 и s/d=2 139,0 125,4 110,0 102,1 81,0 73,2
> # Nu при 2 и s/d=1 140,0 120,0 114,7 104,0 87,0 72,0
> # Nu при 2 и s/d=2 137,0 123,3 109,0 101,0 80,1 71,0
>
> Re<-rep(c(26575.0,22145.8,17716.6,14763.9,10334.7,7381.9),4)
> T<-c(rep(1,12),rep(2,12))
> Sd<-rep(c(rep(1,6),rep(2,6)),2)
> Nu<-c(142.3,121.8,114.8,105.8,88.4,74.4,
+ 139.0,125.4,110,102.1,81,73.2,
+ 140.0,120.0,114.7,104.0,87.0,72.0,
+ 137.0,123.3,109.0,101.0,80.1,71.0)
> MyData<-data.frame(Re=Re,T=T,Sd=Sd,Nu=Nu)
> MyData
Re T Sd Nu
1 26575.0 1 1 142.3
2 22145.8 1 1 121.8
3 17716.6 1 1 114.8
4 14763.9 1 1 105.8
5 10334.7 1 1 88.4
6 7381.9 1 1 74.4
7 26575.0 1 2 139.0
8 22145.8 1 2 125.4
9 17716.6 1 2 110.0
10 14763.9 1 2 102.1
11 10334.7 1 2 81.0
12 7381.9 1 2 73.2
13 26575.0 2 1 140.0
14 22145.8 2 1 120.0
15 17716.6 2 1 114.7
16 14763.9 2 1 104.0
17 10334.7 2 1 87.0
18 7381.9 2 1 72.0
19 26575.0 2 2 137.0
20 22145.8 2 2 123.3
21 17716.6 2 2 109.0
22 14763.9 2 2 101.0
23 10334.7 2 2 80.1
24 7381.9 2 2 71.0
> summary(lm1<-lm(log(Nu)~log(Re)+log(T)+log(Sd)))
Call:
lm(formula = log(Nu) ~ log(Re) + log(T) + log(Sd))
Residuals:
Min 1Q Median 3Q Max
-0.050707 -0.004867 0.009426 0.012959 0.023083
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.19233 0.10997 -1.749 0.09564 .
log(Re) 0.50427 0.01140 44.253 < 2e-16 ***
log(T) -0.02296 0.01438 -1.597 0.12605
log(Sd) -0.04180 0.01438 -2.907 0.00873 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.02442 on 20 degrees of freedom
Multiple R-squared: 0.9899, Adjusted R-squared: 0.9884
F-statistic: 656.5 on 3 and 20 DF, p-value: < 2.2e-16
> confint(lm1)
2.5 % 97.5 %
(Intercept) -0.42172064 0.037066329
log(Re) 0.48049934 0.528038653
log(T) -0.05296572 0.007039156
log(Sd) -0.07180710 -0.011802225
> summary(lm2<-lm(log(Nu)~log(Re)+log(Sd)))
Call:
lm(formula = log(Nu) ~ log(Re) + log(Sd))
Residuals:
Min 1Q Median 3Q Max
-0.057637 -0.005429 0.004838 0.016099 0.031042
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.20029 0.11384 -1.759 0.0931 .
log(Re) 0.50427 0.01181 42.706 <2e-16 ***
log(Sd) -0.04180 0.01490 -2.805 0.0106 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.0253 on 21 degrees of freedom
Multiple R-squared: 0.9887, Adjusted R-squared: 0.9876
F-statistic: 915.8 on 2 and 21 DF, p-value: < 2.2e-16
> confint(lm2)
2.5 % 97.5 %
(Intercept) -0.43702165 0.03645041
log(Re) 0.47971327 0.52882473
log(Sd) -0.07279929 -0.01081003
> oldpar<-par(mfrow=c(2,2))
> plot(lm2)
> par(oldpar)
> exp(-0.20029)
[1] 0.8184934
И диагностические графики.
Сравниваем расчетные и фактические значения
Код:
> Nu2<-0.8184934*Re^0.50427*Sd^-0.04180
> Nu2
[1] 139.36198 127.12040 113.59153 103.61395 86.55759 73.04935 135.38211 123.49012
[9] 110.34760 100.65496 84.08569 70.96322 139.36198 127.12040 113.59153 103.61395
[17] 86.55759 73.04935 135.38211 123.49012 110.34760 100.65496 84.08569 70.96322
> Nu2-Nu
[1] -2.93801638 5.32040349 -1.20846824 -2.18605295 -1.84240704 -1.35065290 -3.61789408
[8] -1.90988108 0.34760260 -1.44504382 3.08569478 -2.23678471 -0.63801638 7.12040349
[15] -1.10846824 -0.38605295 -0.44240704 1.04934710 -1.61789408 0.19011892 1.34760260
[22] -0.34504382 3.98569478 -0.03678471
Примерно так...
Если же делать с помощью нелинейного МНК
Код:
> summary(nlm1<-nls(Nu~A*Re^n*Sd^k*T^m,start=list(A=1,n=1,k=1,m=0)))
Formula: Nu ~ A * Re^n * Sd^k * T^m
Parameters:
Estimate Std. Error t value Pr(>|t|)
A 0.79250 0.10084 7.859 1.53e-07 ***
n 0.50806 0.01296 39.198 < 2e-16 ***
k -0.03333 0.01434 -2.325 0.0307 *
m -0.02086 0.01434 -1.455 0.1611
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.631 on 20 degrees of freedom
Number of iterations to convergence: 16
Achieved convergence tolerance: 3.872e-07
> summary(nlm2<-nls(Nu~A*Re^n*Sd^k,start=list(A=1,n=1,k=1)))
Formula: Nu ~ A * Re^n * Sd^k
Parameters:
Estimate Std. Error t value Pr(>|t|)
A 0.78664 0.10265 7.664 1.63e-07 ***
n 0.50808 0.01330 38.195 < 2e-16 ***
k -0.03333 0.01471 -2.265 0.0342 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.7 on 21 degrees of freedom
Number of iterations to convergence: 16
Achieved convergence tolerance: 1.913e-07
Получаем несколько иные коэффициенты в пределах доверительного интервала. В обоих случаях толщина у нас не при делах.
Код:
> (Nu3<-0.78664*Re^0.50808*Sd^-0.03333)
[1] 139.23948 126.92047 113.31649 103.29129 86.17087 72.62981 136.05956 124.02189
[9] 110.72860 100.93235 84.20292 70.97111 139.23948 126.92047 113.31649 103.29129
[17] 86.17087 72.62981 136.05956 124.02189 110.72860 100.93235 84.20292 70.97111
> Nu3-Nu
[1] -3.06051548 5.12046869 -1.48350586 -2.50870593 -2.22912754 -1.77019073 -2.94043534
[8] -1.37811227 0.72859768 -1.16764914 3.20292157 -2.22889389 -0.76051548 6.92046869
[15] -1.38350586 -0.70870593 -0.82912754 0.62980927 -0.94043534 0.72188773 1.72859768
[22] -0.06764914 4.10292157 -0.02889389
> confint(nlm2)
2.5% 97.5%
A 0.59839995 1.031156973
n 0.48048555 0.535954713
k -0.06394209 -0.002731841
Относительная ошибка аппроксимации составляет для первого и второго решения соответственно (в процентах):
Код:
> mean(abs(Nu2-Nu)/Nu)*100
[1] 1.825308
> mean(abs(Nu3-Nu)/Nu)*100
[1] 1.877169