|
[Sponsors] |
September 9, 2018, 16:00 |
Runge kutta 4th order: Cylinder
|
#1 |
Member
South Yorkshire
Join Date: May 2018
Posts: 33
Rep Power: 8 |
Hi,
I wrote a UDF to calculate the forces of a 2D cylinder subjected to an incoming flow, determine the cylinder's response (acceleration, velocity, and displacement) and prescribe a rigid motion on the cylinder. The response of the cylinder is determined by solving the dynamic equation: F = my'' + cy' + ky The Newmark-beta numerical method seems to work but I want to compare it with another method called runge-kutta 4th order. Considering the same modelling conditions, this numerical method crashes after ~40 iterations due to extremely high displacement values. I can't figure it out where is the problem. Maybe something related to C, where I am only starting to use it. Here is my implementation of the RK4 method: Code:
/*Runge Kutta 4th order*/ #if!RP_HOST func1 = vy; func2 = (fy*water_depth - c*vy - k*y) / total_mass; k0 = dtime * func1; m0 = dtime * func2; k1 = dtime * (func1 + m0*0.5); m1 = dtime * (func2 - (c/total_mass)*0.5*m0 - (k/total_mass)*0.5*k0); k2 = dtime * (func1 + m1*0.5); m2 = dtime * (func2 - (c/total_mass)*0.5*m1 - (k/total_mass)*0.5*k1); k3 = dtime * (func1 + m2); m3 = dtime * (func2 - (c/total_mass)*m2 - (k/total_mass)*k2); y = y + (k0 + 2*k1 + 2*k2 + k3)/6; vy = vy + (m0 + 2*m1 + 2*m2 + m3)/6; func1 = func2 = 0.0; k0 = k1 = k2 = k3 = 0.0; m0 = m1 = m2 = m3 = 0.0; #endif node_to_host_real_4(fy,cg[1],y,vy); Code:
#include "udf.h" #include "mem.h" #define PI 3.141592654 #define zoneID 8 FILE *out; static real y = 0.0; static real dy = 0.0; static real vy = 0.0; static real ay = 0.0; static real current_time = -1; DEFINE_CG_MOTION(cylinderRK4_mesh_MC_3,dt,vel,omega,time,dtime) { real ctime = RP_Get_Real("flow-time"); real ctimestep = RP_Get_Integer("time-step"); if (current_time < ctimestep) { current_time = ctimestep; real k0,k1,k2,k3; real m0,m1,m2,m3; real func1,func2; /*Define variables*/ real cg[3], vcg[3]; real fy,fvy; real Keff,Reff,vprev; real time1, niter; real cg_pos; /*Cylinder variables*/ real density = 998.2; real length = 1; real water_depth = 1; real diameter = 0.01; real fn = 14; real mass = density*pow((0.5*diameter),2)*PI*length; real ad_mass = 1000*PI*pow((0.5*diameter),2)*water_depth; real total_mass = mass + ad_mass; real k = 4*pow(PI*fn,2)*total_mass; real c = 2 * 0.04 * sqrt(k*total_mass); fy = 0.0; fvy = 0.0; int i; real ctime = RP_Get_Real("flow-time"); real ctimestep = RP_Get_Integer("time-step"); #if !RP_HOST Thread *tc,*thread; Domain *d = Get_Domain(1); face_t f; thread = DT_THREAD(dt); tc = Lookup_Thread(d,zoneID); NV_S(vel, =, 0.0); NV_S(omega, =, 0.0); real NV_VEC(A); #endif /*Force calculation*/ #if !RP_HOST begin_f_loop(f,tc) { if (PRINCIPAL_FACE_P(f,tc)) { fvy = F_STORAGE_R_N3V(f,tc,SV_WALL_SHEAR)[1]*-1; F_AREA(A,f,tc); /*Force calculation with a depth of 1m*/ fy += F_P(f,tc)*A[1] + fvy; } } end_f_loop(f,tc) #endif #if!RP_HOST for (i=0;i<3;i++) { cg[i]=DT_CG(dt)[i]; vcg[i] = DT_VEL_CG(dt)[i]; } #endif #if RP_NODE fy = PRF_GRSUM1(fy); /*cg_pos = PRF_GRSUM1(cg[1]);*/ #endif /*Runge Kutta 4th order*/ #if!RP_HOST func1 = vy; func2 = (fy*water_depth - c*vy - k*y) / total_mass; k0 = dtime * func1; m0 = dtime * func2; k1 = dtime * (func1 + m0*0.5); m1 = dtime * (func2 - (c/total_mass)*0.5*m0 - (k/total_mass)*0.5*k0); k2 = dtime * (func1 + m1*0.5); m2 = dtime * (func2 - (c/total_mass)*0.5*m1 - (k/total_mass)*0.5*k1); k3 = dtime * (func1 + m2); m3 = dtime * (func2 - (c/total_mass)*m2 - (k/total_mass)*k2); y = y + (k0 + 2*k1 + 2*k2 + k3)/6; vy = vy + (m0 + 2*m1 + 2*m2 + m3)/6; func1 = func2 = 0.0; k0 = k1 = k2 = k3 = 0.0; m0 = m1 = m2 = m3 = 0.0; #endif node_to_host_real_4(fy,cg[1],y,vy); vel[0] = 0.0; vel[1] = vy; } vel[1] = vy; } Can anyone help me with my code? Thank you. |
|
September 9, 2018, 20:35 |
|
#2 |
Senior Member
|
Is there any particular reason that you exclude pressure contribution from the force calculation? Is time step the same for the Newmark-beta and the Runge-Kutta scheme? The latter might need a smaller time step.
|
|
September 10, 2018, 05:02 |
|
#3 | |
Member
South Yorkshire
Join Date: May 2018
Posts: 33
Rep Power: 8 |
Quote:
F_P(f,tc)*A[1] Before this, I considered a fixed cylinder and obtained the same force values from my UDF and Fluent (report - force in Y direction). The pressure reference value is zero. Yeah, unfortunately, I tried with the same and smaller time step values for the Runge-Kutta scheme with the same results. The weird thing is that I implemented this scheme in matlab at first and works fine. |
||
September 10, 2018, 06:37 |
|
#4 |
Senior Member
|
You are right. I overlooked that line completely.
I think your implementation is right. I suggest that you can reduce the amplitude of oscillatory motion by setting the initial value of position to the stationary point at the first time step, i.e., at first time step. Also, you can monitor the values of (y,y',y'') with R-K scheme during the Newmark-beta simulation, and check out where the difference come from. |
|
September 11, 2018, 06:37 |
|
#5 |
Member
South Yorkshire
Join Date: May 2018
Posts: 33
Rep Power: 8 |
Quote:
After implementing your suggestion I confirmed that both numerical schemes are correctly implemented. Both achieve similar displacements. I think the main problem with my UDF is that it is called at the first iteration of a given time step. I confirm this adding N_ITER inside the if loop. This means that the forces are being calculated with "rough" estimations of the flow. This means that the cylinder will move as a function of a force that is not converged. Is it possible to possible to call the UDF at the last iterations at each time step? Regards, |
|
January 24, 2020, 03:47 |
Runge kutta
|
#7 |
New Member
FDE
Join Date: Feb 2012
Posts: 4
Rep Power: 14 |
Hi,
I need a UDF code to solve the ODE equation (for example: y'=y(y-1), y(0)=...) in FLUENT. I know that the ADJUST macro must be used. But, I have some problems in writing this code. Thanks. |
|
January 27, 2020, 05:55 |
|
#8 |
Senior Member
Alexander
Join Date: Apr 2013
Posts: 2,363
Rep Power: 34 |
show your code
__________________
best regards ****************************** press LIKE if this message was helpful |
|
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
4th order tensor when deriving shear stress tensor in NS equations | TurbJet | Main CFD Forum | 2 | February 7, 2018 11:11 |
Runge Kutta 4th Order Source Code | sugu | Main CFD Forum | 4 | October 26, 2012 03:15 |
Runge Kutta Method | CFDtoy | Main CFD Forum | 12 | May 22, 2005 13:00 |
Navier Stokes - Runge Kutta | CFDtoy | Main CFD Forum | 3 | July 7, 2004 14:13 |
Diagonally Dominate Runge Kutta Method | Anthony Iannetti | Main CFD Forum | 0 | January 23, 2001 21:27 |