# The problem of inner iteration when calling the subroutine

 Register Blogs Members List Search Today's Posts Mark Forums Read

 October 14, 2010, 21:29 The problem of inner iteration when calling the subroutine #1 New Member   Join Date: Oct 2010 Posts: 6 Rep Power: 6 To facilitate the discussion，I open a new theme. Moving the problem description to below: I am simulating the starting behavior of a domain including Impeller, using a Fortran subroutine.The subroutine is used to update the angular velocity of the domain.And I calculate the angular velocity taking into account the torque generated by the blade parts and the inertia of the impeller, and this value is used by the cfx solver to set the angular velocity of the blade. The subroutine is used to calculate angular velocity of the rotational domain in the timestep n by using the value of previous time step.The problem is, when the calculation is proceeding,although the cycle is set to update angular velocity at the start of a new timestep, the subroutine is called many times per inner iteration,not once every timestep. That is angular velocity has been updated dozens of times within each iteration. I want to what is the reason of it? And how to solve this problem? The error message shown below: ERROR #001100279 has occurred in subroutine ErrAction. | | Message: | | Angular velocity can depend only on time for transient simulations.

 October 14, 2010, 21:37 #2 Senior Member   Michael P. Owen Join Date: Mar 2009 Posts: 195 Rep Power: 8 The error doesn't have anything to do with you problem description. The angular velocity specified for a rotating domain cannot be a function of anything except time. This is a well known limitation of CEL. You will have to be more crafty about how you do this. Also, FORTRAN is not required.

 October 14, 2010, 21:44 #3 Senior Member   Michael P. Owen Join Date: Mar 2009 Posts: 195 Rep Power: 8 By the way, if you have R12.1 you can do this using the 6DOF rigid body solver. Otherwise you'll have to us technique like the one I developed here: http://www.youtube.com/watch?v=md8vXAUn7g0 http://www.youtube.com/watch?v=9uAobJm0r1I These were done in v11, before the 6DOF solver, without FORTRAN.

 October 14, 2010, 22:42 #4 New Member   Join Date: Oct 2010 Posts: 6 Rep Power: 6 Thanks for your reply.But i have not the R12.1,also i can not open the links you give me.could you give me some description about the techniques you have developed in other ways.Or could you give me some other suggestions?

 October 14, 2010, 23:26 #5 New Member   Join Date: Oct 2010 Posts: 6 Rep Power: 6 Here i apology for posting new topic behind Filippo's

 October 15, 2010, 09:42 #6 Senior Member   Michael P. Owen Join Date: Mar 2009 Posts: 195 Rep Power: 8 The technique I developed, and I'm sure I'm not the first, is to use an Additional Variable to keep track of the angular momentum of the rotating object. You can use CEL expressions to calculate the aero/hydrodynamic torque on the object, include any external torques acting on it, and create a Source for the angular momentum additional variable that represents the net torque. Finally you use CEL expressions to physically rotate the mesh. Note that this only works for transient runs, and that you do NOT use Domain Motion > Rotating and specify an angular velocity; instead, you are physically rotating the mesh.

 October 16, 2010, 05:34 #7 Super Moderator   Glenn Horrocks Join Date: Mar 2009 Location: Sydney, Australia Posts: 10,668 Rep Power: 84 Not sure about this, but can't you do this without the additional variable and moving mesh? Simply a CEL callback function to read the angular velocity (eg ave(Angular Velocity)@RotatingDomain)? Callback functions are evaluated at the start of a timestep so in effect give you the angular velocity of the last timestep. Then you can get the torques from the model and from external influences and work out the angular velocity of the next timestep. And I understand this works for rotating frame of reference too. But I admit I have not tried it, just heard of it being done that way. Hopefully it works, it would be much simpler, easier and more efficient than the method you describe.

October 16, 2010, 15:56
#8
Senior Member

Michael P. Owen
Join Date: Mar 2009
Posts: 195
Rep Power: 8
Quote:
 Originally Posted by ghorrocks Not sure about this, but can't you do this without the additional variable and moving mesh? Simply a CEL callback function to read the angular velocity (eg ave(Angular Velocity)@RotatingDomain)? Callback functions are evaluated at the start of a timestep so in effect give you the angular velocity of the last timestep. Then you can get the torques from the model and from external influences and work out the angular velocity of the next timestep. And I understand this works for rotating frame of reference too. But I admit I have not tried it, just heard of it being done that way. Hopefully it works, it would be much simpler, easier and more efficient than the method you describe.
