CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

Diffusion Equations coupled by Boundary Conditions

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

Reply
 
LinkBack Thread Tools Display Modes
Old   March 1, 2016, 12:39
Default Diffusion Equations coupled by Boundary Conditions
  #1
New Member
 
Michael
Join Date: Mar 2016
Posts: 4
Rep Power: 2
Amiga-Freak is on a distinguished road
Hi, I'm new here and not sure whether my question is in the right place.

I have the following problem:

I'd like to solve two diffusion equations, e.g.
\partial_t u - D_1\Delta u = 0,\;\;u = f(w)\; on\; \partial \Omega
and
\partial_t w - D_2\Delta w = 0,\;\; \partial_nw = g(u)\; on\; \partial \Omega

which means the PDEs are coupled by their boundary conditions.
The solution of u depends on w by a Dirichlet BC.
And w depends on u by a Neumann BC.

At the moment I first solve the equation for u, then compute g(u) and afterwards solve the equation for w, which allows me to compute f(w) for the next iteration, etc. ....

This works very nice for the stationary case, i.e. for \partial_t u = \partial_t w = 0
For the time dependent case however....the method is highly unstable.

I'm searching now for a method to solve both equations simultaneously, somehow - instead of one after another.

But I have no clue how to do this. If one searches for literature about this issue, it seems nobody has ever had the problem to solve PDEs coupled by boundary conditions.

Any hint? Any literature? Ideas?

BTW: I'm not working with any commercial software, but a own FEM software written in Fortran. Therefore please don't recommend commercial software.

Greetings,
Amiga-Freak
Amiga-Freak is offline   Reply With Quote

Old   March 1, 2016, 13:06
Default
  #2
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,692
Rep Power: 33
FMDenaro will become famous soon enoughFMDenaro will become famous soon enough
Quote:
Originally Posted by Amiga-Freak View Post
Hi, I'm new here and not sure whether my question is in the right place.

I have the following problem:

I'd like to solve two diffusion equations, e.g.
\partial_t u - D_1\Delta u = 0,\;\;u = f(w)\; on\; \partial \Omega
and
\partial_t w - D_2\Delta w = 0,\;\; \partial_nw = g(u)\; on\; \partial \Omega

which means the PDEs are coupled by their boundary conditions.
The solution of u depends on w by a Dirichlet BC.
And w depends on u by a Neumann BC.

At the moment I first solve the equation for u, then compute g(u) and afterwards solve the equation for w, which allows me to compute f(w) for the next iteration, etc. ....

This works very nice for the stationary case, i.e. for \partial_t u = \partial_t w = 0
For the time dependent case however....the method is highly unstable.

I'm searching now for a method to solve both equations simultaneously, somehow - instead of one after another.

But I have no clue how to do this. If one searches for literature about this issue, it seems nobody has ever had the problem to solve PDEs coupled by boundary conditions.

Any hint? Any literature? Ideas?

BTW: I'm not working with any commercial software, but a own FEM software written in Fortran. Therefore please don't recommend commercial software.

Greetings,
Amiga-Freak

numerical instability appears if you do not use a time step within the stability region. What criterion have you used?
FMDenaro is offline   Reply With Quote

Old   March 1, 2016, 13:38
Default
  #3
New Member
 
Michael
Join Date: Mar 2016
Posts: 4
Rep Power: 2
Amiga-Freak is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
numerical instability appears if you do not use a time step within the stability region. What criterion have you used?
Well, there you have me....
I have not given much thought of this at all, so far
Amiga-Freak is offline   Reply With Quote

Old   March 1, 2016, 15:29
Default
  #4
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,692
Rep Power: 33
FMDenaro will become famous soon enoughFMDenaro will become famous soon enough
What about the functions for the bc.s ?
FMDenaro is offline   Reply With Quote

Old   March 1, 2016, 15:57
Default
  #5
New Member
 
Michael
Join Date: Mar 2016
Posts: 4
Rep Power: 2
Amiga-Freak is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
What about the functions for the bc.s ?
What do you like to know about them?
Amiga-Freak is offline   Reply With Quote

Old   March 1, 2016, 17:07
Default
  #6
New Member
 
Join Date: Apr 2014
Posts: 10
Rep Power: 4
pivalse is on a distinguished road
It seems that your instability (I am saying this without knowing anything about the way you implement this numerically...) is due to the fact that while the equation system is elliptical (everything happens instantaneously everywhere at each time-step), you introduce a time-lag between the first equation and the second.

Did you write the solver yourself?
In that case, writing it to solve the two equations at the same time should be straight-forward. A more elegant solution would be to write the solver as implicit. Compared to Navier-Stokes equation, your set is relatively simple and, above all, linear. This means that the implicit solver, even in 2D and 3D, remains quite simple.
The advantage is that the boundary conditions are imposed at the same of the solution.

If you did not write the code, you may "trick" the solver into solving both equations (they are the same) by using a vector that is twice as big {u | w} instead of each vector {u} and {w}
pivalse is offline   Reply With Quote

Old   March 1, 2016, 17:40
Default
  #7
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,692
Rep Power: 33
FMDenaro will become famous soon enoughFMDenaro will become famous soon enough
The set of equations is parabolic in the unsteady case and classical numerical integrations, such as the Crank-Nicolson method, can be used.

The only care is to provide BC.s that ensure that the time derivatives do not blow-up. For example, the Neumann BC.s for w can be critical, it is simple to see that by integration over the whole computational domain V:

d/dt Int [V] w dV = D2 Int[S] g(u) dS


Therefore, if g(u) is such that Int[S] g(u) dS > 0 then the integral of the w will blow-up in time
FMDenaro is offline   Reply With Quote

Old   March 1, 2016, 20:53
Default
  #8
New Member
 
Join Date: Apr 2014
Posts: 10
Rep Power: 4
pivalse is on a distinguished road
Thanks for the correction, FMDenaro. I guess I overlooked the derivative over time.
I agree, parabolic --> Crank-Nicholson

Please, disregard my previous comment
pivalse is offline   Reply With Quote

Old   March 2, 2016, 05:26
Default
  #9
New Member
 
Michael
Join Date: Mar 2016
Posts: 4
Rep Power: 2
Amiga-Freak is on a distinguished road
Quote:
Originally Posted by pivalse View Post
Did you write the solver yourself?
No, I'm a PhD student and the solver was written by other PhD students before me, several years ago.

Quote:
..., you may "trick" the solver into solving both equations (they are the same) by using a vector that is twice as big {u | w} instead of each vector {u} and {w}
That was actually the first idea, I had. But my problem is that I don't know how to impose the boundary conditions - because they are dependent on the solutions u and w.
Maybe It's easy and I'm just dumb, but I don't see how to do this.
Amiga-Freak is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
sliding mesh problem in CFX Saima CFX 45 September 22, 2015 10:53
Radiation in semi-transparent media with surface-to-surface model? mpeppels CFX 10 June 16, 2015 15:48
GETVAR Error in Multiband Monte Carlo Radiation Simulation with Directional Source silvan CFX 3 June 16, 2014 09:49
Moving mesh Niklas Wikstrom (Wikstrom) OpenFOAM Running, Solving & CFD 122 June 15, 2014 06:20
Radiation interface hinca CFX 15 January 26, 2014 18:11


All times are GMT -4. The time now is 08:02.