Показать сообщение отдельно
Старый 04.03.2015, 21:39   #19
Hogfather
Platinum Member
 
Аватар для Hogfather
 
Регистрация: 22.07.2010
Адрес: Санкт-Петербург
Сообщений: 3,281
По умолчанию

Код:
> # Нужно вывести зависимость критерия Нуссельта при течении воды в трубе от 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
---------
DNF is not an option
Hogfather вне форума   Ответить с цитированием
Реклама