Glenn,

I'm afraid that won't work on a couple of accounts. One, the angular velocity of a rotating domain cannot be a function of any variable except time, and you cannot call an integrative function of a non-field variable.

 October 19, 2010, 17:46 #9 Senior Member   Edmund Singer P.E. Join Date: Aug 2010 Location: Minneapolis, MN Posts: 474 Rep Power: 10 What Glenn describes is correct. This can be done entirely within CEL as he describes. But there are advantages to putting this in as an additional variable. In the event your motion doesn't simply rotate a subdomain, but instead it moves the mesh, in the event of a remesh incident, you will need the information of angular velocity to be interpolated onto the new mesh file so that you can correct keep the momentum.

 October 19, 2010, 18:13 #10 Senior Member   Michael P. Owen Join Date: Mar 2009 Posts: 195 Rep Power: 8 Try it then.

 October 19, 2010, 18:15 #11 Senior Member   Edmund Singer P.E. Join Date: Aug 2010 Location: Minneapolis, MN Posts: 474 Rep Power: 10 Have already. Done it often.

 October 19, 2010, 18:29 #12 Super Moderator   Glenn Horrocks Join Date: Mar 2009 Location: Sydney, Australia Posts: 10,668 Rep Power: 84 Can you post the CCL of an example?

 October 19, 2010, 18:30 #13 Senior Member   Edmund Singer P.E. Join Date: Aug 2010 Location: Minneapolis, MN Posts: 474 Rep Power: 10 Glenn: Which case are you refering to? CEL or AV case.

 October 19, 2010, 18:35 #14 Super Moderator   Glenn Horrocks Join Date: Mar 2009 Location: Sydney, Australia Posts: 10,668 Rep Power: 84 The CEL case.

 October 19, 2010, 19:14 #15 Senior Member   Michael P. Owen Join Date: Mar 2009 Posts: 195 Rep Power: 8 I think that there are certain cases where it can be done entirely with CEL without using an additional variable. If you have a surface of constant radius on which you can evalute areaAve(Mesh Velocity), you can turn this into the angular velocity, and then reference that in the CEL expressions for defining the mesh motion. But again, this involves the method I already described, physically moving the mesh via CEL, and does NOT involve the use of Domain Motion > Rotating, which I repeat will NOT work, because the angular velocity of a rotating frame of reference can only be a function of time.

 October 19, 2010, 20:02 #16 Super Moderator   Glenn Horrocks Join Date: Mar 2009 Location: Sydney, Australia Posts: 10,668 Rep Power: 84 You can get the angular velocity with something like maxVal(Mesh Velocity)@rotatingdomain where the maximum radius point is known. This can be worked back to the angular velocity for any geometry.

October 19, 2010, 20:22
#17
Senior Member

Michael P. Owen
Join Date: Mar 2009
Posts: 195
Rep Power: 8
Quote:
 Originally Posted by ghorrocks You can get the angular velocity with something like maxVal(Mesh Velocity)@rotatingdomain where the maximum radius point is known. This can be worked back to the angular velocity for any geometry.
