UDF to compute forces and moments (torques) on surface zones and print results to file
|
Currently, Fluent cannot print force/moment data to file, only to screen. SYNOPSIS ========= * UDF and SCHEME interface have been added to compute surface forces/moments and print results to file and/or to the GUI. * The UDF is a DEFINE_ON_DEMAND type, so that it is always executed at the end of an iteration. * This UDF has been tested in SERIAL AND PARALLEL and matches FLUENT'S output exactly. CONTENTS ========== Instructions, 1 UDF and 1 Scheme file... sample case/data files available from author INSTRUCTIONS ============== 1. Compile UDF 2. Load case/data 3. Read scheme file fm.scm This adds a component under Solve->Monitors->Force and Moment Data.... 4. Load the UDF 5. Set up case. 6. Under Solve->Monitors->Force and Moment Data.... Pick Print/ Write to send output of Force and Moment computation to GUI and/or files (force.out and moment.out), respectively. Data files are always appended! NOTE: Computation is done upon execution of the EXECUTE_ON_DEMAND function "force_and_moment" Pick surfaces (pick walls only!) Enter coordinates for moment center Click OK 6. Iterate 7. Click EXECUTE_ON_DEMAND... (force_and_moment) to see tresults (all quantities in SI units). OUTPUT FORMAT +++++++++++++++ Forces: Time Fp(x) Fp(y) Fp(z) Fv(x) Fv(y) Fv(z) Ft(x) Ft(y) Ft(z) F = Force p = pressure v = viscous t = total (x) = x-component (y) = y-component (z) = z-component This data is sent to "force.out" if the Check Box "Write to File" is checked in the GUI panel Moments: Same as for forces.... Quantities are in N-m. Output sent to file "moment.out". NOTE: output file have been hardwired into UDF. This is easily modified. /*************************/ /* UDF SOURCE CODE fm.c */ /*************************/ /* FLUENT UDF to calculate forces and moments about specified point */ /* Uses scheme file: fm.scm */ /* Last tested on FLUENT v6.0.20 */ /* AUTHOR: William Wangard, Ph.D. */ /* Last Modified: January 23, 2003 */ #include "udf.h" void fm_cross(real a[], real b[], real c[]); #define FORCE_OUTPUT_FILE "force.out" #define MOMENT_OUTPUT_FILE "moment.out" #define MAX_THREADS 200 #define pi 3.1415926 DEFINE_ON_DEMAND(force_and_moment) /* use define function hooks panel in fluent */ { int i, n, n_threads, thread_n_id, thread_id[MAX_THREADS]; cxboolean fm_write_p, fm_print_p; Domain *domain = Get_Domain(1); Thread *thread = NULL; FILE *force_output_file; FILE *moment_output_file; real r[ND_3], origin[ND_3], xc[ND_3]; real Amag, Atot=0., A_thread_part=0.0; real Force[ND_3], Force_viscous[ND_3], Force_pressure[ND_3], Force_viscous_tot[ND_3], Force_pressure_tot[ND_3], Force_tot[ND_3], A[ND_3]; real Moment[ND_3], Moment_viscous[ND_3], Moment_pressure[ND_3], Moment_tot[ND_3], Moment_viscous_tot[ND_3], Moment_pressure_tot[ND_3]; face_t f; /* get true or false value from solver - panel inputs for print and write */ #if !RP_NODE fm_write_p = RP_Get_Boolean("fm/write?"); fm_print_p = RP_Get_Boolean("fm/print?"); #endif host_to_node_boolean_1(fm_write_p); host_to_node_boolean_1(fm_print_p); /* get number of surfaces highlighted in panel */ #if !RP_NODE n_threads= RP_Get_List_Length("fm/surfaces"); for(n=0; n<n_threads; n++) { thread_id[n] = RP_Get_List_Ref_Int("fm/surfaces",n); } #endif host_to_node_int_1(n_threads); for(n=0; n<n_threads; n++) { host_to_node_int_1(thread_id[n]); } if (!(force_output_file = fopen(FORCE_OUTPUT_FILE,"a"))) Error("unable to open file"); if (!(moment_output_file = fopen(MOMENT_OUTPUT_FILE,"a"))) Error("unable to open file"); /* Initialize forces and moments */ N3V_S(Force,=,0.0); N3V_S(Force_tot,=,0.0); N3V_S(Force_viscous,=,0.0); N3V_S(Force_pressure,=,0.0); N3V_S(Force_viscous_tot,=,0.0); N3V_S(Force_pressure_tot,=,0.0); N3V_S(Moment,=,0.0); N3V_S(Moment_tot,=,0.0); N3V_S(Moment_viscous,=,0.0); N3V_S(Moment_pressure,=,0.0); N3V_S(Moment_viscous_tot,=,0.0); N3V_S(Moment_pressure_tot,=,0.0); N3V_S(origin,=,0.0); N3V_S(xc,=,0.0); N3V_S(r,=,0.0); N3V_S(A,=,0.0); /* get real value from solver - panel inputs for moment center */ #if !RP_NODE origin[0] = RP_Get_Real("moment/centerx"); origin[1] = RP_Get_Real("moment/centery"); origin[2] = RP_Get_Real("moment/centerz"); #endif for(i=0; i<3; i++) { host_to_node_real_1(origin[i]); } Message0("n ...Computing moments around (%10.3e, %10.3e, %10.3e)n", origin[0], origin[1], origin[2]); #if RP_NODE || !PARALLEL for (n=0; n<n_threads;++n) { /* get id's of threads highlighted in panel */ thread_n_id = thread_id[n]; /* get data for wall n */ thread = Lookup_Thread(domain, thread_n_id); A_thread_part = 0.0; begin_f_loop_int(f,thread) { F_AREA(A,f,thread); Amag = N3V_MAG(A); A_thread_part += Amag; } end_f_loop_int(f,thread); /* Message(" partition %i, thread %i, Area = %10.3en", myid, thread_n_id, A_thread_part); */ Atot += A_thread_part; } Atot = PRF_GRSUM1(Atot); if(rp_axi) Atot = 2*pi*PRF_GRSUM1(Atot); Message0(" ...Total area of selected threads is %10.3e (m^2)n", Atot); /* DO MOMENTS AND FORCES */ /* calculate Mx, My and Mz for each of the selected surfaces (fm/surfaces) */ if (fm_write_p || fm_print_p) { for (n=0; n<n_threads;++n) { /* get id's of threads highlighted in panel */ thread_n_id = thread_id[n]; /* get data for wall n */ thread = Lookup_Thread(domain, thread_n_id); begin_f_loop_int(f,thread) { F_AREA(A,f,thread); F_CENTROID(xc,f,thread); N3V_VV(r,=,xc,-,origin); N3V_S(Moment_viscous,=,0.0); N3V_S(Moment_pressure,=,0.0); /* Get wall shear force and pressure force on face */ N3V_VS(Force_viscous,=,F_STORAGE_R_N3V(f,thread,SV_WALL_SHEAR),*,-1.0); N3V_VS(Force_pressure, =, A, *, F_P(f,thread)); if (rp_axi) { Force_pressure[0] = 2*pi*Force_pressure[0]; Force_viscous[0] = 2*pi*Force_viscous[0]; Force_pressure[1] = 0.0; Force_viscous[1] = 0.0; } /* Compute viscous and pressure moments */ if (rp_axi) { Moment_viscous[2] = origin[1]*Force_viscous[0]; Moment_pressure[2] = origin[1]*Force_pressure[0]; } else { fm_cross(Moment_viscous,r,Force_viscous); fm_cross(Moment_pressure,r,Force_pressure); } /* Add forces */ N3V_VV(Force, =, Force_viscous, +, Force_pressure); N3V_VV(Moment,=, Moment_viscous, +, Moment_pressure); /* Accumulate viscous/pressure forces */ N3V_V(Force_viscous_tot, +=, Force_viscous); N3V_V(Force_pressure_tot, +=, Force_pressure); N3V_V(Force_tot, +=, Force); N3V_V(Moment_viscous_tot, +=, Moment_viscous); N3V_V(Moment_pressure_tot, +=, Moment_pressure); N3V_V(Moment_tot, +=, Moment); } end_f_loop_int(f,thread); } /* Accumulate terms over all COMPUTE Nodes */ /* This function also passes global summation value to each compute node */ for(i=0;i<3;i++) { Force_viscous_tot[i] = PRF_GRSUM1(Force_viscous_tot[i]); Force_pressure_tot[i] = PRF_GRSUM1(Force_pressure_tot[i]); Force_tot[i] = PRF_GRSUM1(Force_tot[i]); Moment_viscous_tot[i] = PRF_GRSUM1(Moment_viscous_tot[i]); Moment_pressure_tot[i] = PRF_GRSUM1(Moment_pressure_tot[i]); Moment_tot[i] = PRF_GRSUM1(Moment_tot[i]); } } /* PRINT RESULTS to file ? */ if (fm_print_p) { Message0(" time Fp(x) Fp(y) Fp(z) Fv(x) Fv(y) Fv(z) Ft(x) Ft(y) Ft(z)n"); Message0("%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3en", RP_Get_Real("flow-time"), Force_pressure_tot[0], Force_pressure_tot[1], Force_pressure_tot[2], Force_viscous_tot[0], Force_viscous_tot[1], Force_viscous_tot[2], Force_tot[0], Force_tot[1], Force_tot[2]); Message0(" time Mp(x) Mp(y) Mp(z) Mv(x) Mv(y) Mv(z) Mt(x) Mt(y) Mt(z)n"); Message0("%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3en", RP_Get_Real("flow-time"), Moment_pressure_tot[0], Moment_pressure_tot[1], Moment_pressure_tot[2], Moment_viscous_tot[0], Moment_viscous_tot[1], Moment_viscous_tot[2], Moment_tot[0], Moment_tot[1], Moment_tot[2]); } /* Write data to file ? */ if (fm_write_p) { /* Equivalent to Message0 for file output */ #if RP_NODE if (I_AM_NODE_ZERO_P) #endif { fprintf(force_output_file,"%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3en", RP_Get_Real("flow-time"), Force_pressure_tot[0], Force_pressure_tot[1], Force_pressure_tot[2], Force_viscous_tot[0], Force_viscous_tot[1], Force_viscous_tot[2], Force_tot[0], Force_tot[1], Force_tot[2]); fprintf(moment_output_file,"%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3en", RP_Get_Real("flow-time"), Moment_pressure_tot[0], Moment_pressure_tot[1], Moment_pressure_tot[2], Moment_viscous_tot[0], Moment_viscous_tot[1], Moment_viscous_tot[2], Moment_tot[0], Moment_tot[1], Moment_tot[2]); } } #endif /* RP_NODE || !PARALLEL */ fclose(force_output_file); fclose(moment_output_file); } void fm_cross(real a[], real b[], real c[]) { /* Returns a = b x c for 3-vectors ONLY */ a[0] = b[1]*c[2] - b[2]*c[1]; a[1] = b[2]*c[0] - b[0]*c[2]; a[2] = b[0]*c[1] - b[1]*c[0]; return; } /********** END OF UDF *****************/ ;;;;;;;;;;;;;;;;;;;;;;;; ; SCHEME FILE fm.scm ;;;;;;;;;;;;;;;;;;;;;;; ;; Panel for moment output data ;; works with udf fm.c ;; (define (make-new-rpvar name default type) (if (not (rp-var-object name)) (rp-var-define name default type #f))) ;; ;; Create rpvars for udf, assign initial (ie: default) values ;; types: real, boolean, integer, thread-list ;; () is empty list, #f is false (off) ;; (make-new-rpvar 'fm/print? #f 'boolean) (make-new-rpvar 'fm/write? #f 'boolean) (make-new-rpvar 'fm/surfaces '() 'thread-list) (make-new-rpvar 'moment/centerx 0 'real) (make-new-rpvar 'moment/centery 0 'real) (make-new-rpvar 'moment/centerz 0 'real) ;; ;; Create a panel for the udf ;; (define gui-define-output-data-panel ;; "panel" is a variable name ;; panel #f - no panel made yet (let ((panel #f) ;; declare variables - assign values later (fm/box) (fm/print?) (fm/write?) (fm/surfaces) (moment/centerx) (moment/centery) (moment/centerz) ) ;;; ;;; Function which is called when panel is opened ;;; sets panel variables to values in make-new-rpvar section ;;; default: all options off, no surfaces selected ;;; (define (update-cb . args) ;update panel fields (cx-set-toggle-button fm/print? (rpgetvar 'fm/print?)) (cx-set-toggle-button fm/write? (rpgetvar 'fm/write?)) (cx-set-real-entry moment/centerx (rpgetvar 'moment/centerx)) (cx-set-real-entry moment/centery (rpgetvar 'moment/centery)) (cx-set-real-entry moment/centerz (rpgetvar 'moment/centerz)) ;; get list of surfaces (face threads) (cx-set-list-items fm/surfaces (map thread-name (sort-threads-by-name (get-face-threads)))) ) ;;; ;;; Function which is called when OK button is hit. ;;; assigns variable values that were set in panel after clicking "ok" ;;; (define (apply-cb . args) (rpsetvar 'fm/print? (cx-show-toggle-button fm/print?)) (rpsetvar 'fm/write? (cx-show-toggle-button fm/write?)) (rpsetvar 'moment/centerx (cx-show-real-entry moment/centerx)) (rpsetvar 'moment/centery (cx-show-real-entry moment/centery)) (rpsetvar 'moment/centerz (cx-show-real-entry moment/centerz)) (rpsetvar 'fm/surfaces (map thread-name->id (cx-show-symbol-list-selections fm/surfaces))) ) (lambda args ;;; if panel does not exist, make panel (if (not panel) (let ((table)) (set! panel (cx-create-panel "Force and Moment Settings" apply-cb update-cb)) ;; create table w/o title in panel to enable row by column format (set! table (cx-create-table panel "" 'border #f 'below 0 'right-of 0)) ;;; ;;; Create fm inputs. ;;; (set! fm/box (cx-create-table table "Force and Moment Options" 'row 0)) (set! fm/print? (cx-create-toggle-button fm/box "Print Results" 'col 1 'row 1)) (set! fm/write? (cx-create-toggle-button fm/box "Write Results to file" 'col 1 'row 2)) (set! fm/surfaces (cx-create-symbol-list fm/box "Boundary Zones" 'visible-lines 10 'multiple-selections #t 'row 3 'col 1)) (cx-set-list-items fm/surfaces (map thread-name (sort-threads-by-name (get-face-threads)))) (set! moment/centerx (cx-create-real-entry fm/box "x-moment-center" 'row 4 'col 1)) (set! moment/centery (cx-create-real-entry fm/box "y-moment-center" 'row 4 'col 2)) (set! moment/centerz (cx-create-real-entry fm/box "z-moment-center" 'row 4 'col 3)) )) ;; tell solver to display the panel (cx-show-panel panel)))) ;;; ;;; Add new panel to end of User-Defined menu. ;;; no need to modify here except panel name and menu location ;;; (cx-add-item "Monitors" "Force and Moment Data..." #U #f cx-client? gui-define-output-data-panel) ;;;;; END OF SCHEME FILE |
||
|