Fluent:UDF to compute forces and momentson surface zones and print to file

2017-03-21 by:CAE仿真在線來源:互聯(lián)網(wǎng)

Problem:Currently, Fluent cannot print force/moment data to file, only to screen. Resolution: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.3e\n", 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.3e\n",
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.3e\n",
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.3e\n",
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.3e\n",
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

開放分享：優(yōu)質(zhì)有限元技術(shù)文章,助你自學(xué)成才