Yeah, I agree. The AV is definitely unnecessary. The mesh motion is not.

 October 20, 2010, 10:36 #18 Senior Member   Edmund Singer P.E. Join Date: Aug 2010 Location: Minneapolis, MN Posts: 474 Rep Power: 10 I see where the confusion is. I glossed over the Domain Motion > Rotating statement. I concur with what has been talked about. For thread completeness, here is the CEL you inquired about. This has wall motion rotating around X0,Y0,Z0, along a hinge point pointing in the global Z direction, but not coincedent. The new displacements are updated using New X, New Y and New Z. LIBRARY: CEL: &replace EXPRESSIONS: Ifin = 127.00 [kg mm^2] New X = ((x-X0)-Total Mesh Displacement X)*cos(alpha)+((y-Y0)-Total Mesh Displacement Y)*sin(alpha)+X0 New Y = -((x-X0)-Total Mesh Displacement X)*sin(alpha)+((y-Y0)-Total Mesh Displacement Y)*cos(alpha)+Y0 New Z = ((z)-Total Mesh Displacement Z) + Z0 X0 = 44.4355 [mm] Y0 = 25.181 [mm] Z0 = move alpha = ((angdisp -(mit+omgold)*tstep*1[s^-1])+alphstart)*1[deg] angdisp = initang-newang back max z = minVal(Global Z Coordinate)@Proj Rear Wall back min z = maxVal(Global Z Coordinate)@Gun Muzzle Exit back move = move*max(0, (min(1,(back min z-z)/ (back min z- (back max z) ) ) )) finforx = (forx1+forx2+forx3+forx4+forx5) finfory = (fory1+fory2+fory3+fory4+fory5) fintorq = (finforx*Y0-finfory*X0+(tor1+tor2+tor3+tor4+tor5)) fintorqreal = (tor1+tor2+tor3+tor4+tor5) forx1 = force_x()@Proj Fin Bottom forx2 = force_x()@Proj Fin Top forx3 = force_x()@Proj Fin Hinge forx4 = force_x()@Proj Fin Side Edges forx5 = force_x()@Proj Fin Outside Edge fory1 = force_y()@Proj Fin Bottom fory2 = force_y()@Proj Fin Top fory3 = force_y()@Proj Fin Hinge fory4 = force_y()@Proj Fin Side Edges fory5 = force_y()@Proj Fin Outside Edge front max z = maxVal(Global Z Coordinate)@Proj Front Wall front min z = minVal(Global Z Coordinate)@Proj Front Wall front move = move* (min(1,(front max z-z)/ (front max z- (front min z + 0.1 [m]) ) ) ) initang = atan((ycod-tmdy)/(xcod-tmdx))/1 [deg]+step(-(xcod-tmdx)/1[m]-0.00001)*180 mit = fintorq/Ifin*tstep*57.296*1[s] newang = atan(ycod/xcod)/1 [deg]+step(-xcod/1[m]-0.00001)*180 omgold = -(((areaAve(Mesh Velocity X)@Proj Fin Outside Edge)^2+(areaAve(Mesh Velocity Y)@Proj Fin Outside Edge)^2)^0.5/((areaAve(Global X Coordinate)@Proj Fin Outside Edge-X0)^2+(areaAve(Global Y Coordinate)@Proj Fin Outside Edge-Y0)^2)^0.5)*57.296*1[s] periodmove = front move*step((z-minVal(Global Z Coordinate)@Proj Rear Wall)/1[m])+back move*step((maxVal(Global Z Coordinate)@Gun Muzzle Exit-z)/1[m]) pinit = pstart*zstartramp*(xystartstep)*(1-underfinstartp)+underfinstartp*5e7[Pa] pstart = 4.5e7[Pa] tempstart = 1422[K] tend = 10000[s] tmdx = areaAve(Total Mesh Displacement X)@Proj Fin Outside Edge tmdy = areaAve(Total Mesh Displacement Y)@Proj Fin Outside Edge tor1 = torque_z()@Proj Fin Bottom tor2 = torque_z()@Proj Fin Top tor3 = torque_z()@Proj Fin Hinge tor4 = torque_z()@Proj Fin Side Edges tor5 = torque_z()@Proj Fin Outside Edge tstart = 0.0 [s] tstep = 1e-8 [s] xcod = areaAve(Global X Coordinate)@Proj Fin Outside Edge-X0 xystartstep = step(0.08 - ((x^2+y^2)/1[m^2])^0.5) ycod = areaAve(Global Y Coordinate)@Proj Fin Outside Edge-Y0 zstartramp = max(0,1-z/0.04[m]) zvstartramp = max(0,1-z/0.0142[m]) END END END

 October 23, 2010, 20:37 #19 New Member   Join Date: Oct 2010 Posts: 6 Rep Power: 6 My problem has been solved.Thanks for all the valuable advices

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post pitisrisuk CFX 1 July 2, 2012 03:36 cincaipatron OpenFOAM Installation 9 January 11, 2010 07:37 Yuqing Feng CFX 0 May 5, 2004 23:39 sudha FLUENT 3 April 28, 2004 08:40 sisilxp FLUENT 0 April 23, 2004 23:36

All times are GMT -4. The time now is 20:04.