@@ -184,6 +184,139 @@ def write_mixer(mixer, output_folder):
184184 f .write ("\t \t }\n " )
185185
186186
187+ def write_mixer_force_sign (mixer , output_folder ):
188+ with open (os .path .join (output_folder , "fvModels" ), "a+" ) as f :
189+ f .write (f"\t \t source_pt_x={ mixer .x } ;\n " )
190+ f .write (f"\t \t source_pt_y={ mixer .y } ;\n " )
191+ f .write (f"\t \t source_pt_z={ mixer .z } ;\n " )
192+ f .write (f"\t \t disk_rad={ mixer .rad } ;\n " )
193+ f .write ("\t \t disk_area=pi*disk_rad*disk_rad;\n " )
194+ f .write (f"\t \t power={ mixer .power } ;\n " )
195+ f .write (f"\t \t smear_factor={ float (mixer .smear )} ;\n " )
196+ f .write (f"\t \t startTime = { mixer .start_time } ;\n " )
197+ f .write ("\t \t if (time.value() > startTime)\n " )
198+ f .write ("\t \t {\n " )
199+ f .write ("\t \t \t // Get V1\n " )
200+ f .write ("\t \t \t double source_sign_factor = 1.0;\n " )
201+ f .write ("\t \t \t double V1 = 0;\n " )
202+ f .write ("\t \t \t double V2 = 0;\n " )
203+ f .write ("\t \t \t double rhoV;\n " )
204+ f .write ("\t \t \t double dist_tol = disk_rad*3;\n " )
205+ f .write ("\n " )
206+ f .write ("\t \t \t double dist_n;\n " )
207+ f .write ("\t \t \t double upV = 0;\n " )
208+ f .write ("\t \t \t double uprhoV = 0;\n " )
209+ f .write ("\t \t \t double upVvol = 0;\n " )
210+ f .write ("\t \t \t double downV = 0;\n " )
211+ f .write ("\t \t \t double downrhoV = 0;\n " )
212+ f .write ("\t \t \t double downVvol = 0;\n " )
213+ f .write ("\t \t \t double dist2;\n " )
214+
215+ f .write ("\t \t \t forAll(C,i)\n " )
216+ f .write ("\t \t \t {\n " )
217+ f .write (
218+ "\t \t \t \t dist2 = (C[i].x()-source_pt_x)*(C[i].x()-source_pt_x);\n "
219+ )
220+ f .write (
221+ "\t \t \t \t dist2 += (C[i].y()-source_pt_y)*(C[i].y()-source_pt_y);\n "
222+ )
223+ f .write (
224+ "\t \t \t \t dist2 += (C[i].z()-source_pt_z)*(C[i].z()-source_pt_z);\n "
225+ )
226+ f .write ("\n " )
227+ if mixer .normal_dir == 0 :
228+ f .write ("\t \t \t \t dist_n = (C[i].x()-source_pt_x);\n " )
229+ elif mixer .normal_dir == 1 :
230+ f .write ("\t \t \t \t dist_n = (C[i].y()-source_pt_y);\n " )
231+ elif mixer .normal_dir == 2 :
232+ f .write ("\t \t \t \t dist_n = (C[i].z()-source_pt_z);\n " )
233+ f .write ("\n " )
234+
235+ f .write (
236+ "\t \t \t \t if (dist2 < dist_tol*dist_tol && dist_n < -dist_tol/2) {\n "
237+ )
238+ f .write ("\t \t \t \t \t upVvol += V[i] * alphaL[i];\n " )
239+ f .write (
240+ f"\t \t \t \t \t upV += V[i] * alphaL[i] * UL[i][{ int (mixer .normal_dir )} ];\n "
241+ )
242+ f .write ("\t \t \t \t \t uprhoV += V[i] * alphaL[i] * rhoL[i];\n " )
243+ f .write ("\t \t \t \t }\n " )
244+ f .write (
245+ "\t \t \t \t if (dist2 < dist_tol*dist_tol && dist_n > dist_tol/2) {\n "
246+ )
247+ f .write ("\t \t \t \t \t downVvol += V[i] * alphaL[i];\n " )
248+ f .write (
249+ f"\t \t \t \t \t downV += V[i] * alphaL[i] * UL[i][{ int (mixer .normal_dir )} ];\n "
250+ )
251+ f .write ("\t \t \t \t \t downrhoV += V[i] * alphaL[i] * rhoL[i];\n " )
252+ f .write ("\t \t \t \t }\n " )
253+ f .write ("\t \t \t }\n " )
254+ f .write ("\n " )
255+ f .write ("\t \t \t reduce(uprhoV, sumOp<scalar>());\n " )
256+ f .write ("\t \t \t reduce(downrhoV, sumOp<scalar>());\n " )
257+ f .write ("\t \t \t reduce(upV, sumOp<scalar>());\n " )
258+ f .write ("\t \t \t reduce(downV, sumOp<scalar>());\n " )
259+ f .write ("\t \t \t reduce(downVvol, sumOp<scalar>());\n " )
260+ f .write ("\t \t \t reduce(upVvol, sumOp<scalar>());\n " )
261+ f .write ("\n " )
262+ f .write ("\t \t \t downV /= downVvol;\n " )
263+ f .write ("\t \t \t upV /= upVvol;\n " )
264+ f .write ("\t \t \t downrhoV /= downVvol;\n " )
265+ f .write ("\t \t \t uprhoV /= upVvol;\n " )
266+ f .write ("\n " )
267+ if mixer .sign == "+" :
268+ f .write ("\t \t \t source_sign_factor = -1.0;\n " )
269+ f .write ("\t \t \t if (upV >= 0){\n " )
270+ f .write ("\t \t \t \t V1 = 0.0;\n " )
271+ f .write ("\t \t \t } else {\n " )
272+ f .write ("\t \t \t \t V1 = std::abs(upV);\n " )
273+ f .write ("\t \t \t }\n " )
274+ f .write ("\t \t \t rhoV = uprhoV;\n " )
275+ elif mixer .sign == "-" :
276+ f .write ("\t \t \t source_sign_factor = 1.0;\n " )
277+ f .write ("\t \t \t if (downV <= 0){\n " )
278+ f .write ("\t \t \t \t V1 = 0.0;\n " )
279+ f .write ("\t \t \t } else {\n " )
280+ f .write ("\t \t \t \t V1 = std::abs(downV);\n " )
281+ f .write ("\t \t \t }\n " )
282+ f .write ("\t \t \t rhoV = downrhoV;\n " )
283+ # f.write("\t\t\t}\n")
284+ f .write (
285+ '\t \t \t Foam::Info << "[BIRD:DYNMIX INFO] V1 = " << V1 << endl;\n '
286+ )
287+ f .write ("\t \t \t \n " )
288+ f .write ("\t \t \t // Get V2\n " )
289+ f .write ("\t \t \t V2 = findV2(power, rhoV, disk_area, V1);\n " )
290+ f .write ("\n " )
291+ f .write ("\t \t \t forAll(C,i)\n " )
292+ f .write ("\t \t \t {\n " )
293+ f .write (
294+ "\t \t \t \t double Thrust=0.5*rhoL[i]*(V2*V2 - V1*V1)*disk_area;\n "
295+ )
296+ f .write (
297+ "\t \t \t \t double dist2=(C[i].x()-source_pt_x)*(C[i].x()-source_pt_x);\n "
298+ )
299+ f .write (
300+ "\t \t \t \t dist2 += (C[i].y()-source_pt_y)*(C[i].y()-source_pt_y);\n "
301+ )
302+ f .write (
303+ "\t \t \t \t dist2 += (C[i].z()-source_pt_z)*(C[i].z()-source_pt_z);\n "
304+ )
305+
306+ f .write ("\t \t \t \t double epsilon=pow(V[i],0.33333)*smear_factor;\n " )
307+ f .write (
308+ "\t \t \t \t double sourceterm=alphaL[i]*(Thrust/pow(pi,1.5)/pow(epsilon,3.0))*\n "
309+ )
310+ f .write ("\t \t \t \t \t exp(-dist2/(epsilon*epsilon));\n " )
311+
312+ f .write (
313+ f"\t \t \t \t Usource[i][{ int (mixer .normal_dir )} ] -= source_sign_factor*sourceterm*V[i];\n "
314+ )
315+
316+ f .write ("\t \t \t }\n " )
317+ f .write ("\t \t }\n " )
318+
319+
187320def write_end (output_folder ):
188321 with open (os .path .join (output_folder , "fvModels" ), "a+" ) as f :
189322 f .write ("\t #};\n " )
0 commit comments