1
2 #include <cgx.h>
3 extern Values *value;
4
5
6
write2duns(char * datout,Summen * anz,Nodes * node,Elements * elem,Datasets * lcase,NodeBlocks * nBlock,int bouNr,NodeBlockbou * dunsbou)7 int write2duns( char *datout, Summen *anz, Nodes *node, Elements *elem, Datasets *lcase, NodeBlocks *nBlock, int bouNr, NodeBlockbou *dunsbou )
8 {
9 FILE *handle1, *handle2, *handle3;
10 int i,j,k,n,b=0,edges=0;
11 double dl, p1=-1, pt1=-1,p2=-1, vu,vv,vel,flowangle;
12 char buffer[MAX_LINE_LENGTH];
13 double gasConstant=287.1, gamma, dynVisc, dynVisc0=1.786e-5, temp,tt1, cpGas, wl, pr, massflow, velFactor=1.;
14
15 /* Open the files and check to see that it was opened correctly */
16 i=strlen(datout);
17 strcpy (datout, "duns.grid");
18 handle1 = fopen (datout, "w+b");
19 if (handle1==NULL) { printf ("\nThe output file \"%s\" could not be opened.\n\n",
20 datout); return(-1);}
21 else printf (" file %s opened\n",datout);
22
23 strcpy (datout, "duns.conn");
24 handle2 = fopen (datout, "w+b");
25 if (handle2==NULL) { printf ("\nThe output file \"%s\" could not be opened.\n\n",
26 datout); return(-1);}
27 else printf (" file %s opened\n",datout);
28
29 strcpy (datout, "duns.bou");
30 handle3 = fopen (datout, "w+b");
31 if (handle3==NULL) { printf ("\nThe output file \"%s\" could not be opened.\n\n",
32 datout); return(-1);}
33 else printf (" file %s opened\n",datout);
34
35 printf (" write duns mesh \n");
36
37 /* grid file */
38 fprintf (handle1, "%d\n", anz->b);
39 for (b=0; b<anz->b; b++)
40 {
41 if(nBlock[b].dim==2) fprintf (handle1, "%d %d\n", nBlock[b].i, nBlock[b].j);
42 else fprintf (handle1, "%d %d %d\n", nBlock[b].i, nBlock[b].j, nBlock[b].k);
43 n=0;
44 for (k=0; k<nBlock[b].k; k++)
45 {
46 for (j=0; j<nBlock[b].j; j++)
47 {
48 for (i=0; i<nBlock[b].i; i++)
49 {
50 if(nBlock[b].dim==2) fprintf (handle1, "%le %le\n", node[nBlock[b].nod[n]].nx, node[nBlock[b].nod[n]].ny);
51 else fprintf (handle1, "%le %le %le\n", node[nBlock[b].nod[n]].nx, node[nBlock[b].nod[n]].ny, node[nBlock[b].nod[n]].nz);
52 n++;
53 }
54 }
55 }
56 }
57 fclose(handle1);
58
59 /* connectivity file */
60 fprintf (handle2, "%3d blocks\n", anz->b);
61 for (b=0; b<anz->b; b++)
62 {
63 fprintf (handle2,"B %3d",nBlock[b].geo+1-nBlock[0].geo );
64 if(nBlock[b].dim==2) edges=4;
65 if(nBlock[b].dim==3) edges=6;
66 for(i=0; i<edges; i++)
67 {
68 if(nBlock[b].map[i][0]>-1) fprintf(handle2," b 0%3d %1d%1d%1d", nBlock[b].neighbor[i]+1-nBlock[0].geo, nBlock[b].map[i][0], nBlock[b].map[i][1], nBlock[b].map[i][2]);
69 else fprintf(handle2," s %2d 0 000", nBlock[b].neighbor[i]);
70 }
71 fprintf (handle2,"\n" );
72 }
73 fclose(handle2);
74
75 // write the gas-properties and initialization block only if the temperature is defined
76 sprintf(buffer,"temperature");
77 if(getValuNr(buffer)>-1)
78 {
79 /* boundary file */
80
81 temp=atof(value[getValuNr(buffer)].string);
82 cpGas=1005.+0.0785*(temp-273.15); // sp. heat cap. J/Kg/K
83 gamma=1./(1-(gasConstant/cpGas));
84 wl=0.0245*(1.+0.00225*(temp-273.15));
85 // sutherland:
86 dynVisc=pow((temp/273.),(3./2.)) * (273.+111.)/(temp+111.) * dynVisc0;
87 pr=dynVisc*cpGas/wl;
88 fprintf(handle3," gas-properties\n");
89 fprintf(handle3," compressible\n");
90 fprintf(handle3," gas-constant %f\n",gasConstant);
91 fprintf(handle3," gamma %f\n",gamma);
92 fprintf(handle3," viscosity %e\n",dynVisc);
93 fprintf(handle3," prandtl-number %f\n",pr);
94 fprintf(handle3," turbulence\n");
95 fprintf(handle3," prandtl-number %f\n",pr*1.33);
96
97 /* initialize */
98
99 fprintf(handle3," initialize\n");
100 sprintf(buffer,"ini-velocity-factor");
101 if(getValuNr(buffer)>-1)
102 {
103 velFactor=atof(value[getValuNr(buffer)].string);
104 }
105 sprintf(buffer,"pressure");
106 if(getValuNr(buffer)>-1)
107 {
108 p2=atof(value[getValuNr(buffer)].string);
109 sprintf(buffer,"stagnation-pressure");
110 if(getValuNr(buffer)>-1)
111 {
112 pt1=atof(value[getValuNr(buffer)].string);
113 fprintf(handle3," pressure %f\n", (p2+pt1)*.5);
114 }
115 else fprintf(handle3," pressure %f\n", p2);
116 }
117 fprintf(handle3," temperature %f\n",temp);
118 sprintf(buffer,"u-velocity");
119 if(getValuNr(buffer)>-1)
120 {
121 vu=atof(value[getValuNr(buffer)].string);
122 fprintf(handle3," u-velocity %f\n",vu*velFactor);
123 }
124 sprintf(buffer,"v-velocity");
125 if(getValuNr(buffer)>-1)
126 {
127 vv=atof(value[getValuNr(buffer)].string);
128 fprintf(handle3,"# v-velocity %f\n",vv*velFactor);
129 }
130 fprintf(handle3," v-velocity 0.\n");
131 fprintf(handle3," turbulent-intensity 0.06\n");
132 fprintf(handle3," turbulent-viscosity %e\n", dynVisc*10.);
133 fprintf(handle3,"\n");
134 }
135 vel=sqrt(vu*vu+vv*vv);
136 flowangle=asin(-vv/vel);
137
138 fprintf(handle3," boundary-conditions \n");
139 for (b=0; b<bouNr; b++)
140 {
141 fprintf(handle3,"#\n# set:%s\n", dunsbou[b].name);
142 fprintf(handle3," surface ");
143 for(i=0; i<dunsbou[b].surfs; i++) fprintf(handle3,"%d ", dunsbou[b].surf[i]);
144 fprintf(handle3,"\n %s\n", dunsbou[b].bctype);
145 if(compare(dunsbou[b].bctype,"inviscid-wall",11)==11)
146 {
147 sprintf(buffer,"WallTemperature");
148 if(getValuNr(buffer)>-1)
149 { fprintf(handle3," temperature %s\n",value[getValuNr(buffer)].string); }
150 else fprintf(handle3," adiabatic\n");
151 }
152 if(compare(dunsbou[b].bctype,"viscous-wall",11)==11)
153 {
154 sprintf(buffer,"WallTemperature");
155 if(getValuNr(buffer)>-1)
156 { fprintf(handle3," temperature %s\n",value[getValuNr(buffer)].string); }
157 else fprintf(handle3," adiabatic\n");
158 }
159 if(compare(dunsbou[b].bctype,"subsonic-outflow",14)==14)
160 {
161 sprintf(buffer,"pressure");
162 if(getValuNr(buffer)>-1) fprintf(handle3," pressure %s\n",value[getValuNr(buffer)].string);
163 }
164 if(compare(dunsbou[b].bctype,"subsonic-inflow",14)==14)
165 {
166 sprintf(buffer,"stagnation-temperature");
167 if(getValuNr(buffer)>-1)
168 {
169 tt1=atof(value[getValuNr(buffer)].string);
170 fprintf(handle3," stagnation-temperature %e\n",tt1);
171 sprintf(buffer,"stagnation-pressure");
172 if(getValuNr(buffer)>-1)
173 {
174 pt1=atof(value[getValuNr(buffer)].string);
175 fprintf(handle3," stagnation-pressure %e\n",pt1);
176 }
177 fprintf(handle3,"# velocity %f\n",vel);
178
179 p1=pt1*pow( (temp/tt1), (gamma/(gamma-1)) );
180 // specific massflow per area
181 massflow= vu * p1/gasConstant/temp;
182 fprintf(handle3,"# mass-flow calculation based on p:%e\n",p1);
183 fprintf(handle3,"# mass-flow %e\n",massflow);
184 fprintf(handle3,"# temperature %e\n",temp);
185 }
186 else
187 {
188 fprintf(handle3," velocity %f\n",vel);
189 sprintf(buffer,"temperature");
190 if(getValuNr(buffer)>-1) fprintf(handle3," temperature %s\n",value[getValuNr(buffer)].string);
191 }
192 fprintf(handle3," flow-angle %f\n",flowangle);
193 sprintf(buffer,"temperature");
194 if(getValuNr(buffer)>-1)
195 {
196 fprintf(handle3," turbulent-intensity 0.06\n");
197 fprintf(handle3," turbulent-viscosity %e\n", dynVisc*10.);
198 }
199 }
200 }
201
202 // define precondition only if the pressure is known
203 sprintf(buffer,"pressure");
204 if(getValuNr(buffer)>-1)
205 {
206 fprintf(handle3," precondition\n");
207 fprintf(handle3," inviscid on\n");
208 fprintf(handle3," viscous on\n");
209 sprintf(buffer,"pressure");
210 if(getValuNr(buffer)>-1) fprintf(handle3," pressure %s\n",value[getValuNr(buffer)].string);
211 if(getValuNr(buffer)>-1) fprintf(handle3," velocity %f\n",vel);
212 }
213
214 fclose(handle3);
215
216 return (1);
217 }
218
219