On debugging a Womersley Profile UDF
 

Hi all, I am experiencing difficulties with the attached UDF. I found this UDF on the web. It's designed to calculate the velocity profile at the inlet for a pulsatile flow. I am having trouble compiling this succesfully. When doing so I get the following error message:

cpp -IC:\FLUENT.INC\fluent6.1.22/src -IC:\FLUENT.INC\fluent6.1.22/cortex/src -IC:\FLUENT.INC\fluent6.1.22/client/src -IC:\FLUENT.INC\fluent6.1.22/multiport/src -I. -DUDFCONFIG_H="<udfconfig.h>" C:\libudf\womer-sine.c Error: C:\libudf\womer-sine.c: line 56: Complex: too many arguments supplied (argument 1)

1 file(s) copied.

1 file(s) copied. (system "move user_nt.udf test\ntx86\3d")0 (system "copy C:\FLUENT.INC\fluent6.1.22\src\makefile_nt.udf test\ntx86\3d\makefile") 1 file(s) copied. 0 (chdir "test")() (chdir "ntx86\3d")() womer-sine.c Command line error D2027 : cannot execute 'C:\Program Files\Microsoft Visual Studio\VC98\bin\c1.dll'

Done. ------------------------------------------------------ I would be greatful if somebody could provide some pointers/reference to overcome this problem as I have dealt with complex numbers in c.

Best Regards, David

Here is the code itself ------------------------------------------------------- #include "udf.h"

/*** typedefs and defines you will need ***/ typedef struct DCOMPLEX {double r,i;} dcomplex; #ifndef PI #define PI 3.1415926535 #endif double womervel(double alpha, int nfour, double *ccoef, double *scoef, double r, double t); double Cabs(); dcomplex zbes(int n, dcomplex y); dcomplex Cadd(); dcomplex Csub(); dcomplex Cmul(); dcomplex Complex(); dcomplex Cdiv(); dcomplex RCmul();

DEFINE_PROFILE(womersley, thread, position) { double p[ND_ND]; double r, d, x, y, z, t; face_t f; double scoef[10] = {0.0,0.0}; double ccoef[10] = {4.5128,-3.7549}; begin_f_loop(f, thread) { real q = RP_Get_Real("flow-time"); t = q/20.; F_CENTROID(p, f, thread); x = p[0]; y = p[1]; z = p[2]; r = sqrt(z*z + y*y); d = r/0.0254; F_PROFILE(f, thread, position) = 0.01*womervel(5.3, 2, ccoef, scoef, d, t); } end_f_loop(f, thread) } /************************************************** ******* * returns velocity given: * alpha = Womersley number * nfour = number of fourier coefficients * ccoef = vector of n=0..nfour cosine fourier coefficients * scoef = vector of n=0..nfour sine fourier coefficients * r = radius (must be between 0 and 1) * t = time fraction (must be between 0 and 1) ************************************************** ********/ double womervel(double alpha, int nfour, double *ccoef, double *scoef, double r, double t) { dcomplex zi,z1; double vel; double kt; int k; dcomplex za, zar, zJ0, zJ0r, zJ1, zf10, zq, zvel, zcoef, zJ1J0, zexpt, zqf10, zJ0rJ0; zi = Complex(0.0,1.0); z1 = Complex(1.0,0.0); vel = 2*ccoef[0]*(1-r*r); for (k=1;k<=nfour;k++) { kt = 2.0*PI*k*t; za = RCmul(alpha*sqrt(k)/sqrt(2),Complex(-1.0,1.0)); zar = RCmul(r,za); zJ0 = zbes(0,za); zJ0r = zbes(0,zar); zJ1 = zbes(1,za); zJ1J0 = Cdiv(zJ1,zJ0); zf10 = RCmul( 2.0, Cdiv( zJ1J0,za ) ); zcoef = Complex(ccoef[k],-scoef[k]); zexpt = Complex(cos(kt),sin(kt)); zq = Cmul(zcoef,zexpt); zqf10 = Cdiv(zq,Csub(z1,zf10)); zJ0rJ0 = Cdiv(zJ0r,zJ0); zvel = Cmul(zqf10,Csub(z1,zJ0rJ0)); vel += zvel.r; } return vel; } /* the following routines are from Numerical Recipes in C */ dcomplex Cadd(a,b) dcomplex a,b; { dcomplex c; c.r=a.r+b.r; c.i=a.i+b.i; return c; } dcomplex Csub(a,b) dcomplex a,b; { dcomplex c; c.r=a.r-b.r; c.i=a.i-b.i; return c; } dcomplex Cmul(a,b) dcomplex a,b; { dcomplex c; c.r=a.r*b.r-a.i*b.i; c.i=a.i*b.r+a.r*b.i; return c; } dcomplex Complex(re,im) double re,im; { dcomplex c; c.r=re; c.i=im; return c; } dcomplex Cdiv(a,b) dcomplex a,b; { dcomplex c; double r,den; if (fabs(b.r) >= fabs(b.i)) { r=b.i/b.r; den=b.r+r*b.i; c.r=(a.r+r*a.i)/den; c.i=(a.i-r*a.r)/den; } else { r=b.r/b.i; den=b.i+r*b.r; c.r=(a.r*r+a.i)/den; c.i=(a.i*r-a.r)/den; } return c; } double Cabs(z) dcomplex z; { double x,y,ans,temp; x=fabs(z.r); y=fabs(z.i); if (x == 0.0) ans=y; else if (y == 0.0) ans=x; else if (x > y) { temp=y/x; ans=x*sqrt(1.0+temp*temp); } else { temp=x/y; ans=y*sqrt(1.0+temp*temp); } return ans; } dcomplex RCmul(x,a) double x; dcomplex a; { dcomplex c; c.r=x*a.r; c.i=x*a.i; return c; } /* this is something we wrote to compute the nth order Bessel * function given a complex argument */ dcomplex zbes(int n, dcomplex y) { dcomplex z,zarg,zbes; int i; zarg = RCmul(-0.25,Cmul(y,y)); z = Complex(1.0,0.0); zbes = Complex(1.0,0.0); i = 1; while (Cabs(z)>1e-20 && i<=10000) { z = Cmul(z,RCmul(1.0/i/(i+n),zarg)); if (Cabs(z)<=1.e-20) break; zbes = Cadd(zbes,z); i++; } zarg = RCmul(0.5,y); for (i=1;i<=n;i++) zbes = Cmul(zbes,zarg); return zbes; }

Hello David,
maybe the answer is 6 years too late but the code can be compiled quite well. Your error emerges when you try to interpret with FLUENT. Use the compile fuction and everything will work.

Greetings.

sir can you please share the literature from where/which you developed the above code...deepaknitk2015@gmail.com
thank you so much

Hello Bephi,

maybe mine is also 6 years too late, However, the code compiled well, but I think there is a problem inside.

I have got this numerical error in the first iteration.

Error: Divergence detected in AMG solver: x-momentum

and same error in initializing.

Regards
Mahdi