我必須解決 C 中的有限差分問題（2d 傳熱，非瞬態），但是當我嘗試為 (nx xm) 陣列解決此問題時，對于 nx >70，程式崩潰并回傳 -1073741819 (0xC0000005) . 我該怎么辦？最后，腳本應該計算復合介質中的溫度場。這段代碼試圖證明表面 T1[m-1][:]（它的最后一個元素）的溫度。
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
int main(){
int i, j,k;
int nx = 200; // número de pontos em x
int n = 50;
int m=n; // número de pontos em b1 e b2
srand((unsigned)time(NULL));
// Parametros de simula??o
float b1 = 0.01; // altura do meio 1, m
float b2 = 0.01; // altura do meio 2, m
float a = 0.04; // largura das placas, m
float q = -100000.0; // fonte de calor na face superior [W/m2]
float k1e = 54.0; // condutividade térmica do material 1 [W/mC]
float k2 = 54.0; // condutividade térmica do material 1 [W/mC]
float Tb = 0; //temperatura de referência .C
float hmax = 1000;
float x[nx], y1[m], y2[n];
float dx = a/(nx-1);
float dy = b1/(n-1);
y1[0] = b1; x[0] = 0.0; y2[0] = 0.0;
for (i=1; i<nx; i ){
x[i] = x[i-1] dx;
}
for (j=1; j<m; j ){
y1[j]= y1[j-1] dy;
y2[j]= y2[j-1] dy;
}
float T1[n][nx], T2[m][nx];
float T1old[n][nx], T2old[m][nx];
//=============== Cálculo de hc =================
// Caso 1
//char case = 'CP1h1 '
float hc[nx];
for (i=0; i<nx;i ){
if (x[i]<a/4)
hc[i] = hmax;
else if (x[i] > 3*a/4)
hc[i] = hmax;
else
hc[i] = 0.;
}
//==============================================
int nexp = 1;
float diff = 0.0, maxdiff = 0.0;
float err1 = 10.0;
float err2 = 10.0;
float tol = pow(10.0,-12);
int contador = 0;
float k1=k1e;
for (int k=0; k<nexp; k ){
//=======================================================================================
// Valor sintético de k1
float sigmak1 = 0.00;
float sigma = sigmak1*k1e;
float u = (float)rand()/ RAND_MAX;
float v = (float)rand()/ RAND_MAX;
float eps = (2*M_PI*v)*sqrt(-2*log(u));
float k1 = k1e eps*sigma;
//======================================================================================
float k1dy=k1/dy;
float k1k2=k1/k2;
float twoqdyk1=2*q*dy/k1;
for (i=0; i<nx; i ){
for (j=0; j<m; j ){
T1[j][i] = 0.;
T2[j][i] = 0.;
T1old[j][i] = 0.;
T2old[j][i] = 0.;
}
}
diff = abs(T1old[0][0] - T1[0][0]);
maxdiff = diff;
contador = 0;
while (contador<14000){
contador = contador 1;
for (i=0; i<nx; i ){
for (j=0; j<nx;j ){
diff = abs(T1old[j][i] - T1[j][i]);
T1old[j][i] = T1[j][i];
T2old[j][i] = T2[j][i];
printf("%f \n",diff);
if (diff>=maxdiff){
maxdiff = diff;
}
}
}
err1 = maxdiff;
//printf("%i \n", contador);
for (i=0; i<nx;i ){
for(j=0; j<m;j ){
if (i==0 && j!=0 && j!=m-1){ //condi??o de contorno direita -> exclui j=0 e j=n-1
T1[j][i] = 0.25*(T1old[j-i][i] T1old[j i][i] 2*T1old[i 1][j]);
T2[j][i] = 0.25*(T2old[j-i][i] T2old[j i][i] 2*T2old[i 1][j]);
}
else if (i==nx-1 & j!=0 && j!=m-1){ //condi??o de contorno direita -> exclui j=0 e j=m-1
T1[j][i] = 0.25*(T1old[j-i][i] T1old[j i][i] 2*T1old[i-1][j]);
T2[j][i] = 0.25*(T2old[j-i][i] T2old[j i][i] 2*T2old[i-1][j]);
}
else if(j==0){ //condi??o de contorno inferior -> inclui i=0 e i=nx-1
T1[j][i] = (1/(k1dy hc[i]))*(hc[i]*T2old[n-1][i] (k1dy)*(T1old[j 1][i]));
T2[j][i] = Tb;
}
else if(j==m-1){//condi??o de contorno superior -> inclui i=0 e i=nx-1
T1[j][i] = 0.25*(2*T1old[j-1][i] -(2*twoqdyk1) T1old[j][i-1] T1old[j][i-1]);
T2[j][i] = (-k1k2)*(T1old[0][i]-T1[1][i]) T2old[j-1][i];
}
else{ //pontos internos
T1[j][i] = 0.25*(T1old[j-1][i] T1old[j 1][i] T1old[j][i-1] T1old[j][i-1]);
T2[j][i] = 0.25*(T2old[j-1][i] T2old[j 1][i] T2old[j][i-1] T2old[j][i-1]);
}
}
}
}
}
for (i=1; i<nx;i ){
printf("%f \n ", T1[m-1][i]);
}
return 0;
}
uj5u.com熱心網友回復：
至少這些問題
并非所有警告都已啟用
節省時間，啟用它們。
int 絕對
abs(T1old[0][0] - T1[0][0]);確定int絕對值。這會導致不正確的功能和未定義的行為，其值遠遠超出int范圍。
使用浮點函式： fabsf(T1old[0][0] - T1[0][0]);
log(0.0} 可能是 -INF
下面可能會產生 u == 0.0
float u = (float)rand()/ RAND_MAX;
float eps = (2*M_PI*v)*sqrt(-2*log(u)); // Bad
也許 ？
float u = (rand() 0.5f)/ (RAND_MAX 1u);
為什么float？
代碼將double函式呼叫和常量與float變數混合在一起。建議只使用double整個。
除錯提示
注釋掉，srand((unsigned)time(NULL));直到代碼完全正確運行至少一次。
至少在除錯期間使用%g或%e而不是列印浮點值。%f它資訊量更大，噪音更小。
uj5u.com熱心網友回復：
仔細看這段代碼：
for (i = 0; i < nx; i ) {
for (j = 0; j < nx; j ) {
diff = abs(T1old[j][i] - T1[j][i]);
內回圈j從 0 變化到nx但應該是n.
for (i = 0; i < nx; i ) {
for (j = 0; j < n; j ) { // << use n instead of nx
diff = abs(T1old[j][i] - T1[j][i]);
的定義T1：
float T1[n][nx]
Sou 您在這里所擁有的是訪問具有超出范圍索引的陣列，這會導致未定義的行為。
轉載請註明出處，本文鏈接：https://www.uj5u.com/shujuku/413801.html
標籤：