Mass production of hydrogen is a major issue for the coming decades particularly to decrease greenhouse gas production. The development of fourth-generation high-temperature nuclear reactors has led to renewed interest for hydrogen production. In France, the CEA is investigating new processes using nuclear reactors, such as the Westinghouse hybrid cycle. A recent study was devoted to electrical modeling of the hydrogen electrolyzer, which is the key unit of this process. In this electrochemical reactor, hydrogen is reduced at the cathode and SO2 is oxidized at the anode with the advantage of a very low voltage cell. This paper describes an improved model coupling the electrical and thermal phenomena with hydrodynamics in the electrolyzer, designed for a priori computational optimization of our future pilot cell. The hydrogen electrolyzer chosen here is a filter press design comprising a stack of identical cathode and anode compartments separated by a membrane. In a complex reactor of this type the main coupled physical phenomena involved are forced convection of the electrolyte flows, the plume of evolving hydrogen bubbles that modifies the local electrolyte conductivity, and all the irreversible processes that contribute to local overheating (Joule effect, overpotentials, etc.). The secondary current distribution was modeled with a commercial FEM code, Flux Expert (R), which was customized with specific finite interfacial elements capable of describing the potential discontinuity associated with the electrochemical overpotential. Since the finite element method is not capable of properly describing the complex two-phase flows in the cathode compartment, the Fluent (R) CFD code was used for thermohydraulic computations. In this way each physical phenomenon was modeled using the best numerical method. The coupling implements an iterative process in which each code computes the physical data it has to transmit to the other one: the two-phase thermohydraulic problem is solved by Fluent (R) using the Flux-Expert (R) current density and heat sources; the secondary distribution and heat losses are solved by Flux-Expert (R) using the Fluent (R) temperature field and flow velocities. A set of dedicated library routines was developed for process initiation, message passing, and synchronization of the two codes. The first results obtained with the two coupled commercial codes give realistic distributions for the electrical current density, gas fraction, and velocity in the electrolyzer. This approach allows us to optimize the design of a future experimental device.