Welcome to the TRNSYS Users Forum.
The forum is a place where people can interact and have discussions about different topics involving the use of the TRNSYS software package. Here you can post topics for discussion or questions on using TRNSYS and get advice from other users or TRNSYS experts. This forum is not intended for detailed technical support. Users should contact their distributor’s hotline for such assistance.
Some tips for success on using the forum:
- Follow the Forum Rules posted in Forum Administration.
- There are categories for different types of topics and questions. Post your topic or question into the proper category.
- Before posting a topic or question search the existing topics (and the TRNSYS Users listserv archive or Post archive) to see if a similar topic or question has already been answered.
- Use a descriptive topic name. Don’t use attention getting subjects, they don’t get attention and only annoy people.
- State the version of TRNSYS and which add-ons your are using.
- Include enough specific details for your topic of question to be answered. Just posting “Why am I getting an error?” without describing the specific error and what you are trying to do when you get the error will not receive a response that fixes your issue.
- Remember when people help you, they are doing you a favor. Be patient, help people out by posting good descriptions of what you need help with, and be polite even if a response does not solve your issue.
- Moderators may edit your post for clarity or move your topic to a more appropriate category.
I'm writing new modules using TypeStudio.The mathematical model of my model is a set of first-order equations.The manual describes how to write independent differential equations, but how do I write a set first-order equations? I will give a simplified example set of differential equations: dA/dt=B*exp(C/T)-B*exp(C/T)*A
dT/dt=D*T+E*dA/dt
I tried writing code using the solution of independent differential equations, but the result is incorrect, I will attach this part of my code.
I would really appreciate it if I could receive some feedback!
ai = getDynamicArrayValueLastTimestep(2)
aa = -Frequency_factor*DEXP(-Activation_energy/(Ideal_gas_constant*tBar))
bb = Frequency_factor*DEXP(-Activation_energy/(Ideal_gas_constant*tBar))
Call SolveDiffEq(aa,bb,ai,af,abar)
Call setDynamicArrayValueThisIteration(1,af)
Reaction_rate=Frequency_factor*DEXP(-Activation_energy/(Ideal_gas_constant*tBar))*(1-aBar)
Ds=(1-aBar)*Density_of_S1+aBar*Density_of_S0
Ms=(1-aBar)*Molecular_mass_of_S1+aBar*Molecular_mass_of_S0
CPS=(1-aBar)*Specific_heat_of_S1+aBar*Specific_heat_of_S0
Eff=(1-Porosity_of_sorbent)*Ds*CPs+Porosity_of_sorbent*Vapor_density*Specific_heat_capacity_of_water_vapor
Teff=(1-Porosity_of_sorbent)*Thermal_conductivity_of_sorbent+Porosity_of_sorbent*Steam_thermal_conductivity
U=-(Permeability_of_sorbent/Dynamic_viscosity)*(Partial_pressure_of_steam-Vapor_density*Vapor_density)
ti= getDynamicArrayValueLastTimestep(1)
CC=(Teff-Specific_heat_capacity_of_water_vapor*U*Vapor_density)/Eff
DD=-((Ds/Ms)*Reaction_rate*Reaction_enthalpy)/Eff
Call SolveDiffEq(CC,DD,ti,tf,tbar)
Call setDynamicArrayValueThisIteration(1,tf)
A_Weiss 11/05/2023 6:40 pm
@meiling_j I noticed on line 1 you're using slot 2 in the dynamic storage array to read in ai, but on line 5 you're setting af to slot 1 of the dynamic storage array. Is that why the code isn't working as you expected?
meiling_J Topic starter
12/05/2023 1:13 am
Thank you very much for the answer!
I did what you said, but the result is the same as before, it does not seem to be the problem of the two lines of code.
@meiling_j Where is your initial value for tbar coming from? It is used on lines 2, 3, and 6 of your code, but it isn't defined until the second to last line of the snippet shown here.
What you probably want to do is use successive substitution. Set an initial guess value for tbar, then solve your system of equations to get the new value for tbar. Then, compare the new value of tbar to your initial guess. If the difference is <= some convergence tolerance (say 0.001 C), the system is solved. If not, substitute the new value for the initial guess, solve again using the new guess value, and repeat until the value of tbar is within some convergence tolerance of the guess value.
There are multiple ways to write a repeating loop like this in Fortran, but I would recommend something like the following:
maxiters = 100 iters = 0 tbar_guess = 0 tbar = 20 !C tolerance = 0.001 !C Do While ((DABS( tbar-tbar_guess)> tolerance).and.(iters<maxIters)) iters = iters+1 tbar_guess = tbar <insert equations here - use tbar_guess in place of tbar up until where tbar is calculated> EndDo
meiling_J Topic starter
16/05/2023 8:16 am
@a_weiss Thank you very much for answering me. I think my expression is not clear enough. What I want to ask is how to solve the equations of the two differential equations I mentioned together, that is, equation one: dA/dt=B*exp(C/T)-B*exp(C/T)*A and equation two: dT/dt=D*T+E*dA/dt are binary first-order differential equations. How should I write the code when TRNSYS is developing its own module? Or is there any information in this area for reference?
A_Weiss 16/05/2023 2:07 pm
@meiling_J There's no built-in TRNSYS subroutine for simultaneous solution of multiple differential equations. There might be something in Intel's Math Kernel Library or another external Fortran library you could call, but if so I don't know about it or have any familiarity with it.
The way I would solve the system of equations you listed here is 1)Assign a guess value for T; 2) Call SolveDiffEq(aa, bb, ai, af, abar) to solve for A, using aa = -B*exp(C/T_guess) and bb = B*exp(C/T_guess); 3) Solve dA_dt = B*exp(C/T_guess)-B*exp(C/T_guess)*abar; 4) Call SolveDiffEq(cc, dd, ti, tf, tbar) to solve for T, using cc = D and dd = E*dA_dt; 5) Compare tbar to T_guess, and if they don't agree within some tolerance, set T_guess = tbar, return to step 1, and solve again.
I think the code you have is very close, you just need to assign an initial guess value for tbar, use tbar_guess in place of tbar up until where tbar is solved, and embed the code in a loop that compares and updates the value of tbar_guess until tbar converges with tbar_guess. As long as the guess value is something reasonable and the solution is stable, you should get convergence within a handful of iterations.
The TRNSYS built-in compiler (TypeStudio) doesn't provide guidance on writing TRNSYS code; it only compiles Fortran code into a .dll TRNSYS can read and reports errors if there's a syntax problem or if a subroutine was called that it doesn't recognize. If you write your proforma first in Simulation Studio, then (still in the proforma editing window) select File -> Export As -> Fortran, you’ll get skeleton code for TypeStudio that includes some boilerplate for version signing and reading in parameters and inputs, but that’s about as far as its ‘guidance’ goes.
There is no built-in solution capability or formal documentation in TRNSYS for solving coupled ODEs. However, I’ve written several TRNSYS Types that solve coupled ODEs, and I always use successive substitution as described in the above answer and comments. For example, Type 156 solves coupled ODEs for dT_HX/dt and dT_tank/dt. The code is in the SourceCode/Types folder of the main TRNSYS18 directory. It’s a bit daunting, but here’s the gist:
Line 654: Starts the loop that solves the temperature at each node of the HX.
Lines 655-772: Solves for AA and BB in the equation dT_HX/dt = AA*T_HX + BB. Note that Taverage_TankNode and Taverage_HXNode are used in this block of code. These arrays have initial values set on lines 347 and 351, respectively, and are initialized at each timestep on lines 591-606. They’re basically guess values at this point.
Lines 773-782: Solves the final and average temperature for the heat exchanger node. Here, the analytical solution for the ODEs is spelled out, but calling SolveDiffEq(AA_HX(k), BB_HX(k), Tinitial_HXNode(k), Tfinal_HXNode_New(k), Taverage_HXNode_New(k)) would work as well. Note that the solution is for Taverage_HXNode_New, not Taverage_HXNode!
Lines 785-813: Checks if Taverage_HXNode_New differs from Taverage_HXNode by more than ConvergenceCriteria. If Taverage_HXNode hasn’t converged, it sets Taverage_HXNode = Taverage_HXNode_New for each node, then (see line 812) the code returns to command line 205 (which is actually on line 652, just above where the HX nodal calculations began), and the calculations are carried out again. This can repeat as many times as necessary, until either Taverage_HXNode_New converges with Taverage_HXNode for all nodes or until the maximum number of iterations is reached.
Lines 815-960: Solves for AA and BB in the equation dT_Tank/dt = AA*T_tank + BB. Again, Taverage_TankNode and Taverage_HXNode are used in this block of code.
Lines 962-971: Solves the final and average temperature for the tank node. Note the solution is for Taverage_TankNode_New, not Taverage_TankNode!
Lines 973-1001: Checks if Taverage_TankNode_New differs from Taverage_TankNode by more than ConvergenceCriteria. If Taverage_TankNode hasn’t converged, it sets Taverage_TankNode = Taverage_TankNode_New for each node, then (see line 1000) the code returns to command line 200 (which is actually on line 635, a little ways above where the HX nodal calculations began), and ALL of the calculations are carried out again, including the HX node calculations. This can repeat as many times as necessary, until either Taverage_TankNode_New converges with Taverage_TankNode for all nodes or until the maximum number of iterations is reached.
It looks like a lot of work to implement (I wish I had an example in an actual TRNSYS Type that was shorter!), but the parts that implement successive substitution are relatively short and straightforward. If you’ve tried this scheme and are still having trouble solving your coupled ODEs, you’re welcome to message me at weiss@tess-inc.com with your code, and I can take a look.
TM Abir Ahsan 06/06/2023 5:00 pm
@a_weiss amazing, thank you so much for the detailed explanation. I will follow the instructions and follow up with you on this.
TM Abir Ahsan 17/09/2023 3:21 pm
@a_weiss I was looking at all other types at trnsys source code section and found that tess type 963 (lumped capacitance model) solves differential equation analytically without running it through a successive substitution. I was wondering why? In the last comment your gave very detailed explanation on how you implemented successive substitution (Picard's method I suppose) for the storage tank which is daunting indeed to follow. It would be super helpful if you could give the explanation based on type 963 which is rather simple and concise. I appreciate your help